SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BSmooth alignment problem genefan Epigenetics 9 11-18-2013 11:07 PM
SOAP2 alignment problem logicthief Bioinformatics 0 04-22-2013 12:40 PM
Problem with alignment with pseudogenes alexanderxcy Bioinformatics 1 10-02-2012 03:37 PM
Blast problem in metagenomics alignment rexxi Bioinformatics 3 07-10-2012 06:58 AM
Gene Sequencing Alignment Problem anonya Bioinformatics 0 05-01-2012 11:59 PM

Reply
 
Thread Tools
Old 08-05-2014, 04:55 AM   #1
nayshool@gmail.com
Junior Member
 
Location: israel

Join Date: Aug 2014
Posts: 5
Default Problem with alignment

Hello everyone!

We have received from BGI china a few samples we sent for exome sequence + Basic analysis. We got the fastQ files and BAM files. I have made my own BAM files using BWA and Picard tools using the following commends:

BWA:
bwa mem human_g1k_v37.fasta file1.fq.gz file2.fq.gz file3.fq.gz file4.fq.gz > result.sam

PICARD:
java -jar SortSam.jar I=result.sam O=result.bam SORT_ORDER=coordinate

The size of the files that I got (2Gb) , was about a half from the BAM files BGI (4.2 Gb) sent me.
Here are the samtools flagstat results for both files:

my BAM file:

[gen-biorep@gen-biorep-3 VA6]$ samtools flagstat result.bam
28842836 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
28824861 + 0 mapped (99.94%:-nan%)
28842836 + 0 paired in sequencing
14423043 + 0 read1
14419793 + 0 read2
28593748 + 0 properly paired (99.14%:-nan%)
28814004 + 0 with itself and mate mapped
10857 + 0 singletons (0.04%:-nan%)
165631 + 0 with mate mapped to a different chr
131466 + 0 with mate mapped to a different chr (mapQ>=5)



BGI BAM File:

[gen-biorep@gen-biorep-3 result_alignment]$ samtools flagstat VA6.rmdup.bam
57668392 + 0 in total (QC-passed reads + QC-failed reads)
3039053 + 0 duplicates
57369744 + 0 mapped (99.48%:-nan%)
57668392 + 0 paired in sequencing
28834196 + 0 read1
28834196 + 0 read2
56991680 + 0 properly paired (98.83%:-nan%)
57273460 + 0 with itself and mate mapped
96284 + 0 singletons (0.17%:-nan%)
198072 + 0 with mate mapped to a different chr
181774 + 0 with mate mapped to a different chr (mapQ>=5


why does half of my reads disappear? How can I solve this problem?

Thank you,

Omri

Last edited by nayshool@gmail.com; 08-05-2014 at 05:04 AM.
nayshool@gmail.com is offline   Reply With Quote
Old 08-05-2014, 05:05 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Unless you mislabeled the flagstat output, it looks like your results contain about twice as many alignments as those from BGI, though that could be due to yours containing multiple copies of multimappers while BGI's doesn't. In general, you'd need to look at how they aligned things and that will likely tell you why there's a difference.
dpryan is offline   Reply With Quote
Old 08-05-2014, 05:07 AM   #3
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

BTW, please don't cross post on both here and biostars, it creates duplicate work from the community.
dpryan is offline   Reply With Quote
Old 08-05-2014, 06:21 AM   #4
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

Does either set of BAM file include unmapped read pairs? 'samtools idxstats example.bam' might help answer this.

P.S. Duplicate question was https://www.biostars.org/p/108537/
maubp is offline   Reply With Quote
Old 08-05-2014, 07:22 AM   #5
nayshool@gmail.com
Junior Member
 
Location: israel

Join Date: Aug 2014
Posts: 5
Default Sorry for the duplicate, won't happen again..

I'll try to remove the question from BioStars

Quote:
Originally Posted by dpryan View Post
BTW, please don't cross post on both here and biostars, it creates duplicate work from the community.
nayshool@gmail.com is offline   Reply With Quote
Old 08-05-2014, 07:25 AM   #6
nayshool@gmail.com
Junior Member
 
Location: israel

Join Date: Aug 2014
Posts: 5
Default

thank you. I have mislabled the flagstat ouputs, I have change it. BGI have more reads.

Quote:
Originally Posted by dpryan View Post
Unless you mislabeled the flagstat output, it looks like your results contain about twice as many alignments as those from BGI, though that could be due to yours containing multiple copies of multimappers while BGI's doesn't. In general, you'd need to look at how they aligned things and that will likely tell you why there's a difference.
nayshool@gmail.com is offline   Reply With Quote
Old 08-05-2014, 02:29 PM   #7
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

The .bam from bgi reports the exact same number of read 1 reads as read 2 reads. Yours does not. Are you sure that giving it 4 fastq files at once is the proper way to use bwa mem? I wonder if it is just ignoring those last two files. Maybe you could spot-check to see if read names from those fastqs are in your .bam
swbarnes2 is offline   Reply With Quote
Old 08-06-2014, 07:34 AM   #8
nayshool@gmail.com
Junior Member
 
Location: israel

Join Date: Aug 2014
Posts: 5
Thumbs up problem solved

Quote:
Originally Posted by swbarnes2 View Post
The .bam from bgi reports the exact same number of read 1 reads as read 2 reads. Yours does not. Are you sure that giving it 4 fastq files at once is the proper way to use bwa mem? I wonder if it is just ignoring those last two files. Maybe you could spot-check to see if read names from those fastqs are in your .bam
Well i merged the 4 fastQ into one file and it did the trick. It's appear that BWA-MEM ignored from the last 2 files as you said. I got the right size of BAMs after that . Thank you swbarnes2 very much for your help!
nayshool@gmail.com is offline   Reply With Quote
Old 08-06-2014, 02:36 PM   #9
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Well, don't merge them into 1 file, that's just throwing away real paired information.

You want all of the read 1 data in one fastq, and all the read 2 data in a second. Just make sure that the nth entry in the first fastq has the same read name as the nth read in file 2.
swbarnes2 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 03:27 PM.


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