实操一些数据，随机选择上传至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.
参考：
- https://kb.10xgenomics.com/hc/en-us/articles/360058600992
- https://support.10xgenomics.com/docs/bamtofastq
- https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam#bam-bc-tags
- https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/fastq-input
- https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/summary
- https://kb.10xgenomics.com/hc/en-us/articles/115003646912-How-is-sequencing-saturation-calcul
网友评论