非正常数据读取——之仅有bam文件

生信技能树jimmy

发布于 2020-12-11 10:00:09

2.6K0

发布于 2020-12-11 10:00:09

文章被收录于专栏：单细胞天地

分享是一种态度

我所理解的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 datahttps://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP261877

scRNA-seq data

如图，基本了解到做了三组实验，每组两个重复。但关键作者提供的是bam文件格式，就要想办法转换为cellranger所需要的文件格式。

1、下载数据

还是使用ebi数据库下载，详见实战1的介绍或10X的单细胞转录组原始数据也可以在EBI下载。
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软件，也可见使用ebi数据库直接下载fastq测序数据的改进脚本

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-10xhttps://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

three type fastq

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。如图介绍，介绍