BWA step fails on large genomes, because samtools index requires the "-c" (CSI index) flag

[lucventurini]

Observed while trialing the pipeline together with @hesterjvs on H. vulgare, with its genome over 4GBs.
The fix would be to add a "-c" flag to the samtools index command. However, this will result in the BAM file having a .csi, rather than .bai, index.
The analysis was done using the latest dev branch, with the Docker profile (revision: 435812b [dev]).

Key part of the log output:

ERROR ~ Error executing process > 'bwa (Sample_5_S5_L002_R1_001)'

Caused by:
  Process `bwa (Sample_5_S5_L002_R1_001)` terminated with an error exit status (1)

Command executed:

  bwa aln -t 8 bwa_index/*.fasta Sample_5_S5_L002_R1_001.combined.prefixed.fq.gz -n 0.04 -l 32 -k 2 -f "Sample_5_S5_L002_R1_001.combined.prefixed.fq.sai"
  bwa samse -r "@RG\tID:ILLUMINA-Sample_5_S5_L002_R1_001\tSM:Sample_5_S5_L002_R1_001\tPL:illumina" bwa_index/*.fasta "Sample_5_S5_L002_R1_001.combined.prefixed.fq".sai Sample_5_S5_L002_R1_001.combined.prefixed.fq.gz | samtools sort -@ 8 -O bam - > "Sample_5_S5_L002_R1_001".sorted.bam
  samtools index "Sample_5_S5_L002_R1_001".sorted.bam

Command exit status:
  1

Command output:
  (empty)

Command error:
  [bwa_aln] 17bp reads: max_diff = 2
  [bwa_aln] 38bp reads: max_diff = 3
  [bwa_aln] 64bp reads: max_diff = 4
  [bwa_aln] 93bp reads: max_diff = 5
  [bwa_aln] 124bp reads: max_diff = 6
  [bwa_aln] 157bp reads: max_diff = 7
  [bwa_aln] 190bp reads: max_diff = 8
  [bwa_aln] 225bp reads: max_diff = 9
  [bwa_aln_core] calculate SA coordinate... 37.53 sec
  [bwa_aln_core] write to the disk... 0.00 sec
  [bwa_aln_core] 7289 sequences have been processed.
  [main] Version: 0.7.17-r1188
  [main] CMD: bwa aln -t 8 -n 0.04 -l 32 -k 2 -f Sample_5_S5_L002_R1_001.combined.prefixed.fq.sai bwa_index/Hordeum_vulgare.IBSC_v2.dna.toplevel.fasta Sample_5_S5_L002_R1_001.combined.prefixed.fq.gz
  [main] Real time: 124.051 sec; CPU: 41.598 sec
  [bwa_aln_core] convert to sequence coordinate... 4.41 sec
  [bwa_aln_core] refine gapped alignments... 0.91 sec
  [bwa_aln_core] print alignments... 0.01 sec
  [bwa_aln_core] 7289 sequences have been processed.
  [main] Version: 0.7.17-r1188
  [main] CMD: bwa samse -r @RG\tID:ILLUMINA-Sample_5_S5_L002_R1_001\tSM:Sample_5_S5_L002_R1_001\tPL:illumina bwa_index/Hordeum_vulgare.IBSC_v2.dna.toplevel.fasta Sample_5_S5_L002_R1_001.combined.prefixed.fq.sai Sample_5_S5_L002_R1_001.combined.prefixed.fq.gz
  [main] Real time: 18.397 sec; CPU: 5.824 sec
  [bam_sort_core] merging from 0 files and 8 in-memory blocks...
  [E::hts_idx_push] Region 543480053..543480129 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 "Sample_5_S5_L002_R1_001.sorted.bam": Numerical result out of range

[apeltzer]

Ok you're really hitting borders that I never saw before here 👍

Thanks for that, I'll see what we can do about it - should be an easy fix. Do you know by any chance whether I can use a CSI index for all kinds of references in general? That would make the indexing easier...

[apeltzer]

I fixed this in #147 , waiting for @jfy133 to review.

If you want, @lucventurini, I can add you in here as a contributor too - happy to get these kinds of bug reports and fix them to make the overall pipeline better 👍

[jfy133]

Sorry to make it more complicated, but it might be safer for the moment to make this an option to turn on.

Java-based downstream tools like IGV remains shaky with CSI support (due to htsjdk): samtools/htsjdk#1047

Currently, the majority of the aDNA field don't work on really large genomes (although it's awesome people are starting to), so as large genome support is still hit-and-miss I personally consider it safer to keep CSI generation as an optional but not default.

[apeltzer]

Hm, could also make it dependent on the input files size? e.g. if the genome size is > 4GB, I just run with the -c option, in the other case stay with the default? Would that help?

[jfy133]

Yes, that would be sufficient.

[apeltzer]

Did the changes in the updated PR - this automatically sets itself to using CSI indexing whenever a reference FastA file is larger than 3.5GB

[jfy133]

@lucventurini as we don't have a travis test for the .csi case (and I assume we can't set this one up because file sizes are too big, correct @apeltzer?), could you manually test that commit?

[apeltzer]

Would be possible to test this using this here:

nextflow run apeltzer/eager -r fix-samtools-idx [your_other_options]

[apeltzer]

Yes, @jfy133 we are in the process of getting bigger testcases set up, but for now we're bound to our smaller testcases ...

[lucventurini]

Hi all,
thanks for the amazing bug fixing turnaround! @apeltzer , unfortunately in the past years I have been working with genomes far larger than barley ... it's a whole brave world over here, if you ever want to dip your toes :-)

I will be testing the pipeline immediately and report the result.