SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Bowtie Paired end 100% failed to align rahul.m.dhodapkar Bioinformatics 15 09-12-2013 02:05 PM
Bowtie and reads that failed to align: (100.00%) michy Bioinformatics 7 02-08-2011 06:42 PM
Discarding all reads (and their mates) that fail Picard validation wdt Bioinformatics 0 12-08-2010 04:22 PM
Bowtie to align reads to single chromosome or region? jjw14 Bioinformatics 1 10-07-2010 06:20 PM
bowtie and solexa paired ends wont align michy Bioinformatics 1 09-02-2010 05:46 AM

Reply
 
Thread Tools
Old 05-12-2012, 02:07 PM   #1
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default Bowtie: Reads fail to align

I am attempting to align a transcriptome sequenced with Hiseq to a reference using bowtie. The parameters I am using are:

bowtie -S -p 2 reference -q --phred64-quals

And none of the reads align. They also do not align if I do not include the quality parameter, or any modification such as the --ff suggested in related postings.

I have checked for adapter contamination and found very little, the reads were cleaned using ngs backbone, although I am not using that pipeline for anything downstream of cleaning. It has also been reported by some that they do not align in paired end but do in single end, my reads do not align in either case. Around 30% of the reas align in bwa.

Does anyone have any idea why this is the case?

Thanks!
Sarah
sasignor is offline   Reply With Quote
Old 05-14-2012, 08:29 AM   #2
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

The first ten lines of one of the reads files looks like this, if that helps:

@HWI-ST611_0200:5:1101:1394:2137#0/1
AATATTTCGATTCGCTATATTCCTCAAATCAAGGTTACAATTTATCAGTCCGGACCAACAAAGCCAAGATCAAGCGTGACACGGGCACACATATGGCAGC
+
Z^_ccccc\\S`cfgfa_daK``K[Q^bX[d[bRI^^XYSbBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@HWI-ST611_0200:5:1101:1465:2147#0/1
AATTTTTTCTTAGGTTTGGTCTGCACCACAGTTCGCCTCCAAAATTTTTGATTTTGCGGCCCAACGACTTTTTTTTTAAAGCACGCCTCCAACACAAGCT
+
___cccc]ecgc[Yb[bPbY[^d[bJP^b^^bc]Z_c`_BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@HWI-ST611_0200:5:1101:1348:2147#0/1
AAGTCTTATTGTTATAGGGACTTCACCTGATTTACTCACTTCGAATGAATCAACACATCACAACGCTCACGCCTGCGGCGCAACGTCAACGCGGAGACGC
sasignor is offline   Reply With Quote
Old 05-14-2012, 09:08 AM   #3
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Giving us the reads doesn't help much since without the reference we can not tell if they should map. What would be helpful is the full and exact command line you are using. I would expect a line similar to

bowtie -B -p 2 reference read.file

Also helpful would be the command line to BWT that generates the successful hits.
westerman is offline   Reply With Quote
Old 05-14-2012, 09:17 AM   #4
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

This is the full command in bowtie:

./bowtie -S -p 2 eletrans -q --phred64-quals -1 lb_gun1.pl_illumina.sm_gun1.sfastq -2 lb_gun2.pl_illumina.sm_gun2.sfastq gun2elebowtie.sam

The bwa command is:
./bwa bwasw -t 2 eletrans lb_gun1.pl_illumina.sm_gun1.sfastq > gun2elebwa.sam

(single end alignments in bwa, clearly, however single end does not improve alignments for bowtie)
sasignor is offline   Reply With Quote
Old 05-14-2012, 09:27 AM   #5
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

No error messages?

I always put my options (e.g., -S, -p, -q, --phred64-quals) before the reference. But I am not 100% sure that is required.

How about a 'ls -l' of the *.ebwt files just to make sure that the reference is properly formatted.
westerman is offline   Reply With Quote
Old 05-14-2012, 09:31 AM   #6
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

No error messages. The output looks like this:

# reads processed: 73455581
# reads with at least one reported alignment: 0 (0.00%)
# reads that failed to align: 73455581 (100.00%)
No alignments

Here are the reference files:

-rw-r--r-- 1 sarah staff 7365819 May 10 10:24 eletran.1.ebwt
-rw-r--r-- 1 sarah staff 1312752 May 10 10:24 eletran.2.ebwt
-rw-r--r-- 1 sarah staff 37151 May 10 10:23 eletran.3.ebwt
-rw-r--r-- 1 sarah staff 2625496 May 10 10:23 eletran.4.ebwt
-rw-r--r-- 1 sarah staff 7365819 May 10 10:24 eletran.rev.1.ebwt
-rw-r--r-- 1 sarah staff 1312752 May 10 10:24 eletran.rev.2.ebwt
sasignor is offline   Reply With Quote
Old 05-14-2012, 09:35 AM   #7
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

The fact that the command says 'eletrans' and the reference listed here is 'eletran' is not significant, I had named it differently at one point.
sasignor is offline   Reply With Quote
Old 05-14-2012, 09:36 AM   #8
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Two suggestions.

1) What happens if you just treat your reads as single-end; e.g.,

