CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览

评论8,283

总的来说，Cell Ranger主要的流程有：拆分数据 mkfastq、细胞定量 count、定量组合 aggr、调参reanalyze，还有一些小工具比如mkref、mkgtf、upload、sitecheck、mat2csv、vdj、mkvdjref、testrun

首先是mkfastq 拆分数据

虽然这里用不到(因为我们下载的就是fastq数据)，但是为了流程的完整还是要学习一下

目的：将每个flowcell 的Illumina sequencer's base call files (BCLs)转为fastq文件

特色： 它借鉴了Illumina出品的bcl2fastq，另外增加了：

  • 将10X 样本index名称与四种寡核苷酸对应起来，比如A1孔是样本SI-GA-A1，然后对应的寡核苷酸是GGTTTACT, CTAAACGG, TCGGCGTC, and AACCGTAA ，那么程序就会去index文件中将存在这四种寡核苷酸的fastq组合到A1这个样本
  • 提供质控结果，包括barcode 质量、总体测序质量如Q30、R1和R2的Q30碱基占比、测序reads数等
  • 可以使用10X简化版的样本信息表

它的示意流程：

CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览-图片1

两种使用方式：

  1. # 第一种
  2. $ cellranger mkfastq --id=bcl 
  3.                      --run=/path/to/bcl 
  4.                      --samplesheet=samplesheet-1.2.0.csv
  5. # 第二种
  6. $ cellranger mkfastq --id=bcl 
  7.                      --run=/path/to/bcl 
  8.                      --csv=simple-1.2.0.csv
  9. # 其中id指定输出目录的名称，run指的是下机的原始BCL文件目录
  10. # 重要的就是测序lane、样本名称、index等信息

samplesheet.csv文件就是illumina常规使用的，类似下面这种。它除了需要指定各种ID、name之外，还要根据不同的试剂盒版本调整[Reads]长度

V2试剂盒R1序列长度为26bp(包括16bp的barcode+10bp的UMI)，R2为98bp；
V3试剂盒R1序列长度为28bp(包括16bp的barcode+12bp的UMI)，R2为91bp

CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览-图片2

还有一种10X定制的简单化的csv文件，例如：

  1. Lane,Sample,Index
  2. 1,test_sample,SI-GA-A3
  3. # 其中第一列指定lane ID，第二列是样本名称，第三列是index名称

使用简化版的这个文件，可以识别使用的试剂盒版本，然后自行调整reads的长度信息

最后的结果就是三个文件：I1序列文件以及两个测序文件R1、R2

目录结构如下：

  1. - tiny-bcl/outs/fastq_path/bcl/
  2.   - Sample1
  3.     - Sample1_S1_L001_I1_001.fastq.gz
  4.     - Sample1_S1_L001_R1_001.fastq.gz
  5.     - Sample1_S1_L001_R2_001.fastq.gz

自己分析的数据也要改成这种结构存放，方便后续分析

小Tip--指定fastq文件位置

后续分析需要指定fastq位置，但是这些fastq文件可以由cellranger mkfastq得到，也可以利用s Illumina's bcl2fastq 、公共数据、10X的bamtofastq ,每种情况可能得到的fastq存放位置是不同的，那么如何根据不同情况进行指定呢？

第一种情况：

利用mkfastq或者bcl2fastq生成的文件，大概长这样

CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览-图片3
  1. # 会有这几种选择方式[注意几种参数的设置]
  2. # 1.所有mkfastq生成的样本
  3. --fastqs=MKFASTQ_ID/outs/fastq_path
  4. # 2. 多个flowcell生成的所有样本
  5. --fastqs=MKFASTQ_ID/outs/fastq_path1,MKFASTQ_ID/outs/fastq_path2
  6. # 3.所有bcl2fastq 生成的样本
  7. --fastqs=/PATH/TO/bcl2fastq_output
  8. # 4. 所有lanes上的test_sample1样本
  9. --fastqs=MKFASTQ_ID/outs/fastq_path 
  10. --sample=test_sample1
  11. # 5. lane1上的test_sample1样本
  12. --fastqs=MKFASTQ_ID/outs/fastq_path 
  13. --sample=test_sample1 
  14. --lanes=1
  15. # 6. 将test_sample1和test_sample2各自进行操作
  16. --fastqs=MKFASTQ_ID/outs/fastq_path 
  17. --sample=test_sample1,test_sample2

