Question: Need Coding To Run Bwa Mem In Batch Mode
1
gravatar for HG
5.8 years ago by
HG1.1k
Germany
HG1.1k wrote:

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.

perl awk • 3.8k views
ADD COMMENTlink modified 5.8 years ago • written 5.8 years ago by HG1.1k
4
gravatar for Devon Ryan
5.8 years ago by
Devon Ryan92k
Freiburg, Germany
Devon Ryan92k wrote:
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.

ADD COMMENTlink written 5.8 years ago by Devon Ryan92k

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

ADD REPLYlink written 5.8 years ago by HG1.1k

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.

ADD REPLYlink modified 5.8 years ago • written 5.8 years ago by Devon Ryan92k

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
ADD REPLYlink modified 5.8 years ago • written 5.8 years ago by HG1.1k

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.

ADD REPLYlink written 5.8 years ago by Devon Ryan92k

Thank you so much ....

ADD REPLYlink modified 5.8 years ago • written 5.8 years ago by HG1.1k
2
gravatar for Pierre Lindenbaum
5.8 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum124k wrote:

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
ADD COMMENTlink written 5.8 years ago by Pierre Lindenbaum124k

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

ADD REPLYlink written 5.8 years ago by HG1.1k

Yes, you would need to do that.

ADD REPLYlink written 5.8 years ago by Devon Ryan92k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1148 users visited in the last hour