SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Which bam file is what I need to count reads? wmseq Bioinformatics 15 11-05-2013 01:34 PM
BAM file containing ALL reads, next step? chipnewb Bioinformatics 0 01-31-2013 08:50 AM
How to exclude some id's from the file by grep or any other command M.Verma Bioinformatics 11 01-16-2013 07:40 PM
Why there is more reads in BAM file than I can see in tview? tomjan Bioinformatics 0 04-26-2012 01:41 PM
Getting Reads from Bam file empyrean Bioinformatics 4 10-12-2011 01:57 PM

Reply
 
Thread Tools
Old 01-25-2014, 09:23 AM   #1
papori
Senior Member
 
Location: berd

Join Date: Dec 2010
Posts: 179
Default grep only reads from bam file

Hi all,
I am working on TCGA WGS.. and the files are huge..
Is there a way to copy all the reads directly from the bam file without using samtools and bam2fastx (from tophat)?

Right now i am using:
samtools sort -n *.bam *.bam.sorted
bam2fastx -q -Q -A -P -o out.fq *.bam.sorted.bam

and it is very slow..
i dont mind to stay in binary mode, and i am not using the quality + the order isnt relevant..

Or,
to spped up what i am doing now..

Thanks!
papori is offline   Reply With Quote
Old 01-25-2014, 02:19 PM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,479
Default

Are the reads paired-end and do you care if the output is? If the answer to either of those is "no", then you can skip the sort step. If the answer to both is "yes", then there's really no great way to do that. Realistically, paired-end reads will normally close to each other in a file (since that's how they'll map), so you could just keep them in memory until running into the appropriate mate. You might see if you can just get the fastq files, it might be faster.
dpryan is offline   Reply With Quote
Old 01-26-2014, 04:17 AM   #3
papori
Senior Member
 
Location: berd

Join Date: Dec 2010
Posts: 179
Default

Quote:
Originally Posted by dpryan View Post
Are the reads paired-end and do you care if the output is? If the answer to either of those is "no", then you can skip the sort step. If the answer to both is "yes", then there's really no great way to do that. Realistically, paired-end reads will normally close to each other in a file (since that's how they'll map), so you could just keep them in memory until running into the appropriate mate. You might see if you can just get the fastq files, it might be faster.
If i am interesting only in the first pair, can i use bam2fastx without sorting?
Or, what you meant is that there are no seperator between the 'paired ends', and i will have sinfle file with the 2 paired?

thanks
papori is offline   Reply With Quote
Old 01-26-2014, 06:32 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,479
Default

If you're only interested in the first mate in a pair, then:

Code:
samtools view -f 0x40 sample.bam | awk '{printf("%s\n%s\n+\n%s\n", $1, $10, $11)}' | gzip > sample.fq.gz
will work and probably be a bit faster. If you removed the "-f 0x40", then you'd have both mates of the pair in the same file, though not necessarily next to each other. Of course, if the file contains non-linear alignments (i.e., a single read can be split across multiple lines, using the 0x800 bit in the flag) then this won't work.
dpryan is offline   Reply With Quote
Old 01-27-2014, 12:39 AM   #5
papori
Senior Member
 
Location: berd

Join Date: Dec 2010
Posts: 179
Default

Quote:
Originally Posted by dpryan View Post
If you're only interested in the first mate in a pair, then:

Code:
samtools view -f 0x40 sample.bam | awk '{printf("%s\n%s\n+\n%s\n", $1, $10, $11)}' | gzip > sample.fq.gz
will work and probably be a bit faster. If you removed the "-f 0x40", then you'd have both mates of the pair in the same file, though not necessarily next to each other. Of course, if the file contains non-linear alignments (i.e., a single read can be split across multiple lines, using the 0x800 bit in the flag) then this won't work.
I tried this and i got that the first end is: 130713161918
But when i used:
samtools sort -n *.bam *.bam.sorted
bam2fastx -q -Q -A -P -o out.fq *.bam.sorted.bam
i got: 146521575509

What could be the reason?
The procedure that you offer was much quicker.. but it might miss some reads?
papori is offline   Reply With Quote
Old 01-27-2014, 01:47 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,479
Default

The method I gave will never miss any reads. The difference could only be due to the presence of the "-f 0x40" in the command that I suggested and the lack of that in bam2fastx, which put both reads of a pair in the output. I should note that one down-side to the method I gave is that the reads won't be reverse-complemented as appropriate prior to output. If you don't care about the order, then just skip the sort step and use bam2fastx.
dpryan is offline   Reply With Quote
Old 01-27-2014, 03:32 AM   #7
papori
Senior Member
 
Location: berd

Join Date: Dec 2010
Posts: 179
Default

Quote:
Originally Posted by dpryan View Post
The method I gave will never miss any reads. The difference could only be due to the presence of the "-f 0x40" in the command that I suggested and the lack of that in bam2fastx, which put both reads of a pair in the output. I should note that one down-side to the method I gave is that the reads won't be reverse-complemented as appropriate prior to output. If you don't care about the order, then just skip the sort step and use bam2fastx.
I am interesting in all the reads that are in the bam file, without missing any read.
When i am using bam2fastx without sorting i am getting this error:
Error: couldn't retrieve both reads for pair HWI-ST512_115262631:6:1306:12150:22582. Perhaps the input file is not sorted by name?
(using 'samtools sort -n' might fix this)
papori is offline   Reply With Quote
Old 01-27-2014, 03:57 AM   #8
gt1
Junior Member
 
Location: Cambridge, UK

Join Date: Jul 2013
Posts: 9
Default

Quote:
Originally Posted by papori View Post
I am interesting in all the reads that are in the bam file, without missing any read.
When i am using bam2fastx without sorting i am getting this error:
Error: couldn't retrieve both reads for pair HWI-ST512_115262631:6:1306:12150:22582. Perhaps the input file is not sorted by name?
(using 'samtools sort -n' might fix this)
Apparently bam2fastx requires the input file to be precollated (have pairs next to each other). There are several other programs that perform this collation process on the fly. One of them is bamtofastq in biobambam. Another is SamToFastq in Picard.
gt1 is offline   Reply With Quote
Old 01-27-2014, 04:23 AM   #9
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,800
Default

@papori: Have you checked to see if the original fastq files are available for TCGA WGS data?
GenoMax is offline   Reply With Quote
Old 01-27-2014, 04:28 AM   #10
papori
Senior Member
 
Location: berd

Join Date: Dec 2010
Posts: 179
Default

Quote:
Originally Posted by GenoMax View Post
@papori: Have you checked to see if the original fastq files are available for TCGA WGS data?
As i understood from the FAQ of TCGA only fastq of RNA/WXS available now..
Do you know somthing else?
papori is offline   Reply With Quote
Old 01-27-2014, 04:36 AM   #11
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,800
Default

Quote:
Originally Posted by papori View Post
As i understood from the FAQ of TCGA only fastq of RNA/WXS available now..
Do you know somthing else?
No I don't. I was wondering if the group that submitted the WGS data also had made the sequences available. Have you tried to ask for sequence data to see if it can be made available?
GenoMax is offline   Reply With Quote
Reply

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 08:05 AM.


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