美文网首页Single-cell单细胞转录组单细胞测序
cellranger实战2:非正常数据之仅有bam文件

cellranger实战2:非正常数据之仅有bam文件

作者: 小贝学生信 | 来源:发表于2020-11-22 20:28 被阅读0次

    我所理解的cellranger软件理想原始输入数据就是SRA格式,然后利用sra-tools分为read、barcode+UMI、index三个fastq.gz文件。最后直接利用cellranger即可。但总会发现文献提供的数据格式并非如此,就要花费一些心思了。

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

    1、下载数据

    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

    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

    image.png
    如官网介绍,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。
    如图介绍,介绍
    • CBCRCY表示barcode,一般是16个碱基;
    • UBURUY表示UMI,一般是10个碱基。
    • R一般代表原始测序数据,Y代表质量分数,而B代表校正后的R,可能对应碱基质量分数太低等因素。一般来说RB都是相同的。

    相关文章

      网友评论

        本文标题:cellranger实战2:非正常数据之仅有bam文件

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