我所理解的cellranger软件理想原始输入数据就是SRA格式,然后利用sra-tools分为read、barcode+UMI、index三个fastq.gz文件。最后直接利用cellranger即可。但总会发现文献提供的数据格式并非如此,就要花费一些心思了。
- 文献:https://www.nature.com/articles/s43018-020-00139-8
- 单细胞测序数据:https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP261877
scRNA-seq data
https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP261877
scRNA-seq data
如图,基本了解到做了三组实验,每组两个重复。但关键作者提供的是bam文件格式,就要想办法转换为cellranger所需要的文件格式。
1、下载数据
- 还是使用ebi数据库下载,详见实战1的介绍。
- https://www.ebi.ac.uk/ena/browser/view/PRJNA632854
cat > bam.txt
fasp.sra.ebi.ac.uk:/vol1/run/SRR117/SRR11798249/ICBtreated_Brca2null_rep1.bam
fasp.sra.ebi.ac.uk:/vol1/run/SRR117/SRR11798250/ICBtreated_Brca1null_rep2.bam
fasp.sra.ebi.ac.uk:/vol1/run/SRR117/SRR11798251/ICBtreated_Brca1null_rep1.bam
fasp.sra.ebi.ac.uk:/vol1/run/SRR117/SRR11798257/ICBtreated_Parental_rep2.bam
fasp.sra.ebi.ac.uk:/vol1/run/SRR117/SRR11798258/ICBtreated_Parental_rep1.bam
fasp.sra.ebi.ac.uk:/vol1/run/SRR117/SRR11798259/ICBtreated_Brca2null_rep2.bam
conda activate download
cat fq.txt |while read id
do ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
conda deactivate
rawdata
conda环境配置可见实战1,就是使用下ascp软件
2、bam转fastq
2.1 special bam of cellranger
-
这里的bam序列比对结果文件应该是作者使用cellranger后的产生的比对结果。
cellranger - 经过一番查询后,知道了cellranger产生的bam文件里是带有barcode与UMI的,储存在tag标签里。
https://bioinformatics.stackexchange.com/questions/7096/bam-to-gene-expression-matrix-umi-counts-per-gene-per-cell-10x
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam
samtools view ICBtreated_Parental_rep1.bam | less -SN
samtools view ICBtreated_Parental_rep1.bam | head -3 | tr "\t" "\n" | cat -n
bam tag
相关标签具体解释见结尾,这里知道 CB或者CR代表barcode、UB或者UR代表UMI即可。
2.2 cellranger bamtofastq
- 知道上述知识点后,我一开始想手动提取下bam文件里的barcode与UMI序列,也查了很多linux字符处理方法。但这样之后还要自己组装fastq,并且使header一致。总之很费事。
- 后来知道 cellranger也有
bamtofastq
功能,就试试看和其它软件的转换有什么不同。 - 如下图是一个bam文件的处理结果。其很强大的功能:自动根据bam文件,生成配套的三种文件。而且还修改为规范命名!
I1
,R1
,R2
的含义见实战1
ICBtreated_Brca1null_rep1.bam
I1
代表的index序列,一般是用于区分混合样品的,二代测序通配。这里为空,表示bam文件里没有。而且的确也不需要,就提供一个空序列文件即可。
- 批量处理
cat > bamtofastq.sh
bin=/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger
cat name.list |while read id
do
$bin bamtofastq $id ./fastq/${id}
done
bash bamtofastq.sh
3、cellranger count
- 经过上一步的处理,接下来就比较简单了。
find /home/shensuo/test/fastq/ | grep bam/out | grep -v XX/bam > fq.txt
bin=/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger
db=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/
target=bamtofastq
cat fq.txt |while read id
do
echo $bin count --id=${id:0-39:14} \
--localcores=4 \
--transcriptome=$db \
--fastqs=$id \
--sample=$target \
--expect-cells=3000 \
--nosecondary
done
cat > cellranger.sh
/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger count --id=ICBtreated_Brca1null_rep1 --localcores=10 --transcriptome=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/ --fastqs=/home/shensuo/test/fastq/ICBtreated_Brca1null_rep1.bam/output_0_1_HG3HHDRXX --sample=bamtofastq --expect-cells=13009 --nosecondary
/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger count --id=ICBtreated_Brca1null_rep2 --localcores=10 --transcriptome=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/ --fastqs=/home/shensuo/test/fastq/ICBtreated_Brca1null_rep2.bam/output_0_1_HG3HHDRXX --sample=bamtofastq --expect-cells=21789 --nosecondary
/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger count --id=ICBtreated_Brca2null_rep1 --localcores=10 --transcriptome=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/ --fastqs=/home/shensuo/test/fastq/ICBtreated_Brca2null_rep1.bam/output_0_1_HG3HHDRXX --sample=bamtofastq --expect-cells=18684 --nosecondary
/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger count --id=ICBtreated_Brca2null_rep2 --localcores=10 --transcriptome=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/ --fastqs=/home/shensuo/test/fastq/ICBtreated_Brca2null_rep2.bam/output_0_1_HG3HHDRXX --sample=bamtofastq --expect-cells=16731 --nosecondary
/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger count --id=ICBtreated_Parental_rep1 --localcores=10 --transcriptome=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/ --fastqs=/home/shensuo/test/fastq/ICBtreated_Parental_rep1.bam/output_0_1_HG3HHDRXX --sample=bamtofastq --expect-cells=13292 --nosecondary
/home/shensuo/biosoft/cellranger/cellranger-4.0.0/bin/cellranger count --id=ICBtreated_Parental_rep2 --localcores=10 --transcriptome=/home/shensuo/biosoft/cellranger/test/refdata-gex-mm10-2020-A/ --fastqs=/home/shensuo/test/fastq/ICBtreated_Parental_rep2.bam/output_0_1_HG3HHDRXX --sample=bamtofastq --expect-cells=20718 --nosecondary
nohup bash /home/shensuo/test/cr_out/fq.txt &
#需要使用全路径
耐心等待结果即可。估计一夜肯定是需要的。
此外其中有几个注意点,具体如下
- 一开始想只用while语句运行,发现还是不能尽善尽美。就echo下,再根据实际情况修改,保存为cellranger.sh
- 关于
--localcores=
设置,可根据实际情况。我是设置10,之后也单独尝试了32,24,20。结果还是24还比较适合目前得到服务器环境的最大承受,如果想尽快跑完的话。 - 关于
--expect-cells
设置,主要参考原文献的测序结果细胞数。
paper result
4、导出最终结果
-
如前所述,cellranger count结果很多,主要是需要其中如下图的三个文件(每个样品)
three type result
cat name.list | while read id
do
mkdir ./out/$id
cp $id/outs/filtered_feature_bc_matrix/* ./out/$id
done
接下来就可以,将数据导入到R中,用seurat等包进行下游分析了。
附:cellranger indexed bam
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam
如官网介绍,Chromium cellular and molecular barcode information for each read is stored as TAG fields in this bam(produced by cellranger)
cellular and molecular barcode分别对应我们之前说的barcode与UMI序列.前者用来区分不同GEMs,也就是对细胞做了一个标记;后者用于表示基因文库大小,即每种mRNA一个特定的UMI。
如图介绍,介绍
-
CB
、CR
、CY
表示barcode,一般是16个碱基; -
UB
、UR
、UY
表示UMI,一般是10个碱基。 -
R
一般代表原始测序数据,Y
代表质量分数,而B
代表校正后的R
,可能对应碱基质量分数太低等因素。一般来说R
与B
都是相同的。
网友评论