一个10x单细胞转录组项目从fastq到细胞亚群
最近在安排学徒单细胞分享的时候，有一个学徒提到了GSE168522这个数据集，是很标准的6个10x单细胞转录组样品，如下所示：
GSM5145401 Sample 16_Normal-1
GSM5145402 Sample 17_Normal-2
GSM5145403 Sample 23_Early AD-1
GSM5145404 Sample 28_Early AD-2
GSM5145405 Sample 24_Late AD-1
GSM5145406 Sample 27_Late AD-2
但是作者提供的并不是纯粹的表达量矩阵，而是被拆分成为了不同单细胞亚群的：
GSE168522_B_cell_data.txt.gz 318.7 Kb
GSE168522_CD4_T_cell_data.txt.gz 317.3 Kb
GSE168522_CD8_T_cell_data.txt.gz 320.2 Kb
GSE168522_Gene_expression_of_six_samples_data.txt.gz 277.5 Kb
GSE168522_HSC_data.txt.gz 271.0 Kb
GSE168522_Monocytes_data.txt.gz 322.1 Kb
GSE168522_NK_cell_data.txt.gz 319.8 Kb
这样的文件很明显没办法给我们跑单细胞转录组流程，看了看原文：《Single-cell RNA sequencing reveals B cell–related molecular biomarkers for Alzheimer’s disease》，其实在《单细胞天地》公众号有它的介绍：单细胞测序揭示阿尔兹海默症的B细胞相关标志物
可以看到原文提到的每个10x的单细胞转录组样品细胞数量蛮合理（5到8千），如下所示：
作者的降维聚类分群也是超级简单，就是第一层次而已，免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群：
参考前面的例子：人人都能学会的单细胞聚类分群注释 ，也是简简单单看比例变化：
既然作者没有提供可以直接分析的表达量矩阵，那么我们就需要从其提供的fastq测序数据开始啦。接下来我们就演示这个过程：
构建Linux服务器操作环境
在你的服务器空白用户里面安装自己的conda，每个用户独立操作，安装方法代码如下：
# 首先下载文件，20M/S的话需要几秒钟即可
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件，记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh
# 安装成功后需要更新系统环境变量文件
source ~/.bashrc
接下来 使用conda安装aspera，代码如下：
conda create -n download
conda activate download
conda install -y -c hcc aspera-cli
conda install -y -c bioconda sra-tools
which ascp
## 一定要搞清楚你的软件被conda安装在哪
ls -lh ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh
我们已经多次介绍过conda细节了，这里就不再赘述。
conda管理生信软件一文就够 生信技能树B站软件安装视频 https://www.bilibili.com/video/av28836717
接下来配置cellranger的数据库文件和软件：
使用aspera从ENA数据库下载fastq文件
首先需要下面的地址信息，自己去ENA数据库按图索骥即可下载，构建成为fq.txt 文件 ：
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/009/SRR13911909/SRR13911909_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/009/SRR13911909/SRR13911909_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/010/SRR13911910/SRR13911910_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/010/SRR13911910/SRR13911910_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/011/SRR13911911/SRR13911911_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/011/SRR13911911/SRR13911911_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/012/SRR13911912/SRR13911912_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/012/SRR13911912/SRR13911912_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/013/SRR13911913/SRR13911913_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/013/SRR13911913/SRR13911913_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/014/SRR13911914/SRR13911914_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR139/014/SRR13911914/SRR13911914_2.fastq.gz
然后写脚本批量下载；
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
这个脚本就会调用我们前面使用conda安装好的aspera，去读取fq.txt文件里面的路径，一个个下载，可以看到如下所示从晚上八点半到凌晨一点，我们的6个10x单细胞转录组样品的fastq数据文件就全部下载成功咯 。
~$ ls -lh *gz|cut -d" " -f 5-
7.5G 1月 29 20:25 SRR13911909_1.fastq.gz
31G 1月 29 20:53 SRR13911909_2.fastq.gz
7.5G 1月 29 20:58 SRR13911910_1.fastq.gz
31G 1月 29 21:20 SRR13911910_2.fastq.gz
7.3G 1月 29 21:28 SRR13911911_1.fastq.gz
31G 1月 29 21:53 SRR13911911_2.fastq.gz
7.2G 1月 29 21:58 SRR13911912_1.fastq.gz
30G 1月 29 22:58 SRR13911912_2.fastq.gz
7.5G 1月 29 23:15 SRR13911913_1.fastq.gz
31G 1月 30 00:12 SRR13911913_2.fastq.gz
7.2G 1月 30 00:21 SRR13911914_1.fastq.gz
30G 1月 30 00:49 SRR13911914_2.fastq.gz
运行cellranger
目前单细胞转录组以10X公司为主流，我们也是在单细胞天地公众号详细介绍了cellranger全部使用细节及流程，大家可以自行前往学习，如下：
单细胞实战(一)数据下载 单细胞实战(二) cell ranger使用前注意事项 单细胞实战(三) Cell Ranger使用初探 单细胞实战(四) Cell Ranger流程概览 单细胞实战(五) 理解cellranger count的结果
但是这个两年前的系列笔记是基于V2,V3版本的cellranger，在2020的7月我看到了其更新到了V4，也里面写了一个总结，见：cellranger更新到4啦（全新使用教程）
使用方法超级简单， 就是 构建一个简单的脚本即可，代码如下所示：
$ cat run-cellranger.sh
bin=/home/data/bylin/cellranger-6.1.2/bin/cellranger
db=/home/data/bylin/refdata-gex-GRCh38-2020-A
ls $bin; ls $db
fq_dir=/home/x9/raw
$bin count --id=$1 \
--localcores=4 \
--transcriptome=$db \
--fastqs=$fq_dir \
--sample=$1 \
--expect-cells=5000
接下来就是使用这个脚本去批量处理我们前面的6个10x单细胞转录组数据啦。
前面的 SRR 开头的fastq文件名字，我们统一修改一下：
Early_AD-1_S1_L001_R1_001.fastq.gz
Early_AD-1_S1_L001_R2_001.fastq.gz
Early_AD-2_S1_L001_R1_001.fastq.gz
Early_AD-2_S1_L001_R2_001.fastq.gz
Late_AD-1_S1_L001_R1_001.fastq.gz
Late_AD-1_S1_L001_R2_001.fastq.gz
Late_AD-2_S1_L001_R1_001.fastq.gz
Late_AD-2_S1_L001_R2_001.fastq.gz
Normal-1_Homo_S1_L001_R1_001.fastq.gz
Normal-1_Homo_S1_L001_R2_001.fastq.gz
Normal-2_Homo_S1_L001_R1_001.fastq.gz
Normal-2_Homo_S1_L001_R2_001.fastq.gz
可以看到比较规则的文件名，每个样品都是 _S1_L001_R1_001.fastq.gz 以及 _S1_L001_R2_001.fastq.gz的结尾内容，前面的前缀就是每个样品不一样，也是后面我们跑流程的时候用来区分不同样品的前缀。
所以需要构建一个文本文件id.txt，内容如下所示：
x9@bio1-2288H-V5:/home/data/x9$ cat id.txt
Early_AD-1
Early_AD-2
Late_AD-1
Late_AD-2
Normal-1_Homo
Normal-2_Homo
它只有6行简单的内容，就是我们的fastq文件的前缀。
一个简单的脚本就可以处理全部的6个10x单细胞转录组数据文件：
cat id.txt |while read id;do (nohup bash run-cellranger.sh $id 1>log-$id.txt 2>&1 & );done
可以看到：
drwxrwxr-x 4 x9 x9 4.0K 1月 30 02:55 Early_AD-1
drwxrwxr-x 4 x9 x9 4.0K 1月 30 02:54 Early_AD-2
-rw-rw-r-- 1 x9 x9 140 1月 29 20:41 id
-rw-rw-r-- 1 x9 x9 70 1月 29 20:42 id.txt
drwxrwxr-x 4 x9 x9 4.0K 1月 30 03:00 Late_AD-1
drwxrwxr-x 4 x9 x9 4.0K 1月 30 02:51 Late_AD-2
-rw-rw-r-- 1 x9 x9 95K 1月 30 02:55 log-Early_AD-1.txt
-rw-rw-r-- 1 x9 x9 96K 1月 30 02:54 log-Early_AD-2.txt
-rw-rw-r-- 1 x9 x9 96K 1月 30 03:00 log-Late_AD-1.txt
-rw-rw-r-- 1 x9 x9 95K 1月 30 02:51 log-Late_AD-2.txt
-rw-rw-r-- 1 x9 x9 88K 1月 30 02:44 log-Normal-1_Homo.txt
-rw-rw-r-- 1 x9 x9 9.6K 1月 29 21:09 log-Normal-2_Homo.txt
drwxrwxr-x 4 x9 x9 4.0K 1月 30 02:44 Normal-1_Homo
drwxrwxr-x 5 x9 x9 4.0K 1月 29 21:09 Normal-2_Homo
其中 4.0K 的是文件夹，6个样品输出6个文件夹，然后 95K 的是txt日志，每次运行 cellranger的日志都保留好，以后可以检查。我这里是以 x9 这个用户在运行命令。
可以看到， 基本上是在凌晨三点结束的任务。接下来需要把全部运行成功的文件夹里面的必须内容整理一下，默认多个项目都是同一个文件夹下面运行的，可以看到每个样品对应的文件夹里面的格式都是类似的
其中最重要的 filtered_feature_bc_matrix 文件夹里面的内容， 以及 web_summary.html 这个报表：
仍然是需要写一个代码：
pro=tmp
pro_folder=$HOME/$pro
mkdir -p $pro_folder
ls -lh $pro_folder
mkdir -p $pro_folder/html
ls */outs/web_summary.html |while read id;do (cp $id $pro_folder/html/${id%%/*}.html );done
mkdir -p $pro_folder/matrix
ls -d */outs/filtered_feature_bc_matrix |while read id;do (cp -r $id $pro_folder/matrix/${id%%/*} );done
这样我们的报表和表达量矩阵就都整理完毕：
上游流程就到此为止啦，接下来就是读取每个样品的表达量矩阵去R语言里面跑seurat流程，每个样品都是3个文件组成的表达量矩阵：
有了这样的表达量矩阵，就很容易读取它，走标准降维聚类分群流程啦。
首先是批量读取6个10x的单细胞转录组样品
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
###### step1:导入数据 ######
library(data.table)
dir='matrix/'
samples=list.files( dir )
samples
library(data.table)
sceList = lapply(samples,function(pro){
# pro=samples[1]
print(pro)
sce=CreateSeuratObject( Read10X(file.path(dir,pro)),
project = pro,
min.cells = 5,
min.features = 300 )
return(sce)
})
names(sceList)
samples
sce.all=merge(x=sceList[[1]],
y=sceList[ -1 ],
add.cell.ids = samples)
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)
然后需要进行常规的质量控制，以及 harmony整合：
sce = sce.all
sce <- NormalizeData(sce,
normalization.method = "LogNormalize",
scale.factor = 1e4)
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce)
sce <- RunPCA(sce, features = VariableFeatures(object = sce))
library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj, dims = 1:15,
reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T )
sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
dims = 1:15)
sce.all=sce
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
sce.all=FindClusters(sce.all, #graph.name = "CCA_snn",
resolution = res, algorithm = 1)
}
接下来就是可视化我们给大家的标记基因，并且合理的命名：
sce.all
sce=sce.all
library(ggplot2)
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A',
'CD19', 'CD79A', 'MS4A1' ,
'IGHG1', 'MZB1', 'SDC1',
'CD68', 'CD163', 'CD14',
'TPSAB1' , 'TPSB2', # mast cells,
'RCVRN','FPR1' , 'ITGAM' ,
'C1QA', 'C1QB', # mac
'S100A9', 'S100A8', 'MMP19',# monocyte
'FCGR3A',
'LAMP3', 'IDO1','IDO2',## DC3
'CD1E','CD1C', # DC2
'KLRB1','NCR1', # NK
'FGF7','MME', 'ACTA2', ## fibo
'DCN', 'LUM', 'GSN' , ## mouse PDAC fibo
'MKI67' , 'TOP2A',
'PECAM1', 'VWF', ## endo
'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' )
library(stringr)
genes_to_check=str_to_upper(genes_to_check)
genes_to_check
p <- DotPlot(sce.all, features = unique(genes_to_check),
assay='RNA' ) + coord_flip()
p
ggsave('check_last_markers.pdf',height = 11,width = 11)
如下所示：
可以看到免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类。
我们给名字后如下所示：
其中，18，20，21我第一次没办法给出名字，需要后续慢慢探索，但是总体上跟原文的结果是一致的。我没有给名字的可能是，它的标记基因是 NCOR2, NKX3-1, HLX, PRNP ，但是并不在我之前背诵的基因列表里面。
但是，知识就是这样慢慢增加的！
写在文末
我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。