Mapping hundreds of single reads to reference genome with bwa mem using shell script

Question

[t86dan]

Hi guys,

I have 100+ reads ready to be mapped to the reference genome using the bwa mem algorithm. The read files are in .fastq.gz format (ie. Blank.fastq.gz, TR8.fastq.gz, JE997.fastq.gz, IF898.fastq.gz, U1.fastq.gz, U2.fastq.gz, U3.fastq.gz,... etc.). I am very new to programing and tried running a shell script that lists all the samples and then calls another shell script as follows:

#SamplesScript

./BwaMapping Blank
./BwaMapping TR8
./BwaMapping JE997
./BwaMapping IF898
./BwaMapping U1
./BwaMapping U2
./BwaMapping U3
./BwaMapping U4

...Etc with all the remaining samples (100+ total)

#BwaMapping
p=$1;


 # input files
f1=/Path_to_read_files/${p}.fastq.gz;

REFERENCE=/Path_to_reference/reference.fa;



 # map the reads with bwa mem


bwa mem -t 2 -w 600 -R "@RG\tID:${p}\tSM:${p}\tPL:ILLUMINA" ${REFERENCE} ${f1} > ${p}_bwa_mem.sam &

So when I run ./SamplesScript it calls the other script and does the bwa mem command with each sample. The problem is that it starts to process all of the samples at once, running out of memory and obviously crashing. What do I need to change to the script so that it runs the samples in succession one after the other and not all at the same time?

Answer 1

[WouterDeCoster]

The easiest would probably be to just use a for loop in your shell, for example

for fastqfile in *.fastq.gz
do
bwa mem -t 2 -w (etc etc)
done

[t86dan]

It worked, thank you!

Answer 2

[igor]

In your script, you have & at the end of bwa mem command. That sends the process to the background. Remove that and your commands will wait for the previous ones to finish before they start.

[t86dan]

Oh I didn't know that was causing the problem! Well I will keep it in mind from now on thanks!