SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bwa sampe error Mali Salmon Bioinformatics 14 10-27-2014 11:25 AM
bwa sampe hanging krobison Bioinformatics 6 02-13-2013 12:57 PM
too small output .sai file created by bwa aln? ..leads to bwa sampe hanging? Stina Bioinformatics 12 12-05-2012 07:32 AM
bwa sampe -c option cband Bioinformatics 1 08-07-2012 03:50 AM
BWA sampe multiple hits question sahar Bioinformatics 0 08-03-2010 01:20 PM

Reply
 
Thread Tools
Old 02-14-2016, 09:30 AM   #1
juanc
Junior Member
 
Location: ca

Join Date: Feb 2016
Posts: 7
Default bwa sampe on multiple samples at once

I'm looking to loop through many samples using the bwa sampe command.

Here is what the command looks like for one sample:

bwa sampe ref.fasta read1.sai read2.sai read1.fastq.gz read2.fastq.gz | samtools view -bS - > sample.bam

I have read1.sai, read2.sai, and ref.fasta in the same directory, but not fastq files. I have no problem using a for loop on singleton files. I've tried different things with paired-end reads on the above command but nothing seems to work (I'm not a seasoned programmer). Could someone possibly provide an example for the above command? Apologies if this question is obvious or answered elsewhere.
juanc is offline   Reply With Quote
Old 02-14-2016, 02:21 PM   #2
piet
Member
 
Location: planet earth

Join Date: Aug 2014
Posts: 21
Default

You are using short file names in your command. This will only work, if the files are located in the current working directory. If your FASTQ files reside in another directory, you must pass their full path (either relative to the working directory or absolute).
piet is offline   Reply With Quote
Old 02-15-2016, 08:52 AM   #3
juanc
Junior Member
 
Location: ca

Join Date: Feb 2016
Posts: 7
Default

Thanks for the suggestion. Would it be possible to provide a quick example?
juanc is offline   Reply With Quote
Old 02-15-2016, 09:14 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,495
Default

Are you asking for an example for @piet's comment or for the actual loop?

For @piet's comment:

Replace /provide_path_to/ part with real or relative paths from your computer.

Code:
$ bwa sampe ref.fasta read1.sai read2.sai /provide_path_to/read1.fastq.gz /provide_path_to/read2.fastq.gz | samtools view -bS - > sample.bam
GenoMax is offline   Reply With Quote
Old 02-15-2016, 10:25 AM   #5
juanc
Junior Member
 
Location: ca

Join Date: Feb 2016
Posts: 7
Default

Sorry, I meant for the actual loop.
juanc is offline   Reply With Quote
Old 02-15-2016, 10:54 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,495
Default

See this: https://www.biostars.org/p/98222/

You would need to change the actual command and parsing criteria based on your file names. If you post example file names we can help with the commands.

@piet has another example for grabbing the basenames from files which you can then use in your command: https://www.biostars.org/p/176887/#176889

Last edited by GenoMax; 02-15-2016 at 11:09 AM.
GenoMax is offline   Reply With Quote
Old 02-15-2016, 12:40 PM   #7
juanc
Junior Member
 
Location: ca

Join Date: Feb 2016
Posts: 7
Default

1. fastq naming: taxon86-READ1.fastq.gz + taxon86-READ2.fastq.gz and taxon24-READ1.fastq.gz + taxon24-READ2.fastq.gz, etc.

2. sai file naming: taxon86-read1.sai + taxon86-read2.sai and taxon24-read1.sai + taxon24-read2.sai, etc.

3. sai files are in the current working directory where bwa sampe command is executed

4. fastq files are in a different directory: /Users/juanc/projectscx/reads/taxon*/clean/

(* is used to denote multiple directories under the "reads" directory corresponding to the different sample names)

Last edited by juanc; 02-15-2016 at 05:52 PM.
juanc is offline   Reply With Quote
Old 02-15-2016, 01:31 PM   #8
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,495
Default

