单细胞分析流程之Cell Ranger
相信做单细胞的小伙伴对Cell Ranger这个软件都不陌生，我们今天就来了解一下Cell Ranger的安装和使用方法。
Cell Ranger是10X Genomics为单细胞分析专门打造的分析软件，直接对10X的下机数据进行基因组比对、定量、生成单细胞矩阵、聚类以及其他的分析等。所以Cell Ranger能做的分析有很多，我们今天主要学一下Cell Ranger的安装以及对单细胞RNA-Seq数据的定量。
Cell Ranger的官网：https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger
1 Cell Ranger的下载与安装
1. 首先进入Cell Ranger官网，点击对下方的Download Link链接;
如果是第一次进入下载界面，需要填写一些基本信息，填写完后点击continue即可（如下：）
2. 根据需求下载Cell Ranger，可使用curl或者wget命令下载（在linux系统中运行黑框中的内容即可）;
注：默认下载最新版的Cell Ranger，如果需要选择之前的版本可点击右下方的红框，选择想要的版本（如下图）；
3. 安装包下载完之后直接使用tar命令进行解压即可。
tar -xzvf xxx.tar.gz
这样就完成Cell Ranger的安装啦
2 使用Cell Ranger进行单细胞转录组测序数据（scRNA-Seq）的定量
因为我是做单细胞转录组方向的，所以下面介绍一下常用的Cell Ranger命令---cellrange count。
count也是cellrange中一个很重要的命令，用来对单细胞转录组数据进行基因组比对，细胞定量最终得到用后下游分析的单细胞表达矩阵（默认情况也会对表达矩阵进行聚类）。
在做定量之前，我们首先需要准备2组文件：原始fq文件以及物种的References（其中包括参考基因组序列、gtf文件以及star的索引文件）。
1. 原始fq文件
cellranger的输入文件格式是fq格式，并且文件的命名也是有要求，文件命名格式如下：
**[Sample Name]**S1_L00**[Lane Number]****[Read Type]**_001.fastq.gz
如果fq的文件名格式不对，在运行的过程中会出现错误，所以最开始需要确定文件名的格式以及进行修改。习惯是重新创建一个目录并且用软连接将原始文件链接到新的目录中，这样做的好处是首先不会改变原始文件的名字（害怕修改了文件名后有些文件没有同步，导致最后找不到具体的文件），其实也不会占用很多存储（毕竟我们还要在夹缝中生存），下面就是我使用的风格：
2. 物种的References
第二个需要准备的文件就是物种的References。
好消息就是Cell Ranger官网已经为我们提供了人和小鼠的References，如果大家的样本是人或者小鼠的某些细胞可以直接去Cell Ranger官网进行下载。
下载流程和Cell Ranger软件下载流程一致，其中也是有很多版本的References可供大家选择，下载后解压就可用了；
下载网页：
https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest?
那么问题来了，如果我研究的是其他物种，那怎么构建这个References？
cellranger的mkref就是这么一个功能，可用对其他的物种构建cellranger需要的References格式，只需要准备物种的参考基因组序列和gtf注释文件就可以直接运行。
这里就以拟南芥为例子构建References。
mkdir refdata-cellranger-Arabidopsis-TAIR10 #首先创建存放References的目录，这是我的一个习惯，也推荐大家在运行不同步骤的时候能够创建专门的文件，这样也便于文档管理。
具体命令如下：
cellranger mkref \
--genome=TAIR10 \
--nthreads=10 \
--fasta=TAIR10.fa \
--genes=TAIR10_GFF3_genes.miRBase20.gtf
--genome：生成索引的目录
--fasta：基因组序列
--genes：基因注释文件（gtf格式）
运行完上面的命令就构建完索引啦~
这里还要推荐一个运行脚本的命令，希望能够对大家有帮助~
我们可以使用vi编辑器，将上面的内容存放在一个shell脚本中，然后使用后面运行shell脚本，这样后台在运行的同时，我们仍然可以在当前界面进行其他操作，并且网络不稳定的时候也不会影响我们的运行，所以非常推荐。（脚本名：index_test.sh）~
投后台的命令是：
nohup sh index_test.sh >index_test.sh.o 2>index_test.sh.e &
这样的话中间的输出文件会保存在index_test.sh.o，如果脚本报错就会保存在index_test.sh.e中。我们可以通过查看这两个文件了解运行的进展。可以通过使用jobs命令查看后台运行的命令是不是还在。
References构建完后就会生成TAIR10目录，并且该目录下的文件有：
3. 定量
在所有文件都准备好了以后，就可以使用count对单细胞转录组数据进行定量啦。
具体命令如下（一般使用默认参数）：
cellranger count \
--id=sample_test \
--transcriptome=/xx/ AT \
--fastqs=/xxx/fastq_path \
--localcores=8 \
--localmem=64
参数解释：
id：样本名（唯一性）
transcriptome：上一步创建的索引的目录名
fastqs：下机数据的目录名
localcores：内核
localmem：内存
下面是我的脚本，和上面是同样的脚本格式~
成功运行之后会生成sample_test目录（脚本中id参数后面输入的内容），最终结果都保存在sample_test/outs中。
目录
analysis：cellranger聚类的结果
filtered_feature_bc_matrix：过滤后的单细胞表达矩阵（后续可以对接到seurat中）
raw_feature_bc_matrix：过滤前的单细胞表达数据（一般不怎么使用）
文件
possorted_genome_bam.bam：单细胞比对的bam文件，其中包含了每个reads的信息
web_summary.html：报告网页（单细胞定量后的报告，包括检测到的细胞数、基因数、UMI、分群等等）
END Cell Ranger
以上就是cellranger的下载、安装以及初步的使用流程，希望能够帮到大家啦~
结果解读
01
首先我们了解一下运行完Cell Ranger之后，在哪里可以看到生成的结果。
还记得我们在运行Cell Ranger的时候有个参数--id吗？--id=XXX，这里的XXX就是最终生成的目录，该目录中保存了运行过程中所有的中间文件、日记文件以及最终的结果。如下图：
其中outs目录中即保存的最终结果，也是我们最后需要的。当然如果中间出现了报错，我们也可以通过查看日志文件，例如：_log，查看具体的报错原因，随后进行修改即可。
02
结果目录"outs"
首先我们看一下outs目录下的文件结构，如下图：
这些结果中主要分成了两部分：1. 集群中可以使用的结果（具体的内容可以参考上期文章“单细胞分析流程之Cell Ranger”）；2. 网页版报告。
本期的重点是解读网页报告中的内容。
03
网页报告"web_summary.html"**
为了快速了解和方便的了解Cell Ranger定量之后的结果，我们首先会查看html文件，即web_summary.html，了解初步情况。如下图：
可以看到该网页中主要分成了两部分：Summary和Analysis.
04
"Summary"**
1. 异常结果警告
如果数据中存在异常情况，网页的上面会出现黄色的警告信息。找了一下之前遇到警告信息，如下图：
当遇到这种报错情况的时候我们不要慌，首先看一下是哪些值异常，对数据有无影响以及解决办法。在Detail部分会详细解释这个参数是什么，以及解决办法。例如上图中说到在运行Cell Ranger的时候可以调用--force-cells参数，这个参数的修改需要不断的尝试，所以也没有固定的值😭
当然如果这些报错信息并不影响结果，我们是可以用这个结果继续往后分析的～
2. 细胞和基因数的统计
随后就是查看这次分析中捕获到的细胞数以及基因数的情况，从这里就能大概知道数据的情况。
我也做过好多10X的数据，一般捕获的细胞数都是5,000-10,000，平均的基因数大概是1,200-15,00，大家可以看看自己的数据是否也在这些范围内。如果这些值都是在可接受的范围，那么就可以进入下一步的分析啦~
3. 细胞的选取
随后就是细胞的选取了（也是一个相当重要的图），帮助我们更加直观的筛选细胞（如下图）
先我们先来看一下上方的折线图怎么看：
Y轴是每个细胞中UMI的值，X轴是单个细胞的按照UMI大小的排序（降序），所以这个图中的曲线是下降的趋势。蓝色的线是选取的细胞（和**2\. 细胞和基因数的统计**中的细胞数是一致的），灰色的线是背景。
正常的数据来说会有两个下降的趋势（如下图），第1个下降的趋势：区分完整细胞和背景物质（因为细胞和其他物质相比，真正细胞中会有更多的UMI，而其他物质可能没有或者由于一些污染能捕获到少量的转录本，所以会出现第一个下降的趋势）；第2个下降的趋势：区分细胞的质量，捕获率低或细胞破碎（这类细胞中基因数会很少，导致UMI数也少），而正常的细胞中UMI多且分布比较接近，所以质量好和不好的细胞在UMI上也会存在很大的差异，随后就出现了第2个下降趋势。
当数据出现了这两个下降趋势，且在蓝色区域的线条比较平稳时，也能说明我们的数据质量好~
4. 测序结果统计
继续往下走，下一部分是测序的信息，包括总的reads数目以及一些质控的指标，一般情况下Q30>90%表明质量是相当不错的。
当我们看数据的时候，如果遇到一些指标不太明白是什么意思，大家可以点击左上角的？，随后会列出下列指标的解释。
5. 比对结果统计
报告中除了会给出测序信息以外，也会给出与基因组的比对信息，主要包括Genome、Intergenic、Intronic、 Exonic、Transcriptome、Antisense to Gene（见下图）。
虽然测序和比对结果都是一些常规的质控信息，当我们数据一切正常的时候，看这些指标可能没有那么重要，但是一旦我们的数据比较奇怪的时候，例如发现检测到的细胞数还行，但是基因数特别少，这个时候测序和比对结果就相当重要了！之前遇到一个数据就是检测到的基因数特别少，然后聚类的时候就结果很差，后来就返回去看这些质控信息，惊奇的发现很多reads都是比对到了基因间区！
所以测序的reads根本就没有落在基因上，导致了最终每个细胞检测到的基因非常少，然后再去继续往下找原因。
。所以呀，还是得多看数据，从那以后，数据下来都会先看看这些质控信息是否正常，才会继续往后做（质控也是做科研非常重要的一步呀~）
6. 样本信息
最后一部分就是样本信息啦（如下图）~
这一部分就是在运行Cell Ranger时候的参数信息，例如样本名、Chemistry（运行Cell Ranger时候我们没有设置这个参数，那么就默认选择auto：自动配置，在报告中会给出具体的类型，这个就是3' V3版本）、Reference以及Reference路径等等。这些信息的给出方便后面查找信息。
05
"Analysis"**
****介绍完Summary之后，下面就是Analysis.
1. 分群结果
左图：在TNSE中映射每个细胞UMI的值；右图：TSNE中分群的情况。
Cell Ranger做完定量之后呢，会默认拿已有的结果跑一下基本的分群，所以在看报告的时候我们也可以看一下这里的分群结果，心里大概有个数~
2. 基因差异表达分析
Cell Ranger除了做了分群以外，还找了每个群差异表达的基因，类似于Seurat中的 "FindAllMarkers"。
这里比较好的是，上面Graph-based如果选择K=2，那么这里差异基因列表也会随之变动。所以如果觉得Cell Ranger的分群结果已经很符合自己的预期了，完全可以就用这个结果了，而且还可以自己选择分群的个数（直接网页挑选，人性化呀）
3. 饱和度评估
对 reads 抽样，计算不同抽样条件下检测到的转录本数量占检测到的所有转录本的比例（测序饱和度)，如下图：
曲线末端接近平滑状态说明测序达到饱和，因为继续增加测序量，检测到的转录本也不会有特别大的变化
对 reads 抽样，计算不同抽样条件下检测基因数目的分布，如下图：
同样地，曲线末端接近平滑状态说明测序达到饱和，因为继续增加测序量，每个细胞检测到的基因数也不会有特别大的变化
下游barcodes.tsv.gz/features.tsv.gz/matrix.mtx.gz
cellranger count输出结果中的outs.文件夹有几个是非常重要的信息，我们今天只关注于filtered_feature_bc_matrix文件夹下的内容和possorted_genome_bam.bam文件。
一般来说，我们下游的Seurat分析的输入文件会选择filtered_feature_bc_matrix中的文件，而不选择raw_feature_bc_matrix下的文件，前者是经过过滤的，去掉了低质量的信息。进入filtered_feature_bc_matrix文件夹会发现它下面包含3个文件：分别是barcodes.tsv.gz、features.tsv.gz和matrix.mtx.gz。
barcodes.tsv.gz
AAACCCAAGAGATGCC-1
AAACCCAAGGTCGTAG-1
AAACCCACATCAGTCA-1
AAACCCAGTTTCCCAC-1
AAACCCATCCAAACCA-1
AAACCCATCCCTCTAG-1
AAACGAAAGCTGGTGA-1
AAACGAACAGACACAG-1
AAACGAAGTGAGATAT-1
这个文件当中记载了每个细胞的barcode信息。
features.tsv.gz
ENSMUSG00000051951 Xkr4 Gene Expression
ENSMUSG00000089699 Gm1992 Gene Expression
ENSMUSG00000102331 Gm19938 Gene Expression
ENSMUSG00000102343 Gm37381 Gene Expression
ENSMUSG00000025900 Rp1 Gene Expression
ENSMUSG00000025902 Sox17 Gene Expression
ENSMUSG00000104238 Gm37587 Gene Expression
ENSMUSG00000104328 Gm37323 Gene Expression
这个文件记载了小鼠基因注释文件中包含的基因id与symbol信息，注意，这个文件的来源是小鼠基因组的注释文件。
matrix.mtx.gz
%%MatrixMarket matrix coordinate integer general
%metadata_json: {"software_version": "cellranger-6.0.1", "format_version": 2}
32285 5741 11436472
1 1 4
2 1 1
22 1 1
24 1 8
31 1 1
41 1 1
43 1 1
这个文件主体部分包含三列，第一列为基因，即这个基因在前面features.tsv.gz中的位置；第二列为细胞，即这个细胞对应于barcodes.tsv.gz中的barcodes信息；最后一列代表在这个细胞中检测到的这个基因的reads数。举个例子来说：
例如第一行：1 1 4，就表示barcode为AAACCCAAGAGATGCC-1的细胞中检测到的Xkr4基因的reads数为4。
细心的朋友会发现在前面还有一行：32285 5741 11436472 ，这一行实际上就是一个汇总信息，例如有32285个基因，5741个细胞，11436472个非零数值。而最前面不过是指明软件的相关信息罢了。
思考
实际上在我们进行数据分析时，都觉得这3个文件一个不可少，但实际上真的是这样吗？
-
features.tsv.gz
前面已经说到，这个文件实际上是来源于小鼠基因组的注释文件，所以理论上只要你在使用cellranger count时用的基因组注释文件是一样的，这个文件是不会变的，你可以进入Cell Ranger推荐的参考基因组看是否是这样。
cd cellranger/reference/refdata-gex-mm10-2020-A/genes
#这个文件夹下面你会看到一个小鼠基因组的gtf注释文件，名称应该为genes.gtf
cat genes.gtf | awk '$15=="gene_name"{print$10"\t"$16}' | less -S
#看看这样提取的基因id和name是否和features.tsv.gz一样
"ENSMUSG00000051951"; "Xkr4";
"ENSMUSG00000089699"; "Gm1992";
"ENSMUSG00000102331"; "Gm19938";
"ENSMUSG00000102343"; "Gm37381";
"ENSMUSG00000025900"; "Rp1";
"ENSMUSG00000025902"; "Sox17";
"ENSMUSG00000104238"; "Gm37587";
你会发现，顺序和内容竟然和features.tsv.gz一样的，所以看起来似乎features.tsv.gz也不是那么不可或缺，咱也可以自己做，或者说可以通用。
-
matrix.mtx.gz
这个文件，毫无疑问，是必不可少的，可以说花那么多钱做个single cell RNA sequencing就是为了这个文件。。
-
barcodes.tsv.gz
光听这个文件的内容，感觉这个文件很重要，像某个地区居民的名单一样，丢了岂不麻烦大了？但实际上仔细想想，它真的重要到我们不能丢吗？
我们说，matrix.mtx.gz里面实际上已经包含了单个细胞、单个基因的表达信息了，这是cellranger count已经返给我们的信息，举个形象的例子，小孩子在出生时，当地户籍部门记录了这个小孩的性别信息，当然还有他的名字。但是一年后，这个小朋友改名字了，但是他的性别变了吗？并没有！所以实际上这个barcodes.tsv.gz文件如果我们改了，只不过是给每个细胞新起了一个名字，本身并不会造成细胞RNA信息的变化和混乱。
说到这里，不得不提到possorted_genome_bam.bam文件，这个文件里面实际上包含了每个细胞的barcode信息，就在其中以CB开头的那个字段里。
samtools view possorted_genome_bam.bam | less -S
#部分信息如下
CB:Z:ATTCTTGTCTCCTGTG-1
CB:Z:GTGCTGGTCACTCGAA-1
CB:Z:GCATGATAGCCGGATA-1
CB:Z:GCACGTGGTTGCCTAA-1
你可以把这部分信息提取出来，重复内容合并，然后以任意顺序作为barcodes.tsv.gz就可以进行Seurat分析了。哦对了，得某位大佬指点，cellranger count输出的barcodes.tsv.gz是按字母表顺序的，所以（谁知道它是不是最后随意用字母表顺序输出的呢？）
