Need Coding To Run Bwa Mem In Batch Mode

Question

[HG]

Hi every one i was trying to mapping multiple fastq files with multiple references in a single batch mode. for this i did like this

# make an index of our all resequences file
for file in ./*.fasta
    do
    echo $file 
    bwa index $file
    done

#Now take two fastq file and start bwa mem for mapping

for read in./*.fastq
do
echo $read
bwa mem $file $name $name > $aln.sam
done

The problem here for input two fastq file one after another. for example my .fastq file name like this

ResetH16_S2_L001_R1_001.fastq

ResetH16_S2_L001_R2_001.fastq

Only difference in R1/R2 can anyone help me for this coding how can i take input two file one after another according to name

Thank you advance.

Answer 1

[Devon Ryan]

for R1 in ./*R1*.fastq
do
    echo $R1
    R2=`echo $R1 | sed 's/_R1_/_R2_/'`
    bwa mem $file $R1 $R2 > $aln.sam
done

This is dependent on you already defining $aln and $file, the latter of which you appear to not be doing in the way you likely want in your current post (i.e., you're indexing multiple files but only aligning to the last one). You probably want to nest the loops.

[HG]

yes you are right it is not working as a like .Please let me know the solution. its not generating the .sam file

[Devon Ryan]

The general idea is something like:

for file in ./*.fasta
    do
    echo $file 
    bwa index $file
    for R1 in ./*R1*.fastq
    do
        echo $R1
        R2=`echo $R1 | sed 's/_R1_/_R2_/'`
        bname=`echo $R1 | sed 's/_R1_.\+//'`
        bwa mem $file $R1 $R2 > $bname.$file.sam
    done
done

Or something like that, since I can't really test things perfectly. The s/_R1_.\+// thing, in case you're not familiar with regular expressions, means "search for a portion that starts with '_R1_' and replace everything from there to the end of the input with nothing (i.e., delete it)".

Oh, I should also mention that you likely don't need the repeated "./" stuff.

[HG]

Thanks for help but still some error could you please help me out :

./Ecoli_jj1886.fasta
[bwa_index] Pack FASTA... 0.03 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 1.51 seconds elapse.
[bwa_index] Update BWT... 0.02 sec
[bwa_index] Pack forward-only FASTA... 0.03 sec
[bwa_index] Construct SA from BWT and Occ... 0.44 sec
[main] Version: 0.7.6a-r433
[main] CMD: bwa index ./Ecoli_jj1886.fasta
[main] Real time: 2.686 sec; CPU: 2.040 sec
./ResetH16_S2_L001_R1_001.fastq
map2.sh: 10: map2.sh: cannot create ./ResetH16_S2_L001../Ecoli_jj1886.fasta.sam: Directory nonexistent

[Devon Ryan]

Yeah, the unneeded "./" are mucking things up. Just:

for file in *.fasta
do
    echo $file 
    bwa index $file
    for R1 in *R1*.fastq
    do
        echo $R1
        R2=`echo $R1 | sed 's/_R1_/_R2_/'`
        bname=`echo $R1 | sed 's/_R1_.\+//'`
        bwa mem $file $R1 $R2 > $bname.$file.sam
    done
done

That's a bit more likely to work.

[HG]

Thank you so much ....

Answer 2

[Pierre Lindenbaum]

try to only read the _R1_ files

ls *.fastq | grep _R1_ | while read F
do
R=`echo $F | sed 's/_R1_/_R2_/'`
bwa mem $file ${F} ${R} > $aln.sam
done

[HG]

Do i need to define $aln earlier because its not generation the final aln.sam file

[Devon Ryan]

Yes, you would need to do that.