0%

Make bam index for long chromosomes

最近遇到一个问题,在对一个物种用GATK call variant时,发现在用MarkDuplicates命令对bam文件标记重复,在运一小段时间后会卡住终止,命令如下:

gatk --java-options "-Xmx20G" MarkDuplicates \
    -I A1_sorted_RG.bam \
    -O A1_sorted_RG_marked.bam \
    -M A1.metrics \
    --CREATE_INDEX true \
    --VALIDATION_STRINGENCY SILENT

但换picard软件来标记重复时却正常运行(只标记重复,不建index),然后排除下其他原因后,发现是MarkDuplicates命令在对bam文件建index时报错 GATK是调用samtools来建index的,因此我用samtools来试下,提示如下错误:

[E::hts_idx_push] Region 536882174..536882320 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6
samtools index: failed to create index for "A1_sorted_RG_marked.bam": Numerical result out of range

大体意思是:某条染色体序列太长了,无法储存在bai index中,可以尝试csi index;这是因为samtools对于long chromosomes(2^29)序列尽管会产生一个bai index文件,但是是空的,因此是bai文件是无效的

所以如果下游软件支持csi index的话,则可以使用samtool index的-c参数来生成一个csi index

samtool index -c A1_sorted_RG_marked.bam

但是!GATK现在最新的版本是不支持csi index的。。。。(官方是这么说的,还是因为其是调用samtools的)

后来网上查了下是否还有其他软件能支持long chromosomes序列来建index,比如sambamba,功能不如samtools的多,但看了下都是很常用的功能,下载后就一个文件,设定执行权限后就能使用

./sambamba-0.6.8-linux-static index -p A1_sorted_RG_marked.bam

可惜还是无法生成正常的bai index

PS.后来又查了下似乎biobambam2软件可以对大于2^31长度的染色体序列建bai index,但是还没试过。。。


那么对于没有bai索引的bam文件,则可以使用bcftools软件来call variant (Samtools+bcftools Call SNP)

但samtools会在后续版本中取消mpileup生成VCF and BCF的功能;而vcftools防止samtools版本的更新对于VCF and BCF文件影响也新增了一个mpileup,专门用于后续的call variant分析;所以之前的samtools+vcftools方法现在可以更改为vcftools一个软件即可实现了,简单的命令如下:

bcftools mpileup -Ou -f genome.fa A1_sorted_RG_marked.bam |bcftools call -mv -o A1.vcf
  • -m 参数 对应multiallelic-caller算法,主要用于multiallelic and rare-variant calling
  • -v 参数 作用只输出snp variant位点

本文出自于http://www.bioinfo-scrounger.com转载请注明出处