Code:
./bowtie -S -p 2 -q --phred64-quals eletrans  lb_gun1.pl_illumina.sm_gun1.sfastq,lb_gun2.pl_illumina.sm_gun2.sfastq gun2elebowtie.sam
2) What happens if you use the '--verbose' switch?
westerman is offline   Reply With Quote
Old 05-14-2012, 09:38 AM   #9
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

This is what happens if you do a single end alignment:

./bowtie -S -p 2 eletrans -q --phred64-quals lb_gun1.pl_illumina.sm_gun1.sfastq lb1gun2elebowtie4.sam
# reads processed: 73503003
# reads with at least one reported alignment: 478942 (0.65%)
# reads that failed to align: 73024061 (99.35%)
Reported 478942 alignments to 1 output stream(s)


I will try the --verbose switch
sasignor is offline   Reply With Quote
Old 05-14-2012, 09:38 AM   #10
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Quote:
Originally Posted by sasignor View Post
The fact that the command says 'eletrans' and the reference listed here is 'eletran' is not significant, I had named it differently at one point.
Ah, I had not noticed that. Just why do you think that the change is not significant? With your command line Bowtie should be looking for eletrans.*.ebwt files.
westerman is offline   Reply With Quote
Old 05-14-2012, 09:40 AM   #11
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

Its not significant because I am using the correct name for the reference. I indexed the reference on two occasions and the same ebwt files exist for both eletrans and eletran. Variability in whether or not the command I post is eletrans or eletran depends on what day the command was issued.
sasignor is offline   Reply With Quote
Old 05-14-2012, 09:43 AM   #12
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Another suggestion since it does seem like your single-end works, at least to some extent. Try with switch '-v 3'. That should allow for 3 mismatches. It could be possible that your reads simply do not map very well and that BWA was able to pick up the poor matches better.

BTW: If it appears that I am grasping at straws here then, yes, I am. I do not have a good idea of what is happening but rather am going through some troubleshooting steps that I would take.
westerman is offline   Reply With Quote
Old 05-14-2012, 11:13 AM   #13
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

See all those Bs in your quality string? B is the lowest possible quality. It could be that all that sequence is inaccurate, and that's why it won't map.
swbarnes2 is offline   Reply With Quote
Old 05-14-2012, 12:21 PM   #14
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

Its true that the reads shown here have a lot of B's, however when I check the whole file less than 1% of the base pairs have this quality score.
sasignor is offline   Reply With Quote
Old 05-14-2012, 10:41 PM   #15
arvid
Senior Member
 
Location: Berlin

Join Date: Jul 2011
Posts: 156
Default

What does the k-mer plots and sequence overrepresentation plots in FastQC say about your reads? Maybe you do have some problems with unclipped home-brew barcodes or something similar... Do you have in-depth knowledge of the whole sample and read processing pipeline?
If there is no obvious problem with the reads, what happens if you run them through a de novo assembler and BLAST (a few of) the resulting contigs? Maybe you got the wrong files from the sequencing provider, and have reads for a different species?
arvid is offline   Reply With Quote
Old 05-15-2012, 06:09 PM   #16
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

The verbose option crashes terminal every time I use it (after a few hours, but before the mapping is done), so I don't know about that.

The fastqc plots of kmer and overrepresented sequences are rather odd looking, but as this is a transcriptome it is unclear to me what the expectation should be, although from what I can tell by googling around what I have is not unusual for a transcriptome. I did not use barcodes for this dataset.

I have tried doing some additional quality filtering to see if that makes a difference. I will report back on that.

Just mapping the reads that do map and completing the pipeline does produce contigs that blast appropriately.
sasignor is offline   Reply With Quote
Old 05-15-2012, 06:39 PM   #17
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

Additional qc and the addition of the -v 3 option resulted in .034% of the reads mapping in bowtie, so I still have no idea whats wrong.
sasignor is offline   Reply With Quote
Old 05-16-2012, 12:59 AM   #18
arvid
Senior Member
 
Location: Berlin

Join Date: Jul 2011
Posts: 156
Default

Quote:
Originally Posted by sasignor View Post
Just mapping the reads that do map and completing the pipeline does produce contigs that blast appropriately.
Well, that's not the point - obviously these came from the right species. I rather meant to take the reads that do not map (or all, for simplicity), and feed them through e.g. Trinity. Then take some of the contig and blast them to nt to see what species you have sequenced. If you do it for let's say 10 million reads, this goes really fast and should give you an idea whether contamination is a problem...
arvid is offline   Reply With Quote
Old 05-17-2012, 06:35 AM   #19
sasignor
Member
 
Location: Davis, CA

Join Date: Sep 2009
Posts: 28
Default

This is the output of FastQC for one of the files I am using - the other is comparable. Again - it does look unusual for a genome but as far as I can tell not for a transcriptome, and not for transcriptomes I have successfully aligned in the past.
Attached Images
File Type: png duplication_levels.png (29.8 KB, 17 views)
File Type: png kmer_profiles.png (62.9 KB, 18 views)
File Type: png per_base_gc_content.png (32.2 KB, 18 views)
File Type: png per_base_quality.png (20.4 KB, 16 views)
File Type: png per_base_sequence_content.png (63.6 KB, 17 views)
sasignor is offline   Reply With Quote
Old 05-17-2012, 07:22 AM   #20
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Lots of poly-T/A in there.

Try using bowtie2 (not bowtie) in '--local' mode.
westerman 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 10:25 PM.


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