美文网首页
修改bam文件的染色体信息

修改bam文件的染色体信息

作者: jiarf | 来源:发表于2022-05-03 11:11 被阅读0次

    背景

    在运行wes的数据流程的时候,生成marked.bam文件之后的染色体是1 2 这样的,但是下载的gatk的人的基因组的fa文件的染色体是chr1 chr2这样的信息,所以进行后续的recal和bqsr的时候会报错,这是由于自己之前bwa比对的时候比对的那个基因组文件不是自己生成的,是服务器上别人生成的bwa的index,发现那个染色体是没有chr的,十所以我就那么跑了,但是现在需要给生成的bam文件添加chr信息。

    解决方式1-没有成功

    首先也是查到了几个别人能解决的方案

    samtools view test_marked.bam | awk '$3="chr"$3' | awk -F ' ' '$7=($7=="=" || $7=="*"?$7:sprintf("chr%s",$7))' | tr ' ' '\t' | samtools view -b > output.bam
    

    这里建议大家先提取一个1000行的小的bam文件,因为行数太少的话就只有头文件了,所以这里我用的是1000行,感觉运行速度还行,
    但是这个代码我报错了

    [E::sam_parse1] missing SAM header
    [W::sam_read1] Parse error at line 199
    [main_samview] truncated file.
    

    我也是没有找到解决的方法,大家有能用的话可以试试,另外就是不懂不要瞎改人家的代码,我刚开始也以为直接改第三列的染色体信息就行,但是依旧是报错的,说头文件有问题还是啥,逻辑上也不能单独改第三列的,还要有后面的那么多信息。

    解决方式2-没有成功

    samtools view -H test_marked.bam | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - test_marked.bam > output.bam
    

    报错了

    [E::cram_get_ref] Failed to populate reference for id 0
    

    这个我查到的解决方法是要添加环境

    export REF_PATH="$(pwd)/ref/cache/%2s/%2s/%s:http://www.ebi.ac.uk/ena/cram/md5/%s"
    export REF_CACHE="$(pwd)/ref/cache/%2s/%2s/%s"
    

    但是我加了也应用了,还是报同样的错误,
    后来查到可以用-T

    samtools view -T test_marked.bam | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - test_marked.bam > test.CHR.bam
    

    但是这样也会报错

    Segmentation fault (core dumped)  
    

    未解决

    解决方式3-成功

    samtools view -h test_marked.bam | sed -e '/^@SQ/s/SN\:/SN\:chr/' -e '/^[^@]/s/\t/\tchr/2'|awk -F ' ' '$7=($7=="=" || $7=="*"?$7:sprintf("chr%s",$7))' |tr " " "\t" | samtools view -h -b -@ 10 -S - > test.chr.bam
    

    这个也是把运行前后的bam文件检查了,那一列确实是添加进去了,没有什么错误。

    相关文章

      网友评论

          本文标题:修改bam文件的染色体信息

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