SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Convert Sam file with incorrect fields to Fastq using HTSeq AnneBiton Bioinformatics 1 06-06-2013 07:03 AM
Split Large FASTQ file in small FASTQ files with user defined number of reads Windows deepbiomed Bioinformatics 3 04-04-2013 07:14 AM
How do you pull a bunch of reads by ID? dacotahm Bioinformatics 1 02-13-2013 07:58 AM
Count unique reads in a FASTQ file id0 Bioinformatics 2 10-24-2012 01:01 AM
BWASW more reads in the output SAM file than in the input file nanto Bioinformatics 2 09-18-2012 12:41 AM

Reply
 
Thread Tools
Old 06-28-2013, 03:17 PM   #1
pbm13
Junior Member
 
Location: NH

Join Date: Jan 2013
Posts: 5
Default Use SAM file to pull reads from FASTQ

Hi Folks,

I have a SAM file with unpaired reads (originally from a FASTQ) and I would like to use it to pull the read and its pair from the FASTQ file - does anyone know if there is a script out there to do this?

I have used the Picard tools SamToFastq but to my knowledge there is not a script in Picard or SamTools to do exactly what I described here (or maybe there is and I just haven't found it!).

Thank you!
pbm13 is offline   Reply With Quote
Old 06-28-2013, 03:54 PM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

If the reads are unpaired, how can you pull their mate?

Juts work out the read names you desire, and write a short script to fish those reads out.
swbarnes2 is offline   Reply With Quote
Old 06-28-2013, 10:02 PM   #3
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

^ good point...i have a feeling there is some miswording in the question.

I'd try to truncate the file through some form of filtering (maybe samtools or bamtools) and then use one of the sam/bam to fastq conversion scripts.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 06-29-2013, 05:59 PM   #4
pbm13
Junior Member
 
Location: NH

Join Date: Jan 2013
Posts: 5
Default

Yes, sorry for the unclear wording. The SAM file is a result of mapping paired end reads to a reference. I have a SAM file with mapped mated pairs that I was able to convert to a FASTQ which worked great. But I also have a SAM file with mapped unmated pairs - it is with this file that I would like to use to pull the reads that mapped (but their "mate" did not) and their pair from the original FASTQ files.

Ideally the output would be these pairs in a FASTQ file.
pbm13 is offline   Reply With Quote
Old 06-29-2013, 07:44 PM   #5
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

You should be able to extract those alignments as long as the aligner you used set the flags right. The unmapped mates will have 0x4 set and the mapped mates should have 0x8 set. You might need to name sort the bam first but then you could pull out only those reads with this:
Code:
samtools view -f 0xC -b alignments.bam > singletons.bam
Then convert that bam file into fastq.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Old 06-29-2013, 07:46 PM   #6
sdriscoll
I like code
 
Location: San Diego, CA, USA

Join Date: Sep 2009
Posts: 438
Default

Actually that will also pull out all unaligned reads in addition to your singleton alignments. So more filtering will be necessary. Pairs make this tricky because the SAM annotation of pairs is messy.
__________________
/* Shawn Driscoll, Gene Expression Laboratory, Pfaff
Salk Institute for Biological Studies, La Jolla, CA, USA */
sdriscoll is offline   Reply With Quote
Reply

Tags
fastq, paired end, sam

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 01:11 AM.


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