for loop over bwa aln

Question

[phylofun]

[h.mon]

I can't spot the issue. What happens? Are there error messages?

Unless you reads are shorter than 70bp, you should be using bwa mem instead of bwa aln.

[phylofun]

There is an error message when I run bwa sampe, but I think the issue begins with bwa aln. No error message occurs when I run bwa aln; however, the .sai output file goes to the current working directory instead of following the path in READ1_list, and the sample name is not written for the .sai file (e.g. READ1_listread1.sai instead of sample_name-read1.sai)

[WouterDeCoster]

Put an echo in your loop to see how it looks like:

 for i in $(cat READ1_list | sed s'/\-READ1.fastq.gz//' | sort) 
do
  echo "bwa sampe reference/reference.unaligned.fasta $i\-READ1.sai $i\-READ2.sai $i\-READ1.fastq.gz $i\-READ2.fastq.gz \
  | samtools view -bS - > $i\-pe.bam"
done

Please note that you could also do the following and pipe directly to samtools sort, without samtools view:

bwa sample ref.fa R1 R2 | samtools sort -o alignment.bam

[JC]

I see two problems in your code:

1. The second alignment outputs a read1.sai file, it should be read2.sai
2. The final BAM file requires files named as READ but your alignments are read (file names are generally case-sensitive)

Also it is hard to find the errors without any detail, are you running this on Linux/Mac or Windows?, what are the names of the input files?, what is inside the READ1_list and READ2_list files?

[5heikki]

for i in READ1_list

That shouldn't even work. You should instead e.g.

for i in $(cat READ1_list); do ...

Just try e.g.

for i in READ1_list; do echo "$i"; done

vs.

for i in $(cat READ1_list); do echo "$i"; done

Only the latter works in bash (as intended)..

But anyway, I wouldn't even use a for loop for this but instead just define a function and use it with find and GNU parallel. Also, why are you aligning the reads separately? Surely it would make more sense to align read pairs properly in the same runs?

Answer 1

[Joe]