scRNA-Seq | 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官网，点击对下方的Download Link链接

如果是第一次进入下载界面，需要填写一些基本信息，填写完后点击continue即可.

2. 根据需求下载Cell Ranger，可使用curl或者wget命令下载（在linux系统中运行黑框中的内容即可）;

wget -O cellranger-7.0.1.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-7.0.1.tar.gz?Expires=1669394099&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci03LjAuMS50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2NjkzOTQwOTl9fX1dfQ__&Signature=jE~gdVd5roKryS-0dA8R0RQazPC8ABzDNaUy40Opi9FrZqwUmkGCuNc5~Yn69BwGse81GTKIj5mMxZD0cnwu-A6oCbTFktns-WMGV6GaMa4sOVlHORcvwLjmF5lsb3tbFquwNNvxZe~NuBhUhFAZEBk7tRZjl4d~PqCYSDeRwmM~oR6mKUHRZ1yaKefxgy6R4CwQ~YaTcl7az5-Ucb0gmLkStYHQUuGmpy6M4LOds2kvcKcsDYMAw7Crghaja2Pn3rE-D71JyfHaBQ7EoM8gAijWumF1mYLK8yYCd1Tpb0JjWEzHR49eGX0iCQTZq~bfAcinhRLn90xZD2xfPlxhAA__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"

注：默认下载最新版的Cell Ranger，如果需要选择之前的版本可点击右下方的红框，选择想要的版本（如下图）；

3. 安装包下载完之后直接使用tar命令进行解压即可

tar -xzvf xxx.tar.gz

这样就完成Cell Ranger的安装啦

二、使用Cell Ranger进行单细胞转录组测序数据（scRNA-Seq）的定量

cellrange count 是 cellrange 中一个很重要的命令，用来对单细胞转录组数据进行基因组比对，细胞定量最终得到用于下游分析的单细胞表达矩阵（默认情况也会对表达矩阵进行聚类）。

在做定量之前，我们首先需要准备2组文件：原始fq文件以及物种的References（其中包括参考基因组序列、gtf文件以及star的索引文件）。

1. 原始fq文件

cellranger的输入文件格式是fq格式，并且文件的命名也是有要求，文件命名格式如下：

[Sample Name]S1_L00[Lane Number][Read Type]_001.fastq.gz

It is highly likely that these files were initially processed with bcl2fastq. Once you track the origin of the file, you will rename the files in the following format:

[Sample Name]S1_L00[Lane Number][Read Type]_001.fastq.gz

Where Read Type is one of:

• I1: Dual index i7 read (optional)
• I2: Dual index i5 read (optional)
• R1: Read 1
• R2: Read 2

链接：https://support.10xgenomics.com/single-cell-multiome-atac-gex/software/pipelines/latest/using/fastq-input#gex_rightname

如果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。
首先创建存放References的目录，这是我的一个习惯，也推荐大家在运行不同步骤的时候能够创建专门的文件，这样也便于文档管理。

mkdir refdata-cellranger-Arabidopsis-TAIR10

具体命令如下：

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目录，并且该目录下的文件有：

image

3. 定量

在所有文件都准备好了以后，就可以使用count对单细胞转录组数据进行定量啦。

具体命令如下（一般使用默认参数）：

cellranger count \
--id=sample_test  \
--transcriptome=/xx/ AT  \
--fastqs=/xxx/fastq_path  \
--localcores=8  \
--localmem=64

参数解释：

• id：样本名（唯一性）
• transcriptome：上一步创建的索引的目录名
• fastqs：下机数据的目录名
• localcores：内核
• localmem：内存

下面是我的脚本，和上面是同样的脚本格式~

image

成功运行之后会生成sample_test目录（脚本中id参数后面输入的内容），最终结果都保存在sample_test/outs中。

image

目录

• 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的下载、安装以及初步的使用流程，希望能够帮到大家啦~

三、结果解读

首先我们了解一下运行完Cell Ranger之后，在哪里可以看到生成的结果

还记得我们在运行Cell Ranger的时候有个参数--id吗？--id=XXX，这里的XXX就是最终生成的目录，该目录中保存了运行过程中所有的中间文件、日记文件以及最终的结果。如下图：

image

其中outs目录中即保存的最终结果，也是我们最后需要的。当然如果中间出现了报错，我们也可以通过查看日志文件，例如：_log，查看具体的报错原因，随后进行修改即可。

02 结果目录 outs

首先我们看一下outs目录下的文件结构，如下图：

image

这些结果中主要分成了两部分：
1. 集群中可以使用的结果（具体的内容可以参考上期文章“ 单细胞分析流程之Cell Ranger ” ）；
2. 网页版报告。

本期的重点是解读网页报告中的内容。

03 网页报告 web_summary.html

为了快速了解和方便的了解Cell Ranger定量之后的结果，我们首先会查看html文件，即 web_summary.html，了解初步情况。如下图：

image

可以看到该网页中主要分成了两部分：Summary和Analysis.

04 Summary

1. 异常结果警告

