This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BWA-MEM with an array of files

Hi,

I'm trying to run BWA-MEM with an array of files rather than manually inputting each single one... but it seems BWA is not getting the files correctly. Does anyone have some experience with this, please?

This is my command,

bamlist=$(find $PWD -name "${id}*.fastq.gz" | sort)
bwa mem ${reference} "${bamlist[@]}"

I don't really understand it because echoing it seems to be right,

bwa mem genome/gr37.fasta Sample_23/23_E_2_S214_L006_R1_001.fastq.gz Sample_23/23_E_2_S214_L006_R2_001.fastq.gz Sample_23/23_E_2_S214_L007_R1_001.fastq.gz Sample_23/23_E_2_S214_L007_R2_001.fastq.gz Sample_23/23_E_2_S214_L008_R1_001.fastq.gz Sample_23/23_E_2_S214_L008_R2_001.fastq.gz
alignment array bwa-mem
• 1,358 views • link
updated by LChart LChart • written by garcesj • garcesj

what is the (error) message you get, if any?

• 0 views • link
written by lieven.sterck lieven.sterck

None in particular (that's why I didn't indicate it). It just gives back the BWA-MEM help.

• 0 views • link
written by garcesj • garcesj

aha, but bwa mem only takes at best 1 pair of fastq files as input, not a whole list.

So to achieve this efficiently do what Pierre Lindenbaum suggests below.

• 0 views • link
written by lieven.sterck lieven.sterck

1 answer

solution 1: merge each R1 and R2 files into one. Beware, the order of the fastq files must be the same

find Sample_23 -type f -name "*.fastq.gz" | grep _R1_ | sort | xargs cat > Sample_23/concat.R1.fq.gz
find Sample_23 -type f -name "*.fastq.gz" | grep _R2_ | sort | xargs cat > Sample_23/concat.R2.fq.gz

and then

bwa mem genome/gr37.fasta Sample_23/concat.R1.fq.gz Sample_23/concat.R2.fq.gz

solution 2:

use a workflow manager to align each pair of fastq in parallel and later, merge the bams using samtools merge

solution 3:

don't reinvent the wheel and use an existing workflow like sarek https://github.com/nf-core/sarek

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

For solution 1, because bwa streams the data, you can save some space by using a named pipe:

mkfifo Sample_23/concat.R1.fq.gz
find Sample_23 -type f -name "*.fastq.gz" | grep _R1_ | sort | xargs cat > Sample_23/concat.R1.fq.gz &
...
• 0 views • link
written by LChart LChart

Log in to answer this question.

More posts like this