最近遇到一个问题，在对一个物种用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转载请注明出处