SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
'Properly paired' reads in sam flag from TopHat mapping AdamB Bioinformatics 9 03-08-2012 07:30 AM
How to align SOLiD data? Heisman SOLiD 11 01-29-2012 04:21 PM
how do I output the CS tag for BWA align of SOLID reads? KevinLam Bioinformatics 16 07-23-2011 10:06 PM
How to check if reads are properly paired in mate-pair data? genepool_bee Bioinformatics 2 02-22-2011 01:07 AM
using BWA to align SOLiD fastq files from 1000 Genomes tgenahmet Bioinformatics 1 10-15-2009 05:19 PM

Reply
 
Thread Tools
Old 06-23-2011, 11:46 PM   #1
alig
Member
 
Location: adelaide

Join Date: Sep 2008
Posts: 43
Unhappy BWA align SOLiD PE data gives poor mapping & 0.5% properly paired

Hi,

This is my first attempt at aligning SOLiD paired end sequence which is whole exome data.
F3 file has reads minlength=50, F5-BC file has reads minlength=35

Previously I have only worked with Illumina sequence.

The following are the commands I used to do an assembly using BWA.

Unfortunately as you will see the results from samtools flagstat are disappointing. Only 57% reads mapped, & worse still 0.5% properly paired.

Can anyone tell me where I'm going wrong?
Thank you
alig

PS I've played with the bwa aln parameters a bit with no improvement. I'm wondering whether it's the solid2fastq step or even my reference sequence. Any help will be greatly appreciated.

#Created colorspace index
bwa index -a bwtsw -c hg18_noran_mask.fa

#modified solid2fastq.pl script

95 #if (/^>(\d+)_(\d+)_(\d+)_[FR]3/) {
96 if (/^>(\d+)_(\d+)_(\d+)_[F3|R3|F5-P2]/) {

#also renamed the F5-P2 to R3:

solid_data_F5-P2.csfasta -> solid_data_R3.csfasta
solid_data_F5-P2_QV.qual -> solid_data_R3_QV.qual

#converted to fastq.gz
perl ../../BWA_new/bwa-0.5.9/solid2fastq_mod.pl solid_data solid_data

# bwa aln

./../BWA_new/bwa-0.5.9/bwa aln -n 0.04 -o 1 -e -1 -d 5 -i 5 -l 25 -k 2 -t 4 -M 3 -O 11 -E 4 -c -q 0 ../cs_index/hg18_noran_mask.fa solid_data.read1.fastq.gz > solid_data_read1_bwa_hg18.sai

Repeated for fastq 2 file for paired end reads

#bwa sampe

../../BWA_new/bwa-0.5.9/bwa sampe -n 3 -N 10 ../cs_index/hg18_noran_mask.fa solid_data_read1_bwa_hg18.sai solid_data_read2_bwa_hg18.sai solid_data.read1.fastq.gz solid_data.read2.fastq.gz > solid_data.sam

#samtools

samtools view -bT ../cs_index/hg18_noran_mask.fa solid_data.sam -o solid_data.bam (as no .fai file from indexed RefSeq)

samtools sort solid_data.bam solid_data_sorted

samtools index solid_data_sorted.bam

samtools flagstat solid_data.bam
235101678 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
134185711 mapped (57.08%:- nan%)
235101678 paired in sequencing
117550839 read 1
117550839 read 2
1196926 properly paired (0.51%:- nan%)
89122266 with itself & mate mapped
45063445 singletons (19.1%:- nan%)
2798688 with mate mapped to different chr
945434 with mate mapped to a different chr (mapQ>=5)
alig is offline   Reply With Quote
Old 06-24-2011, 03:35 AM   #2
poisson200
Member
 
Location: united kingdom

Join Date: Feb 2010
Posts: 63
Default

Hello there,
I had a problem mapping paired end data with bowtie. In the bowtie manual it is thought that the paired ends from colourspace are in the forward strand- forward strand --ff orientation most of the time (well it maybe) but my reads F3 and F5-BC are in --fr. I did not realise this at first and my read mapping with bowtie was first at 0.1% with --ff, which rose to ~20% with --fr -X 400.

Obviously, this does not look like your problem but as a remark, my mapping is only 20%, so 80% not mapped as pairs. I think SOLID does the blunder bus approach and has a lot of collateral reads that are junk.

You could try bioscope mapping and bowtie, plus Sasson and michael pre-filter
poisson200 is offline   Reply With Quote
Old 06-24-2011, 03:30 PM   #3
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 700
Default

57% is okay^h^h^h^h great!

Last batch the wet lab guys sent me averaged 40%.

Resulting bam files were still quite usable.

BWA [ using default params except -c, of course ] aligns much less reads than the bioscope pipeline. The BWA looks much cleaner as bioscope aggressively hard clips reads and goes ahead and aligns them, many of which appear to be wishful thinking. [ This is using the default bioscope parameters available with the "demo" on NIH's biowulf].
Richard Finney is offline   Reply With Quote
Old 07-08-2011, 08:44 AM   #4
ryuky
Junior Member
 
Location: china

Join Date: Mar 2011
Posts: 7
Default

you should change the input position, just like this :
bwa sampe -n 3 -N 10 ../cs_index/hg18_noran_mask.fa solid_data_read2_bwa_hg18.sai solid_data_read1_bwa_hg18.sai solid_data.read2.fastq.gz solid_data.read1.fastq.gz>*.bam
I can't figure out why....
ryuky 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 07:57 AM.


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