Make bam index for long chromosomes

軟體開發 · 發表 2018-10-21 17:28:56

摘要：最近遇到一個問題，在對一個物種用GATK call variant時，發現在用MarkDuplicates 命令對bam檔案標記重複，在運一小段時間後會卡住終止，命令如下： gatk --java-options "-Xmx20G" MarkDuplicate...

最近遇到一個問題，在對一個物種用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，比如ofollow,noindex" target="_blank">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 轉載請註明出處

Make bam index for long chromosomes

您可能也會喜歡…