14

用broad出品的软件来处理bam文件几次遇到文件头错误

报错如下:ERROR MESSAGE: SAM/BAM file input.marked.bam is malformed: SAM file doesn't have any read groups defined in the header.  The GATK no longer supports SAM files without read groups !

有些人遇到的是bam的染色体顺序不一样,还有可能是染色体的名字不一样,比如>1和>chr1的区别,虽然很傻,但是遇到这样问题的还不少!
还有一些人是遇到基因组没有dict文件,也是用picard处理一下就好。

大部分人是在GATK遇到的,我是在RNA-SeQC遇到的,不过原理都是一样的。
都是因为做alignment的时候并未添加头信息,比如:
bwa samse ref.fa my.sai my.fastq > my.sam
samtools view -bS my.sam > my.bam
samtools sort my.bam my_sorted
java -jar ReordereSam.jar I=/path/my_sorted.bam O=/path/my_reordered.bam R=/path/ref.fa
通过这个代码可以得到排序好的bam,但是接下来用GATK就会报错
java -jar GenomeAnalysisTK.jar -T DepthOfCoverage -R /paht/ref.fa -I /path/aln_reordered.bam
就是因为没有头信息,group相关信息,解决方法有两种:
bwa samse -r @RG\tID:IDa\tSM:SM\tPL:Illumina ref.fa my.sai my.fastq > my.sam
java -jar AddOrReplaceReadGroups I=my.bam O=myGr.bam LB=whatever PL=illumina PU=whatever SM=whatever
一种是比对的时候就加入头信息,这个需要比对工具的支持。
第二种是用picard工具来修改bam,推荐用这个!虽然我其实并不懂这些头文件信息是干嘛的, 但是broad开发的软件就是需要!希望将来去读PHD能系统性的学习一些基础知识!

 

06

GATK使用注意事项

GATK这个软件在做snp-calling的时候使用率非常高,因为之前一直是简单粗略的看看snp情况而已,所以没有具体研究它。

这些天做一些外显子项目以找snp为重点,所以想了想还是用起它,报错非常多,调试了好久才成功。

所以记录一些注意事项!

GATK软件本身是受版权保护的,所以需要申请才能下载使用,大家自己去broad institute申请即可。

下载软件就可以直接使用,java软件不需要安装,但是需要你的机器上面有java,当然软件只是个开始,重点是你还得下载很多配套数据,https://software.broadinstitute.org/gatk/download/bundle(ps:这个链接可能会失效,下面的文件,请自己谷歌找到地址哈。),而且这个时候要明确你的参考基因组版本了!!! b36/b37/hg18/hg19/hg38,记住b37和hg19并不是完全一样的,有些微区别哦!!!
Continue reading