This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Automatize alignment and sorting

I want to automatize alignment and sorting for several fastq files.

When I run the command without automation for one sample It complies fine:

bwa mem -t 15 -R "@RG\tID:id\tSM:sm\tPL:illumina\tLB:lb" hg19 sample1-R1.fastq sample1-R2.fastq | samtools view -@ 15 -b - | samtools sort -@ 15 - > sample1-R1.bam

Whereas It doesn't happens the same when I try to automatize It.

for i in *R1.fastq ; do bwa mem -t 15 -R "@RG\tID:id\tSM:sm\tPL:illumina\tLB:lb" hg19 $i ${i/R1/R2} | samtools view -@ 15 -b - | samtools sort -@ 15 > ${i/_R1.fastq/.bam} ; done

[M::bwa_idx_load_from_disk] read 0 ALT contigs
[W::bseq_read] the 1st file has fewer sequences.
[W::bseq_read] the 1st file has fewer sequences.
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 15 -R @RG\tID:id\tSM:sm\tPL:illumina\tLB:lb ../10.reference_files/hg19 ../12.xengsort_fastq/CDX1-both.R1.fastq ../12.xengsort_fastq/CDX1-both.R2.fastq
[main] Real time: 1.549 sec; CPU: 1.550 sec

When I check the fastq files they have been edited somehow

-rw-rw-r-- 1 user user         5251  Jun 21 09:26 CDX1-both.R1.fastq
-rw-rw-r-- 1 user user 12140549  Jun 21 08:56 CDX1-both.R2.fastq

Where is the mistake?

bash aligment bwa samtools
• 1,370 views • link
written by jomagrax • jomagrax

2 answers

there seems to be a typo in your cmdline.

in the part where you write the samtools sort to a file it is written ${i/_R1.fastq/.bam} while it should be ${i/-R1.fastq/.bam}

what now happens is that bash can't do the replacement (as the 'regex' does not match) and will just write to $i (without any name change), and thus overwrites your input file

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

That was indeed the problem thank you!!

• 0 views • link
written by jomagrax • jomagrax

When in doubt, echo the commands before you actually execute them.

for i in *R1.fastq ; do echo "bwa mem -t 15 -R "@RG\tID:id\tSM:sm\tPL:illumina\tLB:lb" hg19 $i ${i/R1/R2} | samtools view -@ 15 -b - | samtools sort -@ 15 > ${i/_R1.fastq/.bam}" ; done

This will give the following output:

bwa mem -t 15 -R @RGtID:idtSM:smtPL:illuminatLB:lb hg19 CDX1-both.R1.fastq CDX1-both.R2.fastq | samtools view -@ 15 -b - | samtools sort -@ 15 > CDX1-both.R1.fastq

Obviously you are over-writing the .fastq file, and lieven.sterck already explained to you why that's the case.

• 0 views • link
written by Mensur Dlakic Mensur Dlakic

Log in to answer this question.

More posts like this