SEQanswers

Go Back   SEQanswers > Sequencing Technologies/Companies > SOLiD



Similar Threads
Thread Thread Starter Forum Replies Last Post
Strange Bioanalyzer results wzombie RNA Sequencing 5 01-11-2012 03:24 AM
strange results of samtools liying Bioinformatics 3 09-22-2011 11:02 PM
BWA align SOLiD PE data gives poor mapping & 0.5% properly paired alig Bioinformatics 3 07-08-2011 08:44 AM
Strange DE results SMcTaggart Bioinformatics 0 11-25-2010 04:53 AM
Strange Results After maq2sam-long wjeck Bioinformatics 4 08-13-2009 10:58 AM

Reply
 
Thread Tools
Old 08-01-2010, 05:55 AM   #1
Hit
Member
 
Location: United States, TN

Join Date: Jul 2010
Posts: 15
Thumbs up strange mapping results bwa + SOLiD

Hi,
I used bwa to map a pair-end SOLiD data. samtools flagstat showed 66% mapped reads, but only 0.17% properlly mapped.
Here are some detail and my cmd lines:
my input file: SOLID_F3.csfasta + SOLID_F3_QV.qual
SOLID_R3.csfasta + SOLID_R3_QV.qual
(1) make color-space reference: bwa index -p hg18_cs.fa -a bwtsw -c hg18.fa
(2) convert solid to fastq: perl solid2fastq.pl SOLID_ test
this resulted test.read1.fastq.gz test.read2.fastq.gz test.single.fastq.gz
unzip them
(3) bwa aln -t 4 -c hg18_cs.fa test.read1.fastq > test_aln_sa1.sai
bwa aln -t 4 -c hg18_cs.fa test.read2.fastq > test_aln_sa2.sai
(4) bwa sampe hg18_cs.fa test_aln_sa1.sai test_aln_sa2.sai test.read1.fastq test.read2.fastq > aln.sam
(5) samtools view -bt hg18_cs.fa -o aln.bam aln.sam
[samopen] SAM header is present: 25 sequences.

(6) $ samtools flagstat aln.bam
The resutl shows mapped (66.64%), properly paired (0.15%)

Can anyone tell me where I did wrong?

Thanks!
Hit is offline   Reply With Quote
Old 08-01-2010, 10:33 AM   #2
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Are you using the paired end or mate pair protocol? Was the insert size estimate from BWA around what you expect?
nilshomer is offline   Reply With Quote
Old 08-01-2010, 11:02 AM   #3
Hit
Member
 
Location: United States, TN

Join Date: Jul 2010
Posts: 15
Default

Hi Nil,

The data is paired end, one is 50bp while the other is 25bp. The original file name is actually F3 and F5-P2: I changed F5-P2 to R3 to make it recognizable by solid2fastq.pl provided in bwa. The insert size looks fine for me. Here are some example lines (when using ):
[infer_isize] (25, 50, 75) percentile: (150, 160, 173)
[infer_isize] low and high boundaries: 104 and 219 for estimating avg and std
[infer_isize] inferred external isize from 100442 pairs: 161.585 +/- 17.497
[infer_isize] skewness: 0.147; kurtosis: 0.337; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 285 (7.05 sigma)
[bwa_sai2sam_pe_core] time elapses: 22.68 sec
[bwa_sai2sam_pe_core] changing coordinates of 73 alignments.
[bwa_sai2sam_pe_core] align unmapped mate...
[bwa_paired_sw] 20 out of 79202 Q17 singletons are mated.
[bwa_paired_sw] 316 out of 125430 Q17 discordant pairs are fixed.
[bwa_sai2sam_pe_core] time elapses: 14.42 sec
[bwa_sai2sam_pe_core] refine gapped alignments... 4.12 sec
[bwa_sai2sam_pe_core] print alignments... 1.31 sec
[bwa_sai2sam_pe_core] 89915392 sequences have been processed.
[bwa_sai2sam_pe_core] convert to sequence coordinate...
[infer_isize] (25, 50, 75) percentile: (150, 160, 173)
[infer_isize] low and high boundaries: 104 and 219 for estimating avg and std
[infer_isize] inferred external isize from 97115 pairs: 161.534 +/- 17.614
[infer_isize] skewness: 0.137; kurtosis: 0.339; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 286 (7.05 sigma)
Hit is offline   Reply With Quote
Old 08-01-2010, 11:41 AM   #4
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by Hit View Post
Hi Nil,

