利用cellranger分析单细胞数据

背景

当前的单细胞测序主要采用 illumina 测序平台进行测序，一般为双末端测序，测序完成之后首先需要对 illumina 测序数据进行质控过滤，过滤条件与其他分析类似。需要注意的是，虽然单细胞测序也是双末端测序，但是 reads1 中通常为 barcode umi 序列，reads2 为转录本序列。

单细胞分析流程

单细胞的数据处理主要包括 illumina 数据碱基识别，数据质控过滤，生成 feature-count 矩阵等过程。这些过程都可以使用 cellranger 完成。同样也可以使用第三方的 STARsolo，Alevin等软件完成。如果是 10x genomics 数据，使用 cellranger 软件更加方便。

目前一种测序是使用 cell ranger 进行前期处理，得到表达矩阵，后续分析使用增加灵活的下游处理工具。

一、数据转换 bcl2fastq

代码语言：javascript复制

#案例数据，illuminca测序原始数据
wget https://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-1.2.0.tar.gz
wget https://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-simple-1.2.0.csv
tar -zxvf cellranger-tiny-bcl-1.2.0.tar.gz

#illumina basecalling
cellranger mkfastq --id=output --run=cellranger-tiny-bcl-1.2.0 --csv=cellranger-tiny-bcl-simple-1.2.0.csv --localcores=12 --localmem=32
#查看结果
ll output/outs/fastq_path/H35KCBCXY/test_sample/
cd output/outs/fastq_path/H35KCBCXY/test_sample/
#查看index文件
zcat test_sample_S1_L001_I1_001.fastq.gz | head 
#查看read1
zcat test_sample_S1_L001_R1_001.fastq.gz | head
# 查看read2
zcat test_sample_S1_L001_R2_001.fastq.gz | head #需要的序列

二、数据质控

使用 fastqc 进行质控

代码语言：javascript复制

mkdir qc
#对read2 进行质控
ll pbmc_1k_v3_fastqs/pbmc_1k_v3*_R2_001.fastq.gz
fastqc -t 12 -f fastq -o qc pbmc_1k_v3_fastqs/pbmc_1k_v3*_R2_001.fastq.gz

10x流程中很少做质控和过滤，质控是因为现在测序质量基本比较好，过滤是因为read1和read2内容不一样，所以不像其他测序一样处理，不能按同一个标准来，低质量adapter之类的，要是过滤后可能read1和read2不匹配。cellranger对文件夹中的文件名统一，更正后clean文件名干扰流程。另cellranger中有参数选项--r1-length设置，在分析前可以帮我们切短。

除非我们自己比对，这样不做处理，要是数据过于差就做处理。

三、生成矩阵 count

这里使用 10x Genomics 官方分析软件 Cell Ranger 对原始数据进行数据质量统计，并比对参考基因组。Illumina 双末端测序结果中，使用V3试剂 read1包含 16bp barcode 和 12bp UMI，用于区分不同 RNA 分子，一个 mRNA 分子将会被一个 UMI 标记）；Read2 为 cDNA 序列片段。

Cell Ranger 调用 STAR 软件将 read2 比对到参考基因组，生成 bam 文件，然后使用 GTF 文件中的坐标位置，将比对上的 reads 分类为外显子、内含子或基因间区的 reads。这个过程与传统 bulk RNAseq 类似。

Cell Ranger接着会过滤和校正barcodes与UMIs，Cell barcodes要求与数据库已知的barcode序列完全一致，只允许有一个错配且这个错配只能出现在低质量碱基处。接着这个错误将会被校正，而其他不满足该条件的 barcodes 将会被过滤。UMI 不允许是单寡聚链、不允许含有 N、不允许含有质量值低于<10 的碱基，否则会被过滤。如果某个 UMI 与更高计数的 UMI只有一个错配且它们有相同的 barcode 和 gene id，则它会被校正成较高计数的那个 UMI。只有有效验证过 barcode 和 UMI 的 reads 才用于 UMI couting。

将每个 barcode 的每个 gene id 对应的 UMI 去重，计算 unique UMI 的数量作为该细胞该基因的表达量。最终生成 cell barcode 表达矩阵。

3.1 运行软件

count 是 cellranger 最主要也是最重要的功能：完成细胞和基因的定量，也就是产生了我们用来做各种分析的基因表达矩阵。

代码语言：javascript复制

cellranger count --id=run_count_1kpbmcs --fastqs=pbmc_1k_v3_fastqs --sample=pbmc_1k_v3 
    --transcriptome=/share/home/xiehs/15.singlecell/refdata-gex-GRCh38-2020-A/ 
    --localcores=12 --localmem=32

