Science相关 杂

2021-05-18 SAMFormatException: E

2021-05-18  本文已影响0人  AsuraPrince

问题:SAMFormatException: Error parsing SAM header. @RG line missing SM tag. 

Line:@RG    ID:S143; File /public4/zhangxt/KWL/07Combam/S143.bam; Line number 165

RG line 有问题,缺少SM tag.

计划用addhead策略解决(2021-05-07没有bam文件没有RG);但是报相同错误!!!所以行不通。

网上查询到两种策略:

1.  addreplacerg (https://www.jianshu.com/p/d9d8471624cd);实操发现,他是改变每一行最后的RG标识,生成的bam文件,再次运行报同样错误!!

2.  reheader (https://www.biostars.org/p/249376/  ;samtools view -H your.bam | sed 's,^@RG.*,@RG\tID:None\tSM:None\tLB:None\tPL:Illumina,g' | samtools reheader - your.bam > your.new.bam );可以改变Line:@RG;成功解决问题!!!

2021-05-19

新发现一个问题,就是合并好的文件虽然通过reheader修改@RG行规避了上面的问题,但是具体进行picard MarkDuplicates进行去重复时,一般几个bam文件合并的总bam文件会同时保留原来的RG标识,所以导致与@RG会不一致(因为reheader对每行进行修改,只会修改@RG行;而addreplacerg恰恰相反,智慧修改每行);所以此时需要借助addreplacerg修改每行的RG标识,统一成与@RG一样!!!

#####导致最初始问题的根本原因,是在hisat2比对时,通过--rg-id S99参数添加了RG:ID,却没有定义SM等标识。

有两种解决方案:

1.hisat2比对时,完全不用管RG的问题,不添加任何RG标识,这样多个bam文件合并成总bam后,picard MarkDuplicates实际不会报错,但是gatk HaplotypeCaller.sh时会因为RG问题报错,可以通过addhead策略(2021-05-07没有bam文件没有RG)解决,一步解决问题!!!PS:感觉addhead放在picard MarkDuplicates是更规范的流程,也会解决问题。

2.hisat2比对时,通过hisat2 --rg一次性添加所有后面可能用到的标识(示例:hisat2 --rg ID:C3MF6ACXX.L001 --rg SM:501 --rg PL:ILLUMINA --rg LB:lib-501 --rg PU:C3MF6ACXX.1.NoIndex -p 10 -q -x ./ref/GRCh38.fa -1 501N-1_1.fq.gz -2 501N-1_2.fq.gz -S 501N_hst2.sam);这样最初问题就不会报错了,但是每行RG不统一问题应该还是应该存在,所以samtools merge是不是有参数在合并过程中统一成想要RG

正如所料,是有的!samtools merge -c (Combine @RG headers with colliding IDs [alter IDs to be distinct]),这里会把要合并的bam文件中,名字是相同的合成一个,如Test1.bam(RG为S1)、Test2.bam(RG为S1)、Test3.bam(RG为S1),则合并好的Test1-c.bam只会有S1,但是如果为Test1.bam(RG为S1)、Test2.bam(RG为S1)、Test3.bam(RG为S2),则Test1-c.bam会有两种标识S1和S2,因为它只能合并相同的。

而 -h (Copy the header in FILE to <out.bam> [in1.bam])为的就是在合并的bam文件中区分来源,在合并好的文件中,原来是啥现在还是啥;但是某些情况下为了区分做改变!!!即如Test1.bam(RG为S1)、Test2.bam(RG为S1)、Test3.bam(RG为S1),则合并好的Test1-c.bam只会有S1,S1xxx,S1xxxx,在不指定时,默认用的是和h参数!!! ;

而-r(Attach RG tag (inferred from file names)),更他妈简单粗暴,根本不管每个文件内RG是啥,你文件名叫啥,合并好的文件中我叫啥,如Test1.bam(RG为S1)、Test2.bam(RG为S1)、Test3.bam(RG为S1),则合并好的Test1-c.bam的RG分别是Test1、Test2、Test3),

还是必须还得通过addreplacerg在统一一下合并好的每行RG???(samtools addreplacerg -R S7 -m overwrite_all -o S9-hafter.bam S9-h.bam)

目前重新比对是不可能了,只能通过addreplacerg曲线救国了,麻烦一点,就麻烦一点吧。生信的坑为啥一个都没放过我,挨个踩了个遍!!!

上一篇 下一篇

猜你喜欢

热点阅读