Samtools view Exec format error Help

Question

[aman-coasab]

Hi

I am new to sequence analysis and am working on wheat to identify progenitor specific SNPs. I have GBS sequence reads mapped against the wild progenitor of wheat and in the first step, I have separated mapped and unmapped reads from these bam files. Next, I have realigned the unmapped reads from progenitor genome (diploid) to the bread wheat genome (hexaploid) with BWA using following commands:

Command Used for converting .bam files to .fastq files:

$ samtools bam2fq mapped.bam > mapped.fastq

Command Used for aligning mapped reads against bread wheat genome:

$ bwa aln –t 6 –k 2 –o 0 <IndexedReferenceFile> mapped.fastq >  At3_m2.bam

where,

k — maximum difference in the seed (gaps + mismatches)

o — number of allowed gaps

The new mapping files were generated without any error. In the next step, I need to separate mapped and unmapped reads again, but when I am using the new bam files for next steps it is giving me following error:

[E::hts_hopen] Failed to open file At3_m2.bam
[E::hts_open_format] Failed to open file At3_m2.bam
samtools view: failed to open "At3_m2.bam" for reading: Exec format error

I have tried following command for separating mapped and unmapped reads:

$ samtools view -b -f 4 <file>.bam > <file>_unmapped.bam
$ samtools view -h -f 4 <file>.bam > <file>_unmapped.bam
$ samtools view -b -h -f 4 <file>.bam > <file>_unmapped.bam

Is something wrong with the set of commands I am using or something is wrong with the realigning? Could you please help me with this error. Thank you in advance!

[ATpoint]

Are the reads you have very short < 50-70bp? If not better use bwa memfor alignment and pipe it directly into a sorting command such as:

bwa mem (options...) | samtools sort -o sorted.bam

That saves you from large (uncompressed) sam files and gives you the sorted bam file right away, which you'll need for variant calling later.

[aman-coasab]

Hi ATpoint My read length is greater is around 130bp. I have not tried bwa mem with that. But I will surely try it. Thank you for your help:)

[WouterDeCoster]

I suspect bwa aln writes a sam file, and not a bam file.

[aman-coasab]

Hi WouterDeCoster! I did check with my commands and bwa aln actually created .sai file. This .sai file could be further converted to .sam file using bwa samse. Thank you for your help:)

Answer 1

[John Marshall]

“Exec format error” indicates that samtools has been unable to recognise the format of your file. Either you ran your bwa command with nohup or similar and somehow managed to get some binary characters into it, or there are some other commands you have not shown us that have replaced the file with an index file. What does htsfile At3_m2.bam say?

Instead of redirecting the output of bwa aln, use bwa's built-in facilities to specify the output filename:

bwa aln –t 6 –k 2 –o 0 -f At3_m2.sam <indexedreferencefile> mapped.fastq

This is the easy way to ensure that no spurious text from log messages gets into your output file due to the vagaries of nohup or other job scheduling utilities in use.

Edited to add: In the comments, the OP reminds us that this is the (old) bwa aln command that is being used. That command produces an intermediate .sai file that is internal to BWA, and that in turn is converted to SAM by the bwa samse or bwa sampe commands. So it is unsurprising that samtools produces “Exec format error” for this BWA-only intermediate .sai file.

So instead of using shell redirection, use the built-in output file facilities of both commands:

bwa aln –t 6 –k 2 –o 0 -f At3_m2.sai <indexedreferencefile> mapped.fastq
bwa samse -f At3_m2.sam <indexedreferencefile> At3_m2.sai mapped.fastq

(Or bwa sampe as appropriate.)