I am going to assume that you are in the directory where you have the sai files. Since I don't have your files use the following at your own risk :-)

NOTE: How many of these sample directories are there (I am using only one of these directories I have called it Sample_ONE).

Step 1: Grab a listing of the Read1 files for one sample

Code:
$ ls -1 /Users/juanc/projectscx/reads/sample_ONE/clean/*READ1*.fastq.gz > sample_ONE_read1_list
Step 2: See if the following works

Code:
for i in $(cat sample_ONE_read1_list | sed s'/\-READ1.fastq.gz//'); 
do bwa sampe ref.fasta $i\-READ1.fastq.sai $i\-READ2.fastq.sai /Users/juanc/projectscx/reads/sample_ONE/clean/$i\-READ1.fastq.gz /Users/juanc/projectscx/reads/sample_ONE/clean/$i\-READ2.fastq.gz | samtools view -bS - >  $i\.bam; done
GenoMax is offline   Reply With Quote
Old 02-15-2016, 02:21 PM   #9
juanc
Junior Member
 
Location: ca

Join Date: Feb 2016
Posts: 7
Default

I updated my previous post.

Last edited by juanc; 02-15-2016 at 04:04 PM.
juanc is offline   Reply With Quote
Old 02-16-2016, 04:01 AM   #10
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,495
Default

Disclaimers in the last post still apply

Step 1: Grab a listing of the Read1 files

Code:
$ ls -1r /Users/juanc/projectscx/reads/taxon*/clean/*READ1*.fastq.gz > taxon_read1_list
Step 2: Verify that the taxon_read1_list looks good

Code:
$ more taxon_read1_list
Step 3: I think the following should work
Code:
for i in $(cat taxon_read1_list | sed s'/\-READ1.fastq.gz//'); 
do bwa sampe ref.fasta $i\-READ1.sai $i\-READ2.sai /Users/juanc/projectscx/reads/$i/clean/$i\-READ1.fastq.gz /Users/juanc/projectscx/reads/$i/clean/$i\-READ2.fastq.gz | samtools view -bS - >  $i\.bam; done
GenoMax is offline   Reply With Quote
Old 02-16-2016, 09:29 AM   #11
juanc
Junior Member
 
Location: ca

Join Date: Feb 2016
Posts: 7
Default

So I executed the commands and the output report was telling me that it was looking for the .sai files in the "clean" directories. I decided to do steps 1-4 (below) in those directories and move over output files to a different directory for downstream analysis (not shown). Steps 3-4 were generously provided by @GenoMax.

Step 1: Generate alignments for READ1
Code:
for i in /Users/juanc/projectscx/reads/taxon*/clean/*READ1.fastq.gz;
do bwa aln /Users/juanc/projectscx/mapping/ref.fasta $i > ${i%READ1.fastq.gz}read1.sai; done
Step 2: Generate alignments for READ2
Code:
for i in /Users/juanc/projectscx/reads/taxon*/clean/*READ2.fastq.gz;
do bwa aln /Users/juanc/projectscx/mapping/ref.fasta $i > ${i%READ2.fastq.gz}read2.sai; done
Step 3: Create a list of READ1.fastq files
Code:
ls -1r /Users/juanc/projectscx/reads/taxon*/clean/*READ1*.fastq.gz > taxon_read1_list
Step 4: Join read pairs and convert sam to bam
Code:
for i in $(cat taxon_read1_list | sed s'/\-READ1.fastq.gz//');
do bwa sampe /Users/juanc/projectscx/mapping/ref.fasta $i\-READ1.sai $i\-READ2.sai $i\-READ1.fastq.gz $i\-READ2.fastq.gz | samtools view -bS - >  $i\.bam; done

Last edited by juanc; 02-16-2016 at 09:37 AM.
juanc is offline   Reply With Quote
Reply

Tags
bwa, samtools

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:36 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO