系统学习单细胞转录组测序scRNA-Seq(三)

2016/12/18评论4,078

scRNA的上游处理CellRanger，CellRanger就是基于10X genomics，从fastq生成表达矩阵的工具，然后有了表达矩阵，就可以下游利用Seurat等R包分析
CellRangert软件大小有955M，下载地址：https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest

官方教程：https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger

第一步 BCL2fastq

cellranger可以直接对illumina每个flowcell的BCLs文件(base call files)操作，使用cellranger mkfastq ，它整合了illumina平台的bcl2fastq。虽然我们一般拿到的是fastq文件，但为了从头构建流程，这一步也学一下。

测试文件下载

http://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-1.2.0.tar.gz

然后构建一个csv文件，其中包括lane、sample、index信息

这部分内容只是展示，实际分析我们自己是不会用到的

# BCL文件
wget http://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-1.2.0.tar.gz
tar zxvf cellranger-tiny-bcl-1.2.0
mv cellranger-tiny-bcl-1.2.0 tiny_bcl
# csv文件，例如
cat cellranger-tiny-bcl-simple-1.2.0.csv
# Lane,Sample,Index
# 1,test_sample,SI-GA-A3


cellranger mkfastq --id=tiny-bcl 
                     --run=tiny_bcl 
                     --csv=cellranger-tiny-bcl-simple-1.2.0.csv
# configured to use 56GB of local memory
# 对服务器内存有要求，不足的话跑不完

结果文件会存放在之前我们用id参数命名的tiny-bcl目录中

ls -l tiny-bcl/outs/fastq_path/
#H35KCBCXY
#Reports
#Stats
#Undetermined_S0_L001_I1_001.fastq.gz
#Undetermined_S0_L001_R1_001.fastq.gz
#Undetermined_S0_L001_R2_001.fastq.gz

也可以使用--qc参数进行质控，这样就会生成一个qc_summary.json文件，主要包含这些内容

第二步定量

计数分为两类：

细胞计数：利用输入数据的barcode将每个细胞的reads区分出来
UMI计数：只有比对到转录组可信度高的reads才能用于UMI计数

例如得到的fastq文件如下：

sample_S1_L001_R1_001.fastq.gz
sample_S1_L001_R2_001.fastq.gz
#一般来讲R1比R2小是正常的，因为我们要的read1只有二十几bp

一般来讲，read1 为十几或二十几bp(比如：前16bp是barcode，紧接着是12bp的UMI，这样我们需要的read1的有用部分就28bp)【这部分质控效果可能会比较差，因为我们只需要前28bp的序列，却对整条read1进行质控，后面的序列主要是一些重复的poly T】；read2 为91bp【这部分因为全是需要的转录本，因此质控效果会比较好】

比对要用到基因组和注释文件，因此先要准备这些

wget http://cf.10xgenomics.com/supp/cell-exp/refdata-cellranger-GRCh38-3.0.0.tar.gz
# 大概11G
# 其他的参考数据下载可以看：https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest

cellranger count --id=tiny-bcl 
                                 --transcriptome=refdata-cellranger-GRCh38-3.0.0 
                                 --fastqs=tiny-bcl/outs/fastq_path 
                                 --sample=sample 
                                  --expect-cells=1000
# --localcores指定多线程
# --localmem指定内存

结果会创建一个sample123 的目录，其中包括

Outputs:
- Run summary HTML:                         /opt/tiny-bcl/outs/web_summary.html
- Run summary CSV:                          /opt/tiny-bcl/outs/metrics_summary.csv
- BAM:                                      /opt/tiny-bcl/outs/possorted_genome_bam.bam
- BAM index:                                /opt/tiny-bcl/outs/possorted_genome_bam.bam.bai
- Filtered feature-barcode matrices MEX:    /opt/tiny-bcl/outs/filtered_feature_bc_matrix
- Filtered feature-barcode matrices HDF5:   /opt/sample345/outs/filtered_feature_bc_matrix.h5
- Unfiltered feature-barcode matrices MEX:  /opt/tiny-bcl/outs/raw_feature_bc_matrix
- Unfiltered feature-barcode matrices HDF5: /opt/tiny-bcl/outs/raw_feature_bc_matrix.h5
- Secondary analysis output CSV:            /opt/tiny-bcl/outs/analysis
- Per-molecule read information:            /opt/tiny-bcl/outs/molecule_info.h5
- CRISPR-specific analysis:                 null
- Loupe Cell Browser file:                  /opt/tiny-bcl/outs/cloupe.cloupe

analysis是数据分析目录，其中又包含了clustering聚类(graph-based & k-means聚类)；diffexp 差异分析；pca主成分分析(线性降维)；tsne 非线性降维
cloupe.cloupe 可以作为可视化软件Loupe Cell Browser输入文件
filtered_feature_bc_matrix 是过滤后的barcode信息，其中包含了barcodes.tsv.gz; features.tsv.gz; matrix.mtx.gz
filtered_feature_bc_matrix.h5 是过滤后的barcode信息HDF5 format
metrics_summary.csv 是csv format数据摘要
molecular_infos.h5 是UMI信息
raw_feature_bc_matrix 是原始的barcode信息
possorted_genome_bam.bam 是比对文件[另外还有配套的索引]

定量过程对内存要求较高，真实项目一般需要1-2d

【可选】多样本整合

如果有多个样本，可将它们的定量结果整合，先将每个样本信息写好

#都写在lib.csv中
library_id, molecule_h5
sample1,/sample1/sample1/outs/molecule_info.h5
sample2,/sample2/sample2/outs/molecule_info.h5

然后使用

cellranger aggr --id=sample1vssample2 
                              --csv=lib.csv 
                              --normalize=mapped

输出结果中与上面类似

[可选] 重新调整

对细胞重新分群

因为barcode是区别细胞的，因此需要提供barcodes.csv

# barcodes.csv

Barcode
AAACCTGAGCCACGCT-1
AAACCTGAGTGGCGAC-1
AAACCTGAGCGTTACG-1
...

cellranger reanalyze    --id=tiny-bcl 
                                          --matrix=filtered_feature_bc_matrix.h5 
                                          --barcodes=barcodes.csv

适当调整细胞数

计数之前调整

cellranger count  --id=tiny-bcl 
                                --transcriptome=refdata-cellranger-GRCh38-3.0.0 
                               --fastqs=tiny-bcl/outs/fastq_path 
                               --sample=sample 
                               --force-cells=5000

计数之后强制调整

cellranger reanalyze  --id=tiny-bcl 
                                        --matrix=raw_feature_bc_matrix.h5 
                                        --force-cells=5000

kmeans

先将预调整的cluster数量放在一个文件中

# para.csv
max_clusters,15

然后

cellranger reanalyze    --id=tiny-bcl 
                                          --matrix=filtered_feature_bc_matrix.h5 
                                          --params=para.csv

第三步格式转化

$ cellranger mat2csv filtered_feature_bc_matrix sample_gene_bar.csv

得到的文件就是行为基因，列为细胞的表达矩阵

发表评论