利用cellranger分析单细胞数据

生信喵

背景
       当前的单细胞测序主要采用 illumina 测序平台进行测序，一般为双末端测序，测序完成之后首先需要对 illumina 测序数据进行质控过滤，过滤条件与其他分析类似。需要注意的是，虽然单细胞测序也是双末端测序，但是 reads1 中通常为 barcode+umi 序列，reads2 为转录本序列。
      
       单细胞分析流程
       单细胞的数据处理主要包括 illumina 数据碱基识别，数据质控过滤，生成 feature-count 矩阵等过程。这些过程都可以使用 cellranger 完成。同样也可以使用第三方的 STARsolo，Alevin等软件完成。如果是 10x genomics 数据，使用 cellranger 软件更加方便。
       目前一种测序是使用 cell ranger 进行前期处理，得到表达矩阵，后续分析使用增加灵活的下游处理工具。

一、数据转换 bcl2fastq

1. #案例数据，illuminca测序原始数据
   
2. wget https://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-1.2.0.tar.gz
   
3. wget https://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-simple-1.2.0.csv
   
4. tar -zxvf cellranger-tiny-bcl-1.2.0.tar.gz
   
5. 
   
6. #illumina basecalling
   
7. cellranger mkfastq --id=output --run=cellranger-tiny-bcl-1.2.0 --csv=cellranger-tiny-bcl-simple-1.2.0.csv --localcores=12 --localmem=32
   
8. #查看结果
   
9. ll output/outs/fastq_path/H35KCBCXY/test_sample/
   
10. cd output/outs/fastq_path/H35KCBCXY/test_sample/
    
11. #查看index文件
    
12. zcat test_sample_S1_L001_I1_001.fastq.gz | head
    
13. #查看read1
    
14. zcat test_sample_S1_L001_R1_001.fastq.gz | head
    
15. # 查看read2
    
16. zcat test_sample_S1_L001_R2_001.fastq.gz | head #需要的序列

复制代码

二、数据质控
       使用 fastqc 进行质控

1. mkdir qc
   
2. #对read2 进行质控
   
3. ll pbmc_1k_v3_fastqs/pbmc_1k_v3*_R2_001.fastq.gz
   
4. 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 最主要也是最重要的功能：完成细胞和基因的定量，也就是产生了我们用来做各种分析的基因表达矩阵。

1. cellranger count --id=run_count_1kpbmcs --fastqs=pbmc_1k_v3_fastqs --sample=pbmc_1k_v3 \
   
2.     --transcriptome=/share/home/xiehs/15.singlecell/refdata-gex-GRCh38-2020-A/ \
   
3.     --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/目录下，里面包含下面文件。

1. ├── analysis
   
2. │ ├── clustering
   
3. │ ├── diffexp
   
4. │ ├── pca
   
5. │ ├── tsne
   
6. │ └── umap
   
7. ├── cloupe.cloupe
   
8. ├── filtered_feature_bc_matrix
   
9. │ ├── barcodes.tsv.gz
   
10. │ ├── features.tsv.gz
    
11. │ └── matrix.mtx.gz
    
12. ├── filtered_feature_bc_matrix.h5
    
13. ├── metrics_summary.csv
    
14. ├── molecule_info.h5
    
15. ├── possorted_genome_bam.bam
    
16. ├── possorted_genome_bam.bam.bai
    
17. ├── raw_feature_bc_matrix
    
18. │ ├── barcodes.tsv.gz
    
19. │ ├── features.tsv.gz
    
20. │ └── matrix.mtx.gz
    
21. ├── raw_feature_bc_matrix.h5
    
22. └── 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/ ... test/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。

      
       算法改进

4.3 单细胞亚群分类
      
       t-SNE 图
      
       差异表达基因
      
       饱和度评估

生信喵

五、Loupe Browser 可视化
       Loupe Browser 可视化视频教程：https://v.qq.com/x/page/g3330d5jsaq.html