10X单细胞测序之cellranger介绍

简介

  目前10X单细胞测序算是测序行业最热门的方向之一，它可以在低测序深度的情况下一次性的获得成千上万的细胞及其每个细胞内的基因表达情况，对了解细胞异质性和新的细胞类型非常有利。官网介绍的功能和优势如下：

- 鉴定和识别新细胞类型；
- 分析并了解细胞异质性及其对生物系统的影响；
- 用单细胞RNA-seq进行细胞表型鉴定，无需预先选择靶标即可鉴定新的靶标、生物标志物以及细胞类型和状态；
- 评价同一细胞内的mRNA和细胞表面蛋白表达谱；
- 同时在数以万计的细胞中进行高通量和高分辨率的功能遗传筛选；
- 评估单个CRISPR扰动的综合基因表达表型

  既然10X单细胞优势这么大，那么了解它数据的分析过程就十分有必要。单细胞分析的内容主要包括数据拆分、细胞定量、降维聚类、差异、富集和注释。这次我们主要讨论10X GENOMICS公司为单细胞转录组量身打造的软件——cellranger。这款软件能帮助我们实现分析内容的前两部步，其中还有最重要的一步——定量。

cellranger功能介绍

  cellranger功能强大，像数据拆分cellranger mkfastq、细胞定量cellranger count、组合分析cellranger aggr、二次分析cellranger reanalyze等分析都可以完成。

数据拆分

  在测序过程中经常会出现两个文库上同一个lane，或者一个文库上不同lane的情况，对于这种情况，使用cellranger的mkfastq工具就可以实现数据的拆分。有以下两种运行方式：

cellranger mkfastq --id test \  #指定输出目录的名字
--run run_directory \  #illumina 下机bcl文件夹的路径。
--csv simple.csv

cellranger mkfastq --id test --run  run_directory --samplesheet samplesheet.csv

  该命令其实是对illumina提供的拆分数据的bcl2fastq命令的一个封装，需要样本名称，index等信息，支持两种格式，一种就是illlumina常规的samplesheet.csv文件，还有一种是10X genomics定制的一种简化版的csv格式。第一种如下所示，格式复杂：

image.png

  第二种含有三列信息，一列指定lane ID, 第二列指定样本名称，第三列指定index的名称，10X genomics的每个index代表4条具体的oligo序列。推荐使用第二种简化版的csv文件，因为cell ranger可以识别所用试剂盒版本，然后自动化的调整reads长度。

Lane,Sample,Index
1,test,SI-GA-A3

  拆分好后的目录结果如下：

├── fastq_path
│   ├── H35KCBCXY
│   │   └── test
│   │       ├── test_S1_L001_I1_001.fastq.gz #index序列
│   │       ├── test_S1_L001_R1_001.fastq.gz
│   │       └── test_S1_L001_R2_001.fastq.gz

细胞定量

  如果手里头数据已经是拆分好的fq.gz数据，就可以直接进行该部分分析。cellranger提供count工具实现测序数据中细胞和基因的定量，产生后续分析用到的基因表达矩阵，运行方式如下：

 cellranger count --id=sample345 \  #输出文件夹名
                   --transcriptome=/opt/refdata-cellranger-GRCh38-3.0.0 \  #参考基因组路径
                   --fastqs=/home/jdoe/runs/HAWT7ADXX/outs/fastq_path \  #fq.gz所在路径
                   --sample=mysample \ #样品名
                   --expect-cells=1000 #期望的细胞回收数，默认3000
                   --force-cells  强制收回的细胞数，如果Fraction Reads in Cells过低，每个样本过滤后细胞的reads数占总reads数比例过低，表明检测到的reads大部分不在细胞中，可提高该参数来增加细胞数
                   --r1-length  read1保留的长度，不能低于26bp
                   --r2-length read2保留的长度
                   --localcores cpu数
                   --localmem 内存

  输出文件夹内容：

.outs
├── analysis #数据分析文件夹
│   ├── clustering #聚类，图聚类和k-means聚类
│   ├── diffexp #差异分析
│   ├── pca #主成分分析线性降维
│   └── tsne #非线性降维信息
├── cloupe.cloupe #Loupe Cell Browser 输入文件
├── filtered_feature_bc_matrix #过滤后矩阵信息，后续降维聚类分析使用的路径
│   ├── barcodes.tsv.gz  #过滤（细胞中表达基因的数目在阈值内）后的细胞总数文件
│   ├── features.tsv.gz #所有细胞表达基因的并集
│   └── matrix.mtx.gz #坐标文件，第一列是基因行号，第二列是细胞列号，第三列是基因表达量，仅仅是列出有表达量的基因
├── filtered_feature_bc_matrix.h5  #过滤掉的barcode信息HDF5 format
├── metrics_summary.csv #CSV format数据摘要
├── molecule_info.h5 #UMI信息，aggregate的时候会用到的文件
├── raw_feature_bc_matrix #过滤前矩阵信息
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├──  possorted_genome_bam.bam #比对文件
├──  possorted_genome_bam.bam.bai  #索引文件
├── raw_feature_bc_matrix.h5 #原始barcode信息HDF5 format
├── web_summary.html #网页形式的报告以及可视化
└── *_gene_bar.csv_temp #过程文件