The data is paired end, one is 50bp while the other is 25bp. The original file name is actually F3 and F5-P2: I changed F5-P2 to R3 to make it recognizable by solid2fastq.pl provided in bwa. The insert size looks fine for me. Here are some example lines (when using ):
[infer_isize] (25, 50, 75) percentile: (150, 160, 173)
[infer_isize] low and high boundaries: 104 and 219 for estimating avg and std
[infer_isize] inferred external isize from 100442 pairs: 161.585 +/- 17.497
[infer_isize] skewness: 0.147; kurtosis: 0.337; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 285 (7.05 sigma)
[bwa_sai2sam_pe_core] time elapses: 22.68 sec
[bwa_sai2sam_pe_core] changing coordinates of 73 alignments.
[bwa_sai2sam_pe_core] align unmapped mate...
[bwa_paired_sw] 20 out of 79202 Q17 singletons are mated.
[bwa_paired_sw] 316 out of 125430 Q17 discordant pairs are fixed.
[bwa_sai2sam_pe_core] time elapses: 14.42 sec
[bwa_sai2sam_pe_core] refine gapped alignments... 4.12 sec
[bwa_sai2sam_pe_core] print alignments... 1.31 sec
[bwa_sai2sam_pe_core] 89915392 sequences have been processed.
[bwa_sai2sam_pe_core] convert to sequence coordinate...
[infer_isize] (25, 50, 75) percentile: (150, 160, 173)
[infer_isize] low and high boundaries: 104 and 219 for estimating avg and std
[infer_isize] inferred external isize from 97115 pairs: 161.534 +/- 17.614
[infer_isize] skewness: 0.137; kurtosis: 0.339; ap_prior: 1.00e-05
[infer_isize] inferred maximum insert size: 286 (7.05 sigma)
I am not sure if BWA supports the SOLiD paired end protocol, since it assumes that the two ends are on opposite strands, which is not the case for paired end. That is why it is getting upset. You could email the BWA mailing list, PM the user "lh3" (the author of BWA), or how "lh3" reads this.
nilshomer is offline   Reply With Quote
Old 08-01-2010, 04:20 PM   #5
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Check out what Nils mentioned to confirm if that's the reason. Let us know what
you find.

Also, I'd like to mention that BFast has a development branch called bfast2. It
allows you to perform the mapping step (match) using bwa's algorithm. After
that you'd perform the normal bfast process (localalign, postprocess).

In addition, bfast2 has also a hybrid mode that allows you to perform matches
using different algorithms on each end of your reads. This is particularly
useful for PE(50+25) since bfast match does not perform well with 25bp reads.
__________________
-drd
drio is offline   Reply With Quote
Old 08-01-2010, 06:19 PM   #6
Hit
Member
 
Location: United States, TN

Join Date: Jul 2010
Posts: 15
Default

Thank you, Nils and Drio.

I searched online and got the following information:

As explained in this file:
http://www3.appliedbiosystems.com/cm...cms_080186.pdf

my case should be the second one, where F3 tag is forward and F5-P2 tag is reverse.

However, I saw this post in the BWA mailing list:
http://sourceforge.net/mailarchive/f...e=bio-bwa-help

In Heng's reply, it's said "For Illumina reads, only Forward-reverse pairs are set bit 2; for SOLiD reads, R3-forward-F3-forward or F3-reverse-R3-reverse pairs are set bit 2."