其实从上面的各种设置也能看出来，一开始的样本命名规则是非常重要的

第二种情况：

也是利用mkfastq或者bcl2fastq生成的文件，但是同一个样本的数据放在不同的目录

CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览-图片4
  1. # 1. 将所有SI-GA-A1样本的reads组合
  2. --fastqs=MKFASTQ_ID/outs/fastq_path 
  3. --sample=SI-GA-A1_1,SI-GA-A1_2,SI-GA-A1_3,SI-GA-A1_4
  4. # 2. 只处理SI-GA-A1样本的第一个index样本
  5. --fastqs=MKFASTQ_ID/outs/fastq_path 
  6. --sample=SI-GA-A1_1
第三种情况：

也是利用mkfastq或者bcl2fastq生成的文件，但和Reports、Stats在同一个目录

CellRanger走起(四) Cell Ranger流程概览CellRanger走起(四) Cell Ranger流程概览
  1. # 1. mkfastq得到的所有样本
  2. --fastqs=MKFASTQ_ID/outs/fastq_path
  3. # 2. bcl2fastq得到的所有样本
  4. --fastqs=/PATH/TO/bcl2fastq_output
  5. # 3. test_sample样本的所有lanes
  6. --fastqs=MKFASTQ_ID/outs/fastq_path 
  7. --sample=test_sample
  8. # 4. test_sample样本的lane1
  9. --fastqs=MKFASTQ_ID/outs/fastq_path 
  10. --sample=test_sample 
  11. --lanes=1
第四种情况：

使用 mkfastq or bcl2fastq 得到的fastq文件和Report、Stats不在同一个目录，但命名方式与之前一样，这个目录中只能看到fastq文件

CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览-图片5
  1. # 1.处理所有样本
  2. --fastqs=/PATH/TO/PROJECT_FOLDER
  3. # 2. 处理Mysample样本的所有lanes的数据
  4. --fastqs=/PATH/TO/PROJECT_FOLDER 
  5. --sample=MySample
  6. # 3. 只处理Mysample样本的lane1数据
  7. --fastqs=/PATH/TO/PROJECT_FOLDER 
  8. --sample=MySample 
  9. --lanes=1
第五种情况：

fastq命名方式变了，类似于这样：

CellRanger走起(四) Cell Ranger流程概览CellRanger走起(四) Cell Ranger流程概览

它一般是从demux流程中拆分出来的数据，但是目前被mkfastq取代，没有好的方法，需要知道样本相关的index或者oligos

  1. # 1.所有样本
  2. --fastqs=/PATH/TO/PROJECT_FOLDER
  3. # 2.所有SI-GA-A1样本
  4. --fastqs=/PATH/TO/PROJECT_FOLDER 
  5. --indices=SI-GA-A1
  6. # 3.所有SI-GA-A1样本的lane1数据
  7. --fastqs=/PATH/TO/PROJECT_FOLDER 
  8. --indices=SI-GA-A1 
  9. --lanes=1
  10. # 4.利用oligo
  11. --fastqs=/PATH/TO/PROJECT_FOLDER 
  12. --indices=AACCGTAA,CTAAACGG,GGTTTACT,TCGGCGTC
第六种情况：

