CellRanger初探

作者: 新欣enjoy | 来源:发表于2020-02-06 18:49 被阅读0次

    公共数据库中的SRA 单细胞转录组数据究竟包含了哪些数据?CellRanger怎样利用10x平台下机数据进行下游一系列分析?这篇文章简单记录CellRanger 包括的主要分析步骤,纯理论。

    SRA 原始数据转fastq

    公共数据库的SRA 数据需要借助fastq-dump 转为fastq文件,然后进行质控、CellRanger定量等操作。相较于普通转录组数据,原始SRA数据会得到3个fastq文件,分别是Library 的Index(8bp)文件,包括长度为26bp 的Barcode(16bp)和UMI(10bp)的Read1文件,和测序reads文件。

    conda install -c bioconda sra-tools ## 安装软件
    wkd=/home/project/single-cell/MCC
    cd $wkd/raw/P2586-4
    cat SRR_Acc_List-2586-4.txt |while read i
    do
    time fastq-dump --gzip --split-files -A $i ${i}.sra && echo "** ${i}.sra to fastq done **"
    done
    ### 单细胞数据参数为 --split-files 而不是 --split-3
    

    i7 sample index (library barcode)
    是加到Illumina测序接头上的,保证多个测序文库可以在同一个flow-cell上或者同一个lane上进行混合测序(multiplexed)。它的作用就是在CellRanger的mkfastq 功能中体现出来的,它自动识别样本index名称(例如:SA-GA-A1),将具有相同4种oligo的fq文件组合在一起,表示同一个样本。它保证了一个测序lane上可以容纳多个样本

    Barcode
    是10X特有的,用来区分GEMs,也就是对细胞做了一个标记。一般在拆分混样测序数据(demultiplexing)这个过程后进行操作,当然这也很符合原文的操作。

    UMI
    UMI就是Unique Molecular Identifier,由4-10个随机核苷酸组成,在mRNA反转录后,进入到文库中,每一个mRNA随机连上一个UMI,根据PCR结果可以计数不同的UMI,最终统计mRNA的数量。它的主要作用是,处理PCR 扩增偏差,因为起始文库很小时需要的PCR扩增次数就越多,因为越容易引入扩增误差。

    fastq文件更名
    为什么要更名?CellRanger 定量过程输入文件指定命名格式。
    如何更名?下图格式:

    # 比如,将原来的SRR7692286_1.fastq.gz改成SRR7692286_S1_L001_I1_001.fastq.gz
    # 依次类推,将原来_2的改成R1,将_3改成R2
    cat  SRR_Acc_List-9245-3.txt | while read i ;do (mv ${i}_1*.gz ${i}_S1_L001_I1_001.fastq.gz;mv ${i}_2*.gz ${i}_S1_L001_R1_001.fastq.gz;mv ${i}_3*.gz ${i}_S1_L001_R2_001.fastq.gz);done
    

    fastq 文件质控
    主要查看Q30比例。

    # 以P2586-4为例
    mkdir -p $wkd/qc
    cd $wkd/qc
    find $wkd/raw/P2586-4 -name '*R1*.gz'>P2586-4-id-1.txt
    find $wkd/raw/P2586-4 -name '*R2*.gz'>P2586-4-id-2.txt
    cat P2586-4-id-1.txt P2586-4-id-2.txt >P2586-4-id-all.txt
    
    cat P2586-4-id-all.txt| xargs fastqc -t 20 -o ./
    

    CellRanger 流程

    它主要包括四个主要基因表达分析流程:

    • cellranger mkfastq : 它借鉴了Illumina的bcl2fastq ,可以将一个或多个lane中的混样测序样本按照index标签生成样本对应的fastq文件。即对下机数据base calling files转为fastq文件。
    • cellranger count :利用mkfastq生成的fq文件,进行比对(基于STAR)、过滤、UMI计数。利用细胞的barcode生成gene-barcode矩阵,然后进行样本分群、基因表达分析
    • cellranger aggr :接受cellranger count的输出数据,将同一组的不同测序样本的表达矩阵整合在一起,比如tumor组原来有4个样本,PBMC组有两个样本,现在可以使用aggr生成最后的tumor和PBMC两个矩阵,并且进行标准化去掉测序深度的影响
    • cellranger reanalyze :接受cellranger countcellranger aggr生成的gene-barcode矩阵,使用不同的参数进行降维、聚类。属于定制化分析。

    它的结果主要是包含有细胞信息的BAM, MEX, CSV, HDF5 and HTML文件

    CellRanger 软件安装及测试实战操作参考:CellRanger走起(四) CellRanger流程概览。包括安装详细指南及其他一些小技巧。

    拆分数据 mkfastq、细胞定量 count、定量组合 aggr

    mkfastq 拆分
    目的:将每个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 mkfastq --id=bcl \
                         --run=/path/to/bcl \
                         --samplesheet=samplesheet-1.2.0.csv
    # 第二种
    $ cellranger mkfastq --id=bcl \
                         --run=/path/to/bcl \
                         --csv=simple-1.2.0.csv
    # 其中id指定输出目录的名称,run指的是下机的原始BCL文件目录
    # 重要的就是测序lane、样本名称、index等信息
    

    参考文章CellRanger走起(四) CellRanger流程概览 中解释了几种不同方式输出fastq文件的不同存放目录,可根据实际分析自行查找。

    count 定量
    这个过程是最重要的,它完成细胞与基因的定量,它将比对、质控、定量都包装了起来,内部流程很多,但使用很简单。每个版本要求的参数是不同的,尤其是V2与V3版本存在较大差异,这里先对V2进行了解。

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

    count 定量的重点和难点在于参考基因组索引的构建,可详细参考CellRanger走起(四) CellRanger流程概览,文章中还介绍了怎样自行构建参考基因组索引,及一比对过程的重要细节信息。

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

    1. 得到count结果
    2. 构建Aggregation CSV
    3. 运行aggr
    ### step1
    $ cellranger count --id=LV123 ...
    ... wait for pipeline to finish ...
    $ cellranger count --id=LB456 ...
    ... wait for pipeline to finish ...
    $ cellranger count --id=LP789 ...
    ... wait for pipeline to finish ...
    
    ## step2
    # AGG123_libraries.csv
    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
    # 其中
    # molecule_h5:文件molecule_info.h5 file的路径 
    
    ### step3
    cellranger aggr --id=AGG123 \
                      --csv=AGG123_libraries.csv \
                      --normalize=mapped
    # 结果输出到AGG123这个目录中
    

    cellranger count的结果

    count结果一般放在out目录下,主要有summary和analysis两大类

    Summary
    是一个HTML格式文件,包括实验捕获的细胞数目、检测到的基因数目、测序数据的产出与质量统计、参考基因组的比对情况。

    V2 试剂count summary 结果

    几个指标可以关注一下:

    • 左上部分中,包括了reads数、barcodes数、UMI、index、Q30等统计值

    • 左下是reads比对的比例,包括基因间区、外显子、内含子,如果比对率太低(一般认为外显子的比对率要在60%以上)

    • 右上图是利用barcodes上的UMI标签分布来估计细胞数,绿/蓝色表示细胞,灰色表示背景,其中Y轴表示每个barcode对应的UMI数量,X轴是一定数量的UMI序列所对应barcode的数量,比如上图中有1000个barcodes含有10k个UMI,细胞过滤就是通过这个图来展示的。

    • 首先明确,barcode用来区分细胞,UMI用来区分转录本;其次,barcodes数量时要大于细胞数量的(以保证每个细胞都会有barcode来进行区分)如果图线陡降说明

    另参考文章CellRanger走起(二) CellRanger out 结果 中还记录了一个手动构建物种GTF文件的案例尝试,及相应的检验方法,可以根据需要自行学习。

    总之,无论是从公共数据库下载单细胞SRA数据,还是来自10x 平台的原始下机数据,或者先经过fastq-dump转为fastq文件,然后质控,然后CellRanger count 定量,或者首先自行拆分数据mkfastq,然后构建索引进行细胞定量,及可选的定量组合等。CellRanger的总体流程见本文记录,但是其中的各个细节仍需要自行实践总结。

    参考1:CellRanger走起(二) 使用前注意事项
    参考2:CellRanger走起(三) 使用初探
    参考3:CellRanger走起(四) CellRanger流程概览
    参考4:CellRanger走起(二) CellRanger out 结果

    相关文章

      网友评论

        本文标题:CellRanger初探

        本文链接:https://www.haomeiwen.com/subject/sceuxhtx.html