So I'm guessing my case is: it is solid data but it's forward-reverse, so bwa correctly mapped most of them (66%) but failed to recognize the reads as perporly paired -- this is just my guessing. I would greatly appreciate if anyone can help me figure this out.
Hit is offline   Reply With Quote
Old 08-01-2010, 07:46 PM   #7
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by Hit View Post
Thank you, Nils and Drio.

I searched online and got the following information:

As explained in this file:
http://www3.appliedbiosystems.com/cm...cms_080186.pdf

my case should be the second one, where F3 tag is forward and F5-P2 tag is reverse.

However, I saw this post in the BWA mailing list:
http://sourceforge.net/mailarchive/f...e=bio-bwa-help

In Heng's reply, it's said "For Illumina reads, only Forward-reverse pairs are set bit 2; for SOLiD reads, R3-forward-F3-forward or F3-reverse-R3-reverse pairs are set bit 2."

So I'm guessing my case is: it is solid data but it's forward-reverse, so bwa correctly mapped most of them (66%) but failed to recognize the reads as perporly paired -- this is just my guessing. I would greatly appreciate if anyone can help me figure this out.
Note that the read-rescue code in "bwa sampe" depends on the correct orientation, so it is not just the proper-pair bit being different.
nilshomer is offline   Reply With Quote
Old 08-11-2010, 02:26 AM   #8
tamara_0021
Junior Member
 
Location: UK

Join Date: Aug 2010
Posts: 2
Default

Hi All,
Can anybody tell me how to run the script solid2fastq.pl.

I have a cs.fasta file, I have to calls the SNPs and want to use MAQ for that. However I need the files in fastq. Using BWA I did get a colour space index. But after that if i try to run
perl solid2fastq.pl file.sc.fasta out.cs.fa it is not working
the cs.fasta is a F3 file.
It tell me it cannot open the file.

Many thanks
tamara_0021 is offline   Reply With Quote
Old 08-11-2010, 04:36 AM   #9
Hit
Member
 
Location: United States, TN

Join Date: Jul 2010
Posts: 15
Default

Quote:
Originally Posted by tamara_0021 View Post
Hi All,
Can anybody tell me how to run the script solid2fastq.pl.

I have a cs.fasta file, I have to calls the SNPs and want to use MAQ for that. However I need the files in fastq. Using BWA I did get a colour space index. But after that if i try to run
perl solid2fastq.pl file.sc.fasta out.cs.fa it is not working
the cs.fasta is a F3 file.
It tell me it cannot open the file.

Many thanks
Check it here:
http://seqanswers.com/forums/showthread.php?t=5909
drio's answer.
Hit is offline   Reply With Quote
Old 08-11-2010, 05:20 AM   #10
tamara_0021
Junior Member
 
Location: UK

Join Date: Aug 2010
Posts: 2
Default

ok, thanks, I will try.... but tell me what does this instruction means
<in.title> is the string showed in the `# Title:' line of a
".csfasta" read file. Then <in.title>F3.csfasta is read sequence
file and <in.title>F3_QV.qual is the quality file. If
<in.title>R3.csfasta is present, this script assumes reads are
paired; otherwise reads will be regarded as single-end.

The read name will be <out.prefix>anel_x_y/[12] with `1' for R3
tag and `2' for F3. Usually you may want to use short <out.prefix>


I want to know what is it actually asking ?

Many Thanks for the link
tamara_0021 is offline   Reply With Quote
Old 09-24-2010, 04:37 AM   #11
Brugger
Member
 
Location: Cambridge, UK

Join Date: Mar 2010
Posts: 21
Default

I have posted a patch plus an update perl script to the bwa malinglist so bwa now can handle PE data as well as MP.
Brugger is offline   Reply With Quote
Old 05-09-2011, 10:54 AM   #12
2007lab
Member
 
Location: NYC

Join Date: Mar 2009
Posts: 14
Default

Hi Brugger,
I saw your mail in the bwa mailing list about your patch.
Do you know if the later version of BWA incorporates this?
What I have installed is 0.5.9-r16
2007lab 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 12:01 AM.


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