Problem with alignment

nayshool@gmail.com · 08-05-2014, 05:55 AM

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

dpryan · 08-05-2014, 06:05 AM

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 · 08-05-2014, 06:07 AM

BTW, please don't cross post on both here and biostars, it creates duplicate work from the community.

maubp · 08-05-2014, 07:21 AM

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/

nayshool@gmail.com · 08-05-2014, 08:22 AM

I'll try to remove the question from BioStars

Quote:

Originally Posted by dpryan

BTW, please don't cross post on both here and biostars, it creates duplicate work from the community.

nayshool@gmail.com · 08-05-2014, 08:25 AM

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

Quote:

Originally Posted by dpryan

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.

swbarnes2 · 08-05-2014, 03:29 PM

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

nayshool@gmail.com · 08-06-2014, 08:34 AM

Quote:

Originally Posted by swbarnes2

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!

swbarnes2 · 08-06-2014, 03:36 PM

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.

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

08-05-2014, 05:55 AM	#1
nayshool@gmail.com Junior Member Location: israel Join Date: Aug 2014 Posts: 5	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 06:04 AM.

08-05-2014, 06:05 AM	#2
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	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.

08-05-2014, 06:07 AM	#3
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	BTW, please don't cross post on both here and biostars, it creates duplicate work from the community.

08-05-2014, 07:21 AM	#4
maubp Peter (Biopython etc) Location: Dundee, Scotland, UK Join Date: Jul 2009 Posts: 1,543	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/

08-05-2014, 03:29 PM	#7
swbarnes2 Senior Member Location: San Diego Join Date: May 2008 Posts: 912	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

08-06-2014, 03:36 PM	#9
swbarnes2 Senior Member Location: San Diego Join Date: May 2008 Posts: 912	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.