常见选项参数：

--id <ID>：

--transcriptome：参考序列目录

--fastqs：测序数据目录

--no-bam：不生成 bam 文件，节约存储

--nosecondary：不要次优的比对

--r1-length：截取 reads1

--r2-length：截取 reads2

--localcores：CPU 核心数

--localmem：内存大小

3.2 结果文件

最终结果在 run_count_1kpbmcs/outs/目录下，里面包含下面文件。

代码语言：javascript复制

├── analysis
│ ├── clustering
│ ├── diffexp
│ ├── pca
│ ├── tsne
│ └── umap
├── cloupe.cloupe
├── filtered_feature_bc_matrix
│ ├── barcodes.tsv.gz
│ ├── features.tsv.gz
│ └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
│ ├── barcodes.tsv.gz
│ ├── features.tsv.gz
│ └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5
└── web_summary.html

cell ranger count 分析结果文件

文件名字	描述
web_summary.html	网页简版报告以及可视化
metrics_summary.csv	CSV format 数据摘要
possorted_genome_bam.bam	比对排序后的bam文件
possorted_genome_bam.bam.bai	比对排序后的bam文件索引
filtered_gene_bc_matrices	过滤掉的 barcode 信息的结果文件，后面Seurat, Monocle 的输入文件
filtered_gene_bc_matrices_h5.h5	过滤掉的 barcode 信息 HDF5 format
raw_gene_bc_matrices	原始 barcode 信息结果文件夹
raw_gene_bc_matrices_h5.h5	原始 barcode 信息 HDF5 format
analysis	数据分析文件夹，里面包括聚类，差异分析，pca 以及 tsne
molecule_info.h5	cellranger aggr 整合多样本时用到
cloupe.cloupe	Loupe Browser 输入文件

四、结果解读

4.1 结果统计

使用浏览器打开 web_summary.html 文件，查看结果统计信息。

详细文档：

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/summary

结果统计信息

统计结果包括测序统计，细胞统计，比对统计，样品信息统计。例如图中列出本次测序捕获到 1231 个细胞，每个细胞平均 reads 数为 54104，每个细胞中检测到基因的中位值为 3235。这些信息可以反映出本次实验是否成功。

比对信息中列出与参考序列的比对情况，只有具有较高的比对率才行。

4.2 细胞计数质控（cell QC）

细胞计数质控是单细胞数据分析中非常重要的内容。因为 10xgenomics 是采用液滴型的捕获细胞方法。最终每个 GEM 中可能包含一个细胞，两个细胞，多个细胞以及无细胞。这其中只有一小部分的液滴包含珠状物和一个完整细胞。绝大部分是无细胞。但有些 RNA 会从死细胞或破损细胞中漏出来，这些游离在环境中的 RNA 可能会被空细胞捕获到，最终空细胞也会测序到 reads。

在单细胞分析中需要将这些多细胞以及空细胞都过滤掉，只对单细胞结果进行分析。那么如何判断是否为单细胞呢？

一个简单的判断是根据 reads 数据的多少，例如空细胞 reads 条数少，单细胞正好，多细胞最多。但要考虑液滴大小、扩增效率和测序环节中的波动会导致”背景”和真实细胞最终获得的文库大小变化很大，使得区分哪些文库来源于背景哪些来源于真实细胞变得复杂。

Cell Ranger 3.0 引入了一种改进的细胞计数算法，该算法能够更好地识别低 RNA 含量的细胞群体，特别是当低 RNA 含量的细胞与高 RNA 含量的细胞混合时。

该算法分为两步：

在第一步中，使用之前的 Cell Ranger 细胞计数算法识别高 RNA 含量细胞的主要模式，使用基于每个 barcode 的 UMI 总数的 cutoff 值。Cell Ranger 将期望捕获的细胞数量 N 作为输入(see —expect-cells)。然后将 barcodes 按照各自的 UMI 总数由高到低进行排序，取前 N个 UMI 数值的 99%分位数为最大估算 UMI 总数(m)，将 UMI 数目超过 m/10 的 barcodes 标记的细胞视为真实细胞。

在第二步中，选择一组具有低 UMI 计数的 barcode，这些 barcode 可能表示“空的”GEM 分区，建立 RNA 图谱背景模型。利用 Simple Good-Turing smoothing 平滑算法，对典型空GEM 集合中未观测到的基因进行非零模型估计。最后，将第一步中未作为细胞计数的barcode RNA 图谱与背景模型进行比较，其 RNA 谱与背景模型存在较大差异的 barcode 用于区分包含细胞的 barcode 和空 barcode。