0. 概述
Cell Ranger 是 10X genomics 官方提供的一套针对单细胞 RNA 测序输出结果进行比对、定量、聚类及基因表达分析的分析流程,它包含有与单细胞基因表达分析相关的四个pipelines,分别是:
cellranger mkfastq 流程
:其功能为将 Illumina 测序仪产生的 raw base call (BCL) 文件解析成 FASTQ 文件。
cellranger count 流程
:其功能为将 cellranger mkfastq 产生的或其他来源的 FASTQ 文件进行比对、过滤、barcode 计数以及 UMI 计数,并可以生成 feature-barcode 定量矩阵,随后确定细胞群并进行基因表达分析。
cellranger aggr 流程
:接受cellranger count的输出数据,将同一组的不同测序样本的表达矩阵整合在一起,比如tumor组原来有4个样本,PBMC组有两个样本,现在可以使用aggr生成最后的tumor和PBMC两个矩阵,并且进行标准化去掉测序深度的影响。
cellranger reanalyze 流程
:接受 cellranger count 或 cellranger aggr 的输出文件,对数据重新进行降维、聚类等后续分析。
以上四个pipeline 均将转录组常用比对软件 STAR 封装其中,可以输出带有细胞信息的 BAM、MEX、CSV、HDF5 及 HTML 等格式的结果。
参考:CellRanger官网
1. mkfastq 流程
如果是直接拿到的R1/R2的fastq文件,那么就直接上cellranger count。
illumina下机fastq文件命名规则:
但是如果是I1/R1/R2的数据,那就还要跑个cellranger mkfastq。
2. count流程
2.1命令
/software/cellranger-4.0.0/cellranger count \
--id=cellranger_HH \ #数据结果文件夹的名称
--fastqs=rawdata \ #原始文件路径
--sample=Con-1-10X5 \ # 原始下机数据文件的前缀 (fastq.gz文件名当中S1之前的字段)
--transcriptome=/reference/refdata-gex-GRCh38-2020-A #参考基因组文件目录路径
--description=H2022-xxxx_人-PBMC \ # 对样本信息进行描述
--localcores=16 \ # 将限制 cellranger 一次最多使用 16 个核
--localmem=64 \ # 将限制 cellranger 使用的最大内存量 64(GB)
注意:
> --localcores
用于设置程序的核心数,不设置的话,默认是使用系统所能够使用的最多的核心。(设置是为了避免cellranger耗尽所有资源)
> --localmem
:限制cellrange使用的内存数(设置是为了避免cellranger耗尽所有资源)
> 如果是裂核项目,需加上--include-introns
参数运行。因为添加此参数可以将 reads 映射到内含子区域,提高对含有大量前 mRNA 分子(如细胞核)样品的敏感性。
> 如果遇到结果不符合预期,可根据情况利用 --force-cells=Num
指定细胞数,通过调整细胞数以获得预期的结果。
> --except-cells
:期望得到的细胞数目,默认是3000个,一般大家都设置1000
最好不要去使用
--except-cells
和--force-cells=Num
这两个参数,没必要去人为参与细胞识别的过程。
2.2 结果文件
我们主要关注outs这个文件夹。outs文件夹中的结果文件主要如下
Output Files:
Outs
├── analysis # 数据分析结果文件
│ ├── clustering # 聚类文件夹,图聚类,k-means 聚类
│ ├── diffexp # 差异分析
│ ├── pca # pca 线性降维
│ ├── tsne # tsne 非线性降维
│ └── umap # umap 非线性降维
├── cloupe.cloupe # 用于 Loupe Browser 可视化和分析的输入文件
├── filtered_feature_bc_matrix # 过滤后的 barcode 矩阵信息,用于后续分析的必要文件│
│ ├── barcodes.tsv.gz
│ ├── features.tsv.gz
│ └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5 # 过滤后的 barcode 信息 HDF5 格式
├── metrics_summary.csv # 结果汇总 csv 文件
├── molecule_info.h5 # 实验相关的文库,GEM,barcode表达量,UMI 等信息的HDF5格式的文件,cellranger aggr aggregate 数据的时候会用到的文件
├── possorted_genome_bam.bam # 比对到参考基因组和转录本上并带有 barcode 信息的 reads 信息
├── possorted_genome_bam.bam.bai # possorted_genome_bam.bam 的索引文件
├── raw_feature_bc_matrix # 未经过滤的所有 barcode 的矩阵信息
│ ├── barcodes.tsv.gz
│ ├── features.tsv.gz
│ └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5 # 原始 barcode 信息 HDF5 格式
└── web_summary.html # 结果汇总 html 文件及可视化
3. out结果文件解读
3.1Unfiltered feature-barcode matrices
和Filtered feature-barcode matrices
是怎么来的
这两个文件中都包括了如下三个文件,Unfiltered目录下是所有的barcode信息,包含了细胞相关的barcoed和背景barcode,而Filtered目录下只包含细胞相关的barcode信息。下游分析我们需要的就是filtered_gene_bc_matrices 文件夹下的文件作为输入文件输入到seurat或者scanpy。
├── barcodes.tsv
├── genes.tsv
└── matrix.mtx
- Unfiltered feature-barcode matrices
(1)修剪Reads:针对 3’ 建库数据的基因表达比对,在比对之前会先对 reads 进行修剪。
cDNA 的全长结构中,在 3’ 和 5’ 端分别带有 poly-A 尾和TSO 序列结构(相对于比较长的 RNA 分子,一部分来自短 RNA 分子的 reads 可能仅包含 TSO 和 poly-A 序列的其中一种)。由于这种低复杂度的非模板序列的存在有可能混淆 reads 的映射,所以在比对之前一般会将 poly-A 尾和 TSO 序列分别从 reads 的 3’ 端和 5’ 端切除,这一步骤有助于提高分析的灵敏度和软件分析的效率。
(2)判断 reads 是否比对到了基因组
Cell Ranger 中封装了比对软件 STAR,根据转录本的注释文件 GTF 中的注释信息,使用 STAR 来判断reads 是比对到了外显子、内含子还是基因间区上,或者说来判断 reads 是否比对到了基因组上。
当一条 read 至少要有 50% 碱基序列与基因组上的外显子碱基互补配对,认为其比对到了外显子上;若 reads 未比对上外显子但与内含子相交,则认为其比对到了内含子上;否则为比对到了基因间区。若 reads 比对到了一个单一的外显子位点,但同时比对到了一个或多个非外显子位点,则优先认为该 read 比对到了外显子位点,MAPQ 为 255。
(3)判断 reads 是否比对到了转录本
Cell Ranger 通过检测 reads 比对上的外显子和内含子与转录本的相容性,进一步将 reads 与注释的转录本对齐。如下图所示,reads 根据它们是正义还是反义,以及它们是外显子还是内含子,或者它们的剪接模式是否与该基因相关的转录本注释兼容来分类。
上图中,绿色展示的是基因及基因中所包含的外显子,Transcript 1 和 Transcript 2 为基因经过可变剪切形成的两种转录本所包含的外显子。针对比对到正义链上的reads,如果 reads 比对到了一个外显子上或者比对到两个相邻的外显子上,则该 read 被分类为转录本 read(蓝色);如果 reads 比对到两个不相邻的外显子上,则该 read 被分类为外显子 read(浅蓝色);如果 reads 比对到内含子区域,则该 read 被分类为内含子 read(红色);紫色表示 reads 比对到反义链上。
小知识(单细胞测序和单核细胞测序的区别)
在默认情况下,只有蓝色的转录本 read 会被计入到 UMI 计数中。但在某些情况下,如在实验时输入的为细胞核时,未剪接的转录本有可能产生高水平的内含子序列,为了将这些内含子 read 计入,cellranger count 可以添加一个参数为 include-introns。当使用该参数时,任何比对到单个基因的 reads ---- 包括转录本 read(蓝色)、外显子 read(浅蓝色)和内含子 read(红色)都会计入 UMI 计数中。
此外,只有在基因组上有唯一比对位点的 reads 才被计入到UMI计数中。
(4)进行 UMI 计数
在计算 UMIs 之前,Cell Ranger 会试图矫正 UMI 序列中的测序错误。
在转录本上有唯一比对位点的 reads 根据他们的barcode、UMI 和比对到的基因被分成不同的组。如果两个组的 reads 拥有相同的 barcode 序列并比对到同一个基因上,但是 UMI 序列中有一个碱基不同,那么其中一个 UMI 有可能是因为测序中的碱基替换错误而引入的。在这种情况下,UMI 的reads 数量少的那一组会被更正为 UMI 的reads数量多的那组。
矫正可能的测序错误后进行 UMI 计数。
Cell Ranger 会再次按照 UMI(可能是修正后的)、barcode 和比对到的基因对 reads 进行分组。如果两组或者多组的 reads 拥有相同的 barcods 和 UMI 序列,但是比对到了不同的基因上,那么 reads 计数最高的那组比对到的基因会被进行 UMI 计数,其他的组则被舍弃掉。如果 reads 最高计数相同,则全部的组都被舍弃掉。
经过这两步过滤步骤后,每一个被统计到的barcode、UMI 和 基因都会被保存在未过滤的 feature-barcode 矩阵中,输出在 unfiltered feature-barcode matrix 文件夹中。
- Filtered feature-barcode matrices
Cell Ranger 进行细胞识别的算法直接影响了我们可以获得的高质量细胞的数量及数据质量。在基于 droplet 方法的单细胞技术中,通常认为含有细胞的液滴应该含有更多的RNA,因此其在 UMI 总量上应该与空细胞(背景噪音)存在明显的区分(也就是我们常说的 barcodes 排序图上的拐点)。然而实际上微滴间不同的扩增效率会导致一些较小的细胞跟空细胞在 UMI 总数上是相近的,无法仅通过 UMI 总数很好地区分空细胞和非空细胞。尤其当样本中混杂了不同大小的细胞,例如在肿瘤样本中一般混杂着体型较大的肿瘤细胞和较小的肿瘤浸润淋巴细胞,浸润淋巴细胞则较难与空细胞区分。
为解决这一难题,Cell Ranger 的算法采用了两步法分别基于 UMI 阈值识别高 RNA 含量细胞以及基于表达谱识别低 RNA 含量细胞。
(1)第一步,选取一个 UMI 总数的阈值,所有大于这个阈值的 barcodes 被识别为细胞。这一步保证了高 RNA 含量的 barcodes 被保留。
具体的算法是:将 UMI 计数从高到低进行排序,根据预期细胞数 N(软件默认N=3000),排名前 N 个细胞中的 99 分位 UMI 数值记为 m,将所有 UMI 计数大于 m/10 这一阈值的 barcodes 标记为高质量细胞。
(2)剩余未通过阈值的 barcodes 则进行第二步的筛选,根据与空细胞 RNA 表达谱是否存在显著差异来回收潜在的低 RNA 含量细胞。
此算法基于 Lun et al. 2019 年发表的算法 EmptyDrops。
a. 首先选取一组低 UMI 计数的 barcodes 作为背景集(来代表空细胞),用这部分 barcodes 的表达谱构建一个“背景模型”。
b. 接着将在第一步骤中所有未被标注为高质量细胞的 barcodes 依次和背景模型相比较,那些与背景模型存在显著差异的细胞会被回收到高质量细胞的行列。
例如:
下图是一群高 RNA 含量的 293T细胞和一群低 RNA 含量的 PBMC 细胞的混合样本数据。可以看到在被标记为高质量的部分出现了两个群体(由第一个拐点A 大致分开),在第二个拐点 B 附近的区域同时包含空细胞和高质量细胞,这部分细胞即为从第二步中回收出的细胞,图片中颜色的深浅代表了局部高质量细胞的比例。
所有被回收的高质量细胞的矩阵,会被输出到 filtered_feature_bc_matrix 文件夹中,根据矩阵信息进行下游分析。
3.2 网页报告解读
网页的结果分成了summary和analysis两部分, summary部分包含如下结果:
参考:
Cell Ranger 知多少?(上)
Cell Ranger 知多少?(中)
Cell Ranger 知多少(下)
【单细胞测序】 关于cellranger