数据命名与上面完全不同，因此需要自己重命名，方式就是

  1. # 这个在单细胞实战(二)中介绍过
  2. [Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz
  3. # 其中Read Type
  4. # I1: Sample index read (optional)
  5. # R1: Read 1
  6. # R2: Read 2
CellRanger单细胞转录组分析教程(四) Cell Ranger流程概览-图片6

分析时就可以直接调用了

  1. # 1.所有样本
  2. --fastqs=/PATH/TO/PROJECT_FOLDER
  3. # 2. 某个样本的所有lanes数据
  4. --fastqs=/PATH/TO/PROJECT_FOLDER 
  5. --sample=SAMPLENAME
  6. # 3. 某个样本的某个lane
  7. --sample=SAMPLENAME 
  8. --fastqs=/PATH/TO/PROJECT_FOLDER 
  9. --lanes=1

然后是count 细胞定量

这个过程是最重要的，它完成细胞与基因的定量，它将比对、质控、定量都包装了起来，内部流程很多，但使用很简单

先学会使用

每个版本要求的参数是不同的，尤其是V2与V3版本存在较大差异，这里先对V2进行了解

基本上自己需要输入的参数是：

  1. # 这是示例，不是真实数据 #
  2. cellranger count --id=sample345 
  3.                    --transcriptome=/opt/refdata-cellranger-GRCh38-1.2.0 
  4.                    --fastqs=/home/scRNA/runs/HAWT7ADXX/outs/fastq_path 
  5.                    --sample=mysample 
  6.                    --expect-cells=1000 
  7.                    --nosecondary
  8. # id指定输出文件存放目录名
  9. # transcriptome指定与CellRanger兼容的参考基因组
  10. # fastqs指定mkfastq或者自定义的测序文件
  11. # sample要和fastq文件的前缀中的sample保持一致，作为软件识别的标志
  12. # expect-cells指定复现的细胞数量，这个要和实验设计结合起来
  13. # nosecondary 只获得表达矩阵，不进行后续的降维、聚类和可视化分析(因为后期会自行用R包去做)
它的输出文件有很多
  1. Outputs:
  2. - Run summary HTML:                      /opt/sample345/outs/web_summary.html
  3. - Run summary CSV:                       /opt/sample345/outs/metrics_summary.csv
  4. - BAM:                                   /opt/sample345/outs/possorted_genome_bam.bam
  5. - BAM index:                             /opt/sample345/outs/possorted_genome_bam.bam.bai
  6. - Filtered gene-barcode matrices MEX:    /opt/sample345/outs/filtered_gene_bc_matrices
  7. - Filtered gene-barcode matrices HDF5:   /opt/sample345/outs/filtered_gene_bc_matrices_h5.h5
  8. - Unfiltered gene-barcode matrices MEX:  /opt/sample345/outs/raw_gene_bc_matrices
  9. - Unfiltered gene-barcode matrices HDF5: /opt/sample345/outs/raw_gene_bc_matrices_h5.h5
  10. - Secondary analysis output CSV:         /opt/sample345/outs/analysis
  11. - Per-molecule read information:         /opt/sample345/outs/molecule_info.h5
  12. - Loupe Cell Browser file:               /opt/sample345/outs/cloupe.cloupe
  13.  
  14. Pipestance completed successfully!

从上到下依次来看：

  • web_summary.html：官方说明 summary HTML file
  • metrics_summary.csv：CSV格式数据摘要
  • possorted_genome_bam.bam：比对文件
  • possorted_genome_bam.bam.bai：索引文件
  • filtered_gene_bc_matrices：是重要的一个目录，下面又包含了 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz，是下游Seurat、Scater、Monocle等分析的输入文件
  • filtered_feature_bc_matrix.h5：过滤掉的barcode信息HDF5 format
  • raw_feature_bc_matrix：原始barcode信息
  • raw_feature_bc_matrix.h5：原始barcode信息HDF5 format
  • analysis：数据分析目录，下面又包含聚类clustering（有graph-based & k-means）、差异分析diffexp、主成分线性降维分析pca、非线性降维tsne
  • molecule_info.h5：下面进行aggregate使用的文件
  • cloupe.cloupe：官方可视化工具Loupe Cell Browser 输入文件

一些内置软件和算法

基因组比对—是否在外显子？

利用了 STAR比对工具，这款比对工具比对速度快，灵敏度高，是ENCODE、GATK推荐使用的工具，允许基因的可变剪切。比对完之后，利用GTF文件将reads溯源回外显子区、内含子区、基因间区：如果一条read的50%以上与外显子有交集，那么就认为它在外显区；如果不在外显子区，与内含子有交集，那么就认为它在内含子区；与外显子、内含子都没有交集，那么就认为在基因间区

MAPQ 辅助判断—在外显子的正确率有多少？

如果reads比对到了一个外显子区，同时也比对到了1个或多个的非外显子区，更相信它在外显子区，然后看MAPQ值，值越大越可信，如果MAPQ的值为255的话，那么就可以非常确定它比对到了外显子区

MAPQ即mapping quality，告诉我们这个read比对到参考基因组上某个位置的可信度，它的公式是：-10logP(error)，如果这个值大于30就认为比对发生错误的概率是千分之一

转录组比对—是否特异比对？

如果上面得到的外显子区域reads同时比对上有注释转录本上的外显子，并且在同一条链上，那么认为这个reads也比对到了转录组；如果只比对到单个基因的注释信息，那么认为它是特异比对到转录组的(uniquely /confidently mapped )，这样的reads才会拿来做接下来的UMI 计数

重点和难点在于自主构建参考信息

Cell Ranger为比对和定量提供了参考基因组及注释 pre-built human (hg19, GRCh38), mouse (mm10), and ercc92 reference packages

但是很多时候，我们需要根据自己的需要，自定义一套参考信息，但需要注意以下问题：

  • 参考序列只能有很少的 overlapping gene annotations，因为reads比对到多个基因会导致流程检测的分子数更少(它只要uniquely mapped的结果)
  • FASTA与GTF比对和STAR兼容，GTF文件的第三列（feature type）必须有exon，过滤后的GTF只包含有注释的基因类型
首先利用mkgtf过滤GTF文件

先从 ENSEMBL或UCSC上下载，然后使用mkgtf

  1. cellranger mkgtf input.gtf output.gtf --attribute=key:allowable_value
  2. # 其中键值对可以指定多个，比如
  3. $ cellranger mkgtf Homo_sapiens.GRCh38.ensembl.gtf Homo_sapiens.GRCh38.ensembl.filtered.gtf 
  4.                    --attribute=gene_biotype:protein_coding 
  5.                    --attribute=gene_biotype:lincRNA 
  6.                    --attribute=gene_biotype:antisense 
  7.                    --attribute=gene_biotype:IG_LV_gene 
  8.                    --attribute=gene_biotype:IG_V_gene 
  9.                    --attribute=gene_biotype:IG_V_pseudogene 
  10.                    --attribute=gene_biotype:IG_D_gene 
  11.                    --attribute=gene_biotype:IG_J_gene 
  12.                    --attribute=gene_biotype:IG_J_pseudogene 
  13.                    --attribute=gene_biotype:IG_C_gene 
  14.                    --attribute=gene_biotype:IG_C_pseudogene 
  15.                    --attribute=gene_biotype:TR_V_gene 
  16.                    --attribute=gene_biotype:TR_V_pseudogene 
  17.                    --attribute=gene_biotype:TR_D_gene 
  18.                    --attribute=gene_biotype:TR_J_gene 
  19.                    --attribute=gene_biotype:TR_J_pseudogene 
  20.                    --attribute=gene_biotype:TR_C_gene
  21. # 这样得到的Homo_sapiens.GRCh38.ensembl.filtered.gtf结果中就不包含gene_biotype:pseudogene这部分
然后利用mkref构建参考索引
  1. # 基本使用(单个物种)
  2. cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf
  3. # 可以使用--nthreads使用多线程加速
  4. # 得到的输出结果(保存在--genome名称的目录中)
  5. ls hg19
  6. fasta/  genes/  pickle/  reference.json  star/
  7.  
  8. # 如果对于多个物种组合(本文的数据其实就应该这样组合起来)
  9. cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf 
  10.                    --genome=mm10 --fasta=mm10.fa --genes=mm10-filtered-ensembl.gtf
  11. # 得到的结果
  12.  ls hg19_and_mm10
  13.  fasta/  genes/  pickle/  reference.json  star/
如果要增加基因信息

参考链接：https://kb.10xgenomics.com/hc/en-us/articles/115003327112-How-can-we-add-genes-to-a-reference-package-for-Cell-Ranger-

第一步，在fasta/genome.fa的FASTA基础上增加序列信息；

第二步，在genes/genes.gtf的GTF基础上增加注释信息，注意格式

  1. # 每一行有9列tab分隔信息
  2. # 第一列：Chromosome 指定基因组上染色体或contig位置
  3. # 第二列：Source 这个用处不大
  4. # 第三列：Feature CellRanger软件只取exon的部分
  5. # 第四列：Start 起始位点(1-based)
  6. # 第五列：End 终止位点(1-based)
  7. # 第六列：Score 这个用处不大，建议用"."表示
  8. # 第七列：Strand feature信息在基因组的+或-链
  9. # 第八列：Frame 用处不大，建议“.”
  10. # 第九列：分号分隔的键值对，重点是transcript_id 和gene_id。gene_name可选
  11. 例如：
  12. mylocus    annotation    exon    100    200    .    +    .    gene_id "mygene"; transcript_id "mygene";

第三步，使用cellranger mkref运行更新一下

P.S. 最后得到的参考信息(包括参考基因组、注释信息)文件结构如下：
  1. # 这是官网下载的hg38数据
  2. refdata-cellranger-GRCh38-1.2.0
  3. |-- [ 222]  README.BEFORE.MODIFYING
  4. |-- [4.0K]  fasta
  5. |   `-- [2.9G]  genome.fa
  6. |-- [4.0K]  genes
  7. |   `-- [1.3G]  genes.gtf
  8. |-- [4.0K]  pickle
  9. |   `-- [ 58M]  genes.pickle
  10. |-- [ 424]  reference.json
  11. |-- [4.0K]  star
  12. |   |-- [3.0G]  Genome
  13. |   |-- [8.0G]  SA
  14. |   |-- [1.5G]  SAindex
  15. |   |-- [1.2K]  chrLength.txt
  16. |   |-- [1.9K]  chrName.txt
  17. |   |-- [3.0K]  chrNameLength.txt
  18. |   |-- [2.1K]  chrStart.txt
  19. |   |-- [ 37M]  exonGeTrInfo.tab
  20. |   |-- [ 15M]  exonInfo.tab
  21. |   |-- [526K]  geneInfo.tab
  22. |   |-- [ 909]  genomeParameters.txt
  23. |   |-- [9.1M]  sjdbInfo.txt
  24. |   |-- [7.1M]  sjdbList.fromGTF.out.tab
  25. |   |-- [7.1M]  sjdbList.out.tab
  26. |   `-- [9.4M]  transcriptInfo.tab
  27. `-- [   6]  version
  28.  
  29. 4 directories, 21 files

多个文库的整合 aggr

当处理多个生物学样本或者一个样本存在多个重复/文库时，最好的操作就是先分别对每个文库进行单独的count定量，然后将定量结果利用aggr组合起来

第一步 得到count结果

例如现在分别进行3个定量流程

  1. $ cellranger count --id=LV123 ...
  2. ... wait for pipeline to finish ...
  3. $ cellranger count --id=LB456 ...
  4. ... wait for pipeline to finish ...
  5. $ cellranger count --id=LP789 ...
  6. ... wait for pipeline to finish ...
第二步 构建Aggregation CSV

就像这样：

  1. # AGG123_libraries.csv
  2. library_id,molecule_h5
  3. LV123,/opt/runs/LV123/outs/molecule_info.h5
  4. LB456,/opt/runs/LB456/outs/molecule_info.h5
  5. LP789,/opt/runs/LP789/outs/molecule_info.h5
  6. # 其中
  7. # molecule_h5：文件molecule_info.h5 file的路径
第三步 运行aggr
  1. cellranger aggr --id=AGG123 
  2.                   --csv=AGG123_libraries.csv 
  3.                   --normalize=mapped
  4. # 结果输出到AGG123这个目录中

发表评论

匿名网友

拖动滑块以完成验证
加载失败