CellRanger的前前后后（scRNA）

2022-03-20 本文已影响0人林枫bioinfo

实操一些数据，随机选择上传至SRA数据库SRP278381中的三个bam文件进行下载。

cd ~/scRNA
mkdir raw_data_bam; cd raw_data_bam
nohup wget -c https://sra-download.ncbi.nlm.nih.gov/traces/sra74/SRZ/012492/SRR12492083/MUX7093.bam &
nohup wget -c https://sra-download.ncbi.nlm.nih.gov/traces/sra58/SRZ/012492/SRR12492084/MUX7094.bam &
nohup wget -c https://sra-download.ncbi.nlm.nih.gov/traces/sra52/SRZ/012492/SRR12492085/MUX7095.bam &

下载完毕，得到文件如下：

利用集成到cellranger软件中的二进制文件bamtofastq，完成由bam到fastq的转换。
运行nohup sh ~/script/bamtofastq.sh > /dev/null 2>&1 &，其中bamtofastq.sh中的内容如下：

#!/bin/bash
bamtofastq=~/bioinfo_soft/cellranger/cellranger-6.1.2/lib/bin/bamtofastq   # v1.3.1
mkdir raw_data_fq

for sample in $(ls raw_data_bam)
do
{
    ${bamtofastq} raw_data_bam/${sample} raw_data_fq/$(basename ${sample} .bam) > bamtofastq_${sample}.log 2>&1
} &
done

这里运行命令时，将输出重定向到/dev/null，是因为已将log文件的输出步骤写在了脚本内部，for循环中使用 {} & 实现并行（此处相对简陋，目前是三个示例样本问题不大，样本量大时，同时间的高并发，会对系统造成很大的压力和损耗，处理速度也会越来越慢，应使用更为严谨专业的snakemake或者WDL实现并发，或严格控制shell每次循环的并发数目）。
利用ps --forest查看，成功生成一个子shell，并同时运行三个样本的bam->fastq转换：

运行完成后，文件名按照10x Genomics官方要求修改前缀，最终得到文件的结构如下：

此时可进行常规的cellranger count pipeline，运行nohup sh ~/script/cellranger_count.sh > /dev/null 2>&1 &，其中cellranger_count.sh中的内容如下：

#!/bin/bash
cellranger=~/bioinfo_soft/cellranger/cellranger-6.1.2/cellranger   # v6.1.2
ref=~/bioinfo_soft/cellranger/refdata-gex-GRCh38-2020-A

for sample in $(ls raw_data_fq)
do
{
    fq_path=$(dirname raw_data_fq/${sample}/*/* | uniq)
    ${cellranger} count \
        --id=${sample} \
        --transcriptome=${ref} \
        --fastqs=${fq_path}/ \
        --sample=${sample} \
        --localcores=8 \
        --nosecondary > cellranger_count_${sample}.log 2>&1
} &
done

得到三个结果文件夹：

最重要是的其中的outs文件夹，filtered_feature_bc_matrix文件夹中的三个文件，是作为下游Seurat、Monocle等软件的标准输入：

以上是一些简单的基本流程，下面深入研究下各文件。

首先是bam->fastq，以样本MUX7094为例，输出路径下首先生成名为mux7094_MissingLibrary_1_HTL2MBBXX的文件夹，代表什么呢？
官网给出的解释是[sample_name]_[library_id]_[gem_group]_[flowcell_id]。

Special tags included by 10x Genomics pipelines are required to reconstruct the original FASTQ sequences correctly. If your BAM file lacks the appropriate headers, you will get an error message: WARNING: no @RG (read group) headers found in BAM file.

从官网的描述中可知，bam->fastq的转换，bam文件必须包含特定的注释信息，@RG是关键，提取bam文件中的@RG信息，印证了文件夹命名的含义：

而bam的主体与常规测序bam文件也有所不同，有25列，前11列作为必需字段与标准格式保持一致。后14列有些是cellranger特有的，有些是标准格式的可选字段：

文件夹下的fastq文件的命名规则为[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz，其中Read Type分为I1（Sample index read (optional)），R1（Read 1），R2（Read 2）。本实例样本的建库方式为3' v2：

查看转回的fastq文件，符合文库构成，即I1中是8bp的Sample Index，R1中是26bp的10x Barcode+UMI，R2是98*bp的Insert。

转回了可靠的fastq文件后，用常规cellranger count pipeline对数据进行表达定量，R包Seurat的Read10X()读取filtered_feature_bc_matrix文件夹中的三个文件，生成的对象为稀疏矩阵，CreateSeuratObject()创建Seurat S4对象，进而可以进行各类下游分析……

同时pipeline生成了html网页版的报告，分成Summary和Analysis两个panel（由于我在cellranger count脚本中加了--nosecondary参数，因而Analysis中没有出现降维、差异表达分析等内容，这部分内容我们还是交给Seurat来做）：

这里的一个重点是测序饱和度（Sequencing Saturation），测序饱和度是反映当前测序量与文库复杂度相关性的指标，其大小主要取决于测序深度和文库复杂度。

The fraction of reads originating from an already-observed UMI. This is a function of library complexity and sequencing depth. More specifically, this is the fraction of confidently mapped, valid cell-barcode, valid UMI reads that had a non-unique (cell-barcode, UMI, gene). The formula for calculating this metric is as follows: Sequencing Saturation = 1 - (n_deduped_reads / n_reads)，where n_deduped_reads = Number of unique (valid cell-barcode, valid UMI, gene) combinations among confidently mapped reads，n_reads = Total number of confidently mapped, valid cell-barcode, valid UMI reads.

参考：

CellRanger的前前后后（scRNA）

猜你喜欢

热点阅读