SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to extract read count from bam file lylyf1987 RNA Sequencing 0 04-07-2015 11:19 AM
fastq.gz stats READ-COUNT BASE-COUNT jgibbons1 Bioinformatics 9 10-30-2013 05:24 AM
convert sorted bam to sorted sam for htseq-count ugolino Bioinformatics 3 10-19-2012 06:30 AM
multiBamCov or htseq-count to count read per feature ? NicoBxl Bioinformatics 1 07-03-2012 02:05 AM
small rna seq with total lane read count different between conditions anle Bioinformatics 0 01-25-2012 01:25 PM

Reply
 
Thread Tools
Old 05-12-2015, 08:05 AM   #1
pkMyt1
Junior Member
 
Location: Chapel Hill, NC U.S.A

Join Date: Dec 2013
Posts: 4
Default Sorted BAM read count >2x total FASTQ count

This might be a silly question but it is bugging me and I need the answer. I have a set of paired end FASTQ files that contain about 30 million reads total. After aligning them with BWA and sorting the output with samtools, the resulting BAM file now has about 72 million reads. Why????
pkMyt1 is offline   Reply With Quote
Old 05-12-2015, 08:14 AM   #2
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 701
Default

Show us your alignment and sorting commands.
Richard Finney is offline   Reply With Quote
Old 05-12-2015, 08:28 AM   #3
fanli
Senior Member
 
Location: California

Join Date: Jul 2014
Posts: 198
Default

Multiple alignments...
fanli is offline   Reply With Quote
Old 05-12-2015, 08:33 AM   #4
pkMyt1
Junior Member
 
Location: Chapel Hill, NC U.S.A

Join Date: Dec 2013
Posts: 4
Default

Both come from iterating the files in Python. The FASTQ files are read in in blocks of four lines each which is one read. This example is a MiSeq run so 30 million (15 million in each FASTQ) seems realistic. The BAM count is from

bamfile = pysam.AlignmentFile(o['bamfile'], "rb")
bamfile.count()
or
bamfile_reads = functools.reduce(lambda x, y: x + y, [eval('+'.join(l.rstrip('\n').split('\t')[2:])) for l in pysam.idxstats(o['bamfile'])])

https://www.biostars.org/p/1890/

or simply counting the reads as I iterate the BAM file to do my analysis.

Last edited by pkMyt1; 05-12-2015 at 08:36 AM.
pkMyt1 is offline   Reply With Quote
Old 05-12-2015, 08:39 AM   #5
pkMyt1
Junior Member
 
Location: Chapel Hill, NC U.S.A

Join Date: Dec 2013
Posts: 4
Default

Quote:
Originally Posted by fanli View Post
Multiple alignments...
So....
Would this imply my alignment settings are keeping things I should not?

bwa mem -a -T 25 -L '(100, 100)'
pkMyt1 is offline   Reply With Quote
Old 05-12-2015, 09:47 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

I would imagine that the -a flag is to blame.
dpryan is offline   Reply With Quote
Old 05-12-2015, 11:35 AM   #7
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by pkMyt1 View Post
So....
Would this imply my alignment settings are keeping things I should not?
That depends on the goal of your experiment. What are you trying to do?
Brian Bushnell is offline   Reply With Quote
Old 05-13-2015, 04:43 AM   #8
pkMyt1
Junior Member
 
Location: Chapel Hill, NC U.S.A

Join Date: Dec 2013
Posts: 4
Default

Quote:
Originally Posted by Brian Bushnell View Post
That depends on the goal of your experiment. What are you trying to do?
This is duplex exome sequencing. Very deep but only about 80 kb of capture. I did not want to lose any alignments where one read aligned and the other did not either due to a translocation or simply a sequencing error. This is why I did the -a option. Each read is uniquely tagged so I had been able to filter things in the end. This is the first time I have seen this but it is also the first time I have run a sample that I know contains many chromosomal rearrangements in the way of translocations, duplications, deletions. I will need to try and pull out some of these multiple alignments and have a look at them so I can understand what they are better.
pkMyt1 is offline   Reply With Quote
Old 05-13-2015, 09:43 AM   #9
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

In that case, it sounds like considering all good alignments of the reads is probably best. The reason for all the multiple alignments is presumably that you're targeting a repetitive region.
Brian Bushnell is offline   Reply With Quote
Reply

Tags
bam files, fastq files

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 03:32 PM.


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