如果数据中存在异常情况，网页的上面会出现黄色的警告信息。找了一下之前遇到警告信息，如下图：

image

当遇到这种报错情况的时候我们不要慌，首先看一下是哪些值异常，对数据有无影响以及解决办法。在Detail部分会详细解释这个参数是什么，以及解决办法。例如上图中说到在运行Cell Ranger的时候可以调用--force-cells参数，这个参数的修改需要不断的尝试，所以也没有固定的值😭

当然如果这些报错信息并不影响结果，我们是可以用这个结果继续往后分析的～

2. 细胞和基因数的统计

随后就是查看这次分析中捕获到的细胞数以及基因数的情况，从这里就能大概知道数据的情况。

我也做过好多10X的数据，一般捕获的细胞数都是5,000-10,000，平均的基因数大概是1,200-15,00，大家可以看看自己的数据是否也在这些范围内。如果这些值都是在可接受的范围，那么就可以进入下一步的分析啦~

image

3. 细胞的选取

随后就是细胞的选取了（也是一个相当重要的图），帮助我们更加直观的筛选细胞（如下图）

image

先我们先来看一下上方的折线图怎么看：

Y轴是每个细胞中UMI的值，X轴是单个细胞的按照UMI大小的排序（降序），所以这个图中的曲线是下降的趋势。蓝色的线是选取的细胞（和2. 细胞和基因数的统计中的细胞数是一致的），灰色的线是背景。

正常的数据来说会有两个下降的趋势（如下图），第1个下降的趋势：区分完整细胞和背景物质（因为细胞和其他物质相比，真正细胞中会有更多的UMI，而其他物质可能没有或者由于一些污染能捕获到少量的转录本，所以会出现第一个下降的趋势）；第2个下降的趋势：区分细胞的质量，捕获率低或细胞破碎（这类细胞中基因数会很少，导致UMI数也少），而正常的细胞中UMI多且分布比较接近，所以质量好和不好的细胞在UMI上也会存在很大的差异，随后就出现了第2个下降趋势。

当数据出现了这两个下降趋势，且在蓝色区域的线条比较平稳时，也能说明我们的数据质量好~

image

4. 测序结果统计

继续往下走，下一部分是测序的信息，包括总的reads数目以及一些质控的指标，一般情况下Q30>90%表明质量是相当不错的。

image

当我们看数据的时候，如果遇到一些指标不太明白是什么意思，大家可以点击左上角的？，随后会列出下列指标的解释。

image

5. 比对结果统计

报告中除了会给出测序信息以外，也会给出与基因组的比对信息，主要包括Genome、Intergenic、Intronic、 Exonic、Transcriptome、Antisense to Gene（见下图）。

image

虽然测序和比对结果都是一些常规的质控信息，当我们数据一切正常的时候，看这些指标可能没有那么重要，但是一旦我们的数据比较奇怪的时候，例如发现检测到的细胞数还行，但是基因数特别少，这个时候测序和比对结果就相当重要了！之前遇到一个数据就是检测到的基因数特别少，然后聚类的时候就结果很差，后来就返回去看这些质控信息，惊奇的发现很多reads都是比对到了基因间区！

所以测序的reads根本就没有落在基因上，导致了最终每个细胞检测到的基因非常少，然后再去继续往下找原因。

。所以呀，还是得多看数据，从那以后，数据下来都会先看看这些质控信息是否正常，才会继续往后做（质控也是做科研非常重要的一步呀~）

6. 样本信息

最后一部分就是样本信息啦（如下图）~

image

这一部分就是在运行Cell Ranger时候的参数信息，例如样本名、Chemistry（运行Cell Ranger时候我们没有设置这个参数，那么就默认选择auto：自动配置，在报告中会给出具体的类型，这个就是3' V3版本）、Reference以及Reference路径等等。这些信息的给出方便后面查找信息。

image

05 Analysis

介绍完 Summary之后，下面就是Analysis.

1. 分群结果

image

左图：在TNSE中映射每个细胞UMI的值；右图：TSNE中分群的情况。

Cell Ranger做完定量之后呢，会默认拿已有的结果跑一下基本的分群，所以在看报告的时候我们也可以看一下这里的分群结果，心里大概有个数~

2. 基因差异表达分析

Cell Ranger除了做了分群以外，还找了每个群差异表达的基因，类似于Seurat中的 "FindAllMarkers"。

image

这里比较好的是，上面Graph-based如果选择K=2，那么这里差异基因列表也会随之变动。所以如果觉得Cell Ranger的分群结果已经很符合自己的预期了，完全可以就用这个结果了，而且还可以自己选择分群的个数（直接网页挑选，人性化呀）

3. 饱和度评估

对 reads 抽样，计算不同抽样条件下检测到的转录本数量占检测到的所有转录本的比例（测序饱和度)，如下图：

image

曲线末端接近平滑状态说明测序达到饱和，因为继续增加测序量，检测到的转录本也不会有特别大的变化

对 reads 抽样，计算不同抽样条件下检测基因数目的分布，如下图：

image

