一个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细节了，这里就不再赘述。

• 生信技能树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全部使用细节及流程，大家可以自行前往学习，如下：

但是这个两年前的系列笔记是基于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个文件组成的表达量矩阵：

每个样品都是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([e-mail], 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 ，但是并不在我之前背诵的基因列表里面。

但是，知识就是这样慢慢增加的！

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 [e-mail]

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。