细胞合并

  对于这个功能，官网如此介绍：当进行涉及多个GEM Well的大型研究时，运行cellranger请分别从每个GEM Well收集fastq数据，然后使用cellranger aggr汇集结果。也就是说，需要分样进行cellranger count分析，然后再使用aggr进行合并。

cellranger aggr --id=AGG123 \ #输出文件夹名
                  --csv=AGG123_libraries.csv                   --normalize=mapped #测序文库深度校正的标准化方法，默认mapped
                  --nosecondary  是否跳过二次分析（降维、聚类和可视化）

  csv文件需要两列文件，第一列是GEM well唯一的标识ID，第二列是运行count产生的molecule_info.h5文件，格式如下：

library_id,molecule_h5
LV123,/opt/runs/LV123/outs/molecule_info.h5
LB456,/opt/runs/LB456/outs/molecule_info.h5
LP789,/opt/runs/LP789/outs/molecule_info.h5

  输出结果，目录结构和count基本一致：

Outputs:
- Aggregation metrics summary HTML:         /home/jdoe/runs/AGG123/outs/web_summary.html
- Aggregation metrics summary JSON:         /home/jdoe/runs/AGG123/outs/summary.json
- Secondary analysis output CSV:            /home/jdoe/runs/AGG123/outs/analysis
- Filtered feature-barcode matrices MEX:    /home/jdoe/runs/AGG123/outs/filtered_feature_bc_matrix
- Filtered feature-barcode matrices HDF5:   /home/jdoe/runs/AGG123/outs/filtered_feature_bc_matrix.h5
- Unfiltered feature-barcode matrices MEX:  /home/jdoe/runs/AGG123/outs/raw_feature_bc_matrix
- Unfiltered feature-barcode matrices HDF5: /home/jdoe/runs/AGG123/outs/raw_feature_bc_matrix.h5
- Unfiltered molecule-level info:           /home/jdoe/runs/AGG123/outs/raw_molecules.h5
- Barcodes of cell-containing partitions:   /home/jdoe/runs/AGG123/outs/cell_barcodes.csv
- Copy of the input aggregation CSV:        /home/jdoe/runs/AGG123/outs/aggregation.csv
- Loupe Cell Browser file:                  /home/jdoe/runs/AGG123/outs/cloupe.cloupe

二次分析

  如果第一次count分析结果不理想，如检测到的reads大部分不在细胞中，可以在二次分析中调参数重新分析，并且使用的数据不再是fq.gz数据，速度更快，使用方法：

cellranger reanalyze --id=AGG123_reanalysis \  #输出目录名
                       --matrix=AGG123/outs/filtered_feature_bc_matrix.h5 \   #count分析后的h5文件，如果--force-cells中的数目大于第一次过滤得到的细胞数，需要使用raw_feature_bc_matrix.h5文件
                       --params=AGG123_reanalysis.csv  #重分析的csv文件，包含一些参数及值，不含表头，主要是聚类、PCA分析参数。如果为空，则按照参数的默认值分析。因为10X单细胞后续的聚类差异等分析会使用其他软件（Seurat、Monocle）分析，所以该文件为空就行。
                       --force-cells  3000  #强制收回的细胞数

  输出结果如下：

Outputs:
- Secondary analysis output CSV:          /home/jdoe/runs/AGG123_reanalysis/outs/analysis_csv
- Secondary analysis web summary:         /home/jdoe/runs/AGG123_reanalysis/outs/web_summary.html
- Copy of the input parameter CSV:        /home/jdoe/runs/AGG123_reanalysis/outs/params_csv.csv
- Copy of the input aggregation CSV:      /home/jdoe/runs/AGG123_reanalysis/outs/aggregation_csv.csv
- Loupe Cell Browser file:                /home/jdoe/runs/AGG123_reanalysis/outs/cloupe.cloupe

参考文章

使用cell ranger拆分10X单细胞转录组原始数据
10X单细胞测序分析软件:Cell ranger