Automatize alignment and sorting
2
0
Entering edit mode
23 months ago
jomagrax ▴ 40

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 • 654 views
ADD COMMENT
4
Entering edit mode
23 months ago

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

ADD COMMENT
0
Entering edit mode

That was indeed the problem thank you!!

ADD REPLY
2
Entering edit mode
23 months ago
Mensur Dlakic ★ 27k

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.

ADD COMMENT

Login before adding your answer.

Traffic: 2106 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6