同样地，曲线末端接近平滑状态说明测序达到饱和，因为继续增加测序量，每个细胞检测到的基因数也不会有特别大的变化

下游 barcodes.tsv.gz/features.tsv.gz/matrix.mtx.gz

cellranger count输出结果中的outs.文件夹有几个是非常重要的信息，我们今天只关注于filtered_feature_bc_matrix文件夹下的内容和possorted_genome_bam.bam文件。

image

一般来说，我们下游的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是按字母表顺序的，所以（谁知道它是不是最后随意用字母表顺序输出的呢？）

需了解

将Reads转为 Count Matrix。

1. 参考基因组及注释

目前，大多数scRNA-seq是使用人类或小鼠的组织、器官或细胞培养物进行的。常用的就是UCSC（hg19、hg38、mm10等），和GRC（GRCh37、GRCh38、GRCm38）。二者在主要染色体上是对等的（如hg38的chr1 = GRCh38的chr1），但在一些小的位点上会有细微差异。

基因组注释过程包括定义基因组的转录区域，明确exon和intron，将其分成protein coding, non-coding等。举个栗子，假设我们有一个基因，包含5个转录本组成的基因。其中3个编码（红色）和2个非编码（蓝色）。

在实际操作中，我们通常可以下载GTF或GFF3格式的文件进行注释。每个基因都含有一个ID，而这个ID是唯一的。

Note! 这里也提醒大家在实际操作中，尽量使用ID进行分析操作，而不要使用symbol，当然在展示结果的时候你需要转换回symbol方便阅读。

我们常用的人类和小鼠基因组注释包括RefSeq, ENSEMBL和GENCODE，实际应用中选择最新的版本就可以了，会有更多的已知基因。

2. Full-length scRNA-seq的处理

处理方法与bulk RNA-seq类似。

Full length scRNA-seq的raw data的处理通常分两步进行：比对(read alignment)和计数 (read counting)。常用软件：STAR和hisat2。normalization方法：推荐使用TPM。

Droplet-based scRNA-seq的比对和定量
1️⃣ 首先我们要搞清楚scRNA-seq都有哪些产物。

• cDNA片段 (识别转录本)；
• Cell barcode (CB，识别细胞)；
• Unique Molecular Identifier (UMI，减小PCR扩增带来的bias)。

2️⃣ 典型的scRNA-seq的workflow包括以下几个步骤：👇

• 将cDNAmapping到reference上；
• 计算基因reads；
• 计算细胞reads（用到cell barcode）；
• 计算的RNA数量（UMI去重）。

3. 具体步骤

Read Mapping
处理10x Genomics Chromium scRNAseq数据，我们通常要用到Cell Ranger。这里只介绍一下外显子(exon)的定义，即reads比对到外显子的 50% 以上，就可以定义为外显子。

Cell Ranger Reference
在选择Reference的时候，大家可以按以下table进行选择。👇

Cell Ranger Reference   Species Assembly/Annotation Genes before filtering  Genes after filtering
2020-A  human   GRCh38/GENCODE v32  60668   36601
2020-A  mouse   mm10/GENCODE vM23   55421   32285
3.0.0   human   GRCh38/Ensembl 93   58395   33538
3.0.0   human   hg19/Ensembl 87 57905   32738
3.0.0   mouse   mm10/Ensembl 93 54232   31053
2.1.0   mouse   mm10/Ensembl 84 47729   28692
1.2.0   human   GRCh38/Ensembl 84   60675   33694
1.2.0   human   hg19/Ensembl 82 57905   32738
1.2.0   mouse   mm10/Ensembl 84 47729   27998

UMI计数
这里我们注意一下两点：

• 如果两组或更多的reads具有相同的barcode和UMI，但基因注释却不相同，那么reads最多的基因注释用于UMI计数，舍弃其他组。

• 我们再极端一点，如果两组reads一样的话，这个时候我们可能需要舍弃所有组，因为这个时候基因注释已经不准确了。

细胞过滤
1️⃣ 未经过滤的raw data, feature-barcode matrix会包含很多空的 droplets，在矩阵中并不是0，因为会有来自破碎细胞的RNA。所以，这种数据我们需要进行过滤，而后再进行分析。

2️⃣ 我们通常需要使用Cell Ranger 2.2和Cell Ranger 3.0进行过滤。

3️⃣ 举个栗子

肿瘤样本通常包含大型肿瘤细胞与少量的肿瘤浸润淋巴细胞（TIL），如果你对TIL特别感兴趣，那这个时候可能就要用到EmptyDrops的方法来进行过滤了。

其他方法
当你使用Cell Ranger时，你可能会觉得它不够快，这里我们介绍两个速度快、准确度高的方法，STARsolo和Alevin，这里不做具体介绍了，推荐大家选择STARsolo。

参考：
https://www.jianshu.com/p/3f01016b5302
https://mp.weixin.qq.com/s/iB1o4BhCElU5IExOtXbtmw
https://mp.weixin.qq.com/s/9QA20YQEe1PPkDarB964HQ