SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
RNA-Seq: Comparative Analysis of RNA-Seq Alignment Algorithms and the RNA-Seq Unified Newsbot! Literature Watch 3 07-31-2011 08:08 PM
ChIP-Seq: ZINBA integrates local covariates with DNA-seq data to identify broad and n Newsbot! Literature Watch 0 07-27-2011 04:30 AM
gap alignment and local alignment? mingkunli Illumina/Solexa 3 02-19-2009 12:13 PM

Reply
 
Thread Tools
Old 01-10-2013, 12:19 PM   #1
Eric Fournier
Member
 
Location: Quebec, Canada

Join Date: Jul 2011
Posts: 21
Default Suggested aligner for local alignment of RNA-seq data

Greetings,

I'm trying to analyze the results of Illumina RNA sequencing (~5x150M 100bp PE reads). One problem that we are facing is that for a very large number of our reads, only the first ~50bp are of actual biological material, with the rest consisting of Illumina primers. Would anyone who has faced a similar problem care to suggest an alignment program/parameters to analyze this kind of data? I've tried using bowtie2, but I either get terrible alignment rates using --end-to-end, or I am unable to get any splice junctions using --local.

Thank you very much,
-Eric Fournier
Eric Fournier is offline   Reply With Quote
Old 01-10-2013, 12:21 PM   #2
pbluescript
Senior Member
 
Location: Boston

Join Date: Nov 2009
Posts: 224
Default

You can use something like Trimmomatic to trim off the adapter sequences first.
pbluescript is offline   Reply With Quote
Old 01-14-2013, 12:37 PM   #3
Eric Fournier
Member
 
Location: Quebec, Canada

Join Date: Jul 2011
Posts: 21
Default Trouble with Trimmomatic

I've now spent quite a fair amount of time trying to clean up my sequences with Trimmomatic, but I've been unable to find a set of parameters that gets rid of most of the Illumina adapter sequences.

I've attached an example set of 5 sequences that contain Illumina adapters, as well as the adapter file I've been using (Sequences obtained from UniVec).

When running the R1 sequences through VecScreen, it is quite obvious that the adapters are present:
Code:
Query  53  AGATCGGAAGAGCGGCTCAGCAGGAATGTCGTGACCGATCTCGT  96
           ||||||||||||||| |||||||||||| || ||||||||||||
Sbjct  61  AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGT  18

Query  52  AGATCGGAAGAGCGGTTCAGCA  73
           ||||||||||||||||||||||
Sbjct  61  AGATCGGAAGAGCGGTTCAGCA  40

Query  48  AGATCGGAAGAGCGGTTCAGCAGGAATGACGAGACCGATCTCGTATGCC  96
           |||||||||||||||||||||||||||| ||||||||||||||||||||
Sbjct  61  AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCC  13

Query  52  AGATCGGAAGAGCGGCTCAGCAGGTATGTCGAGACCGATCTCG  94
           ||||||||||||||| |||||||| ||| ||||||||||||||
Sbjct  61  AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCG  19

Query  45  AGATCGGAAGAGCGGCTCAGCAGGTATGCCGAGAGCGATCTCGTATG  91
           ||||||||||||||| |||||||| ||||||||| ||||||||||||
Sbjct  61  AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATG  15
And yet unless I use absurd thresholds (9:3:3), most are left untouched. Am I doing something wrong?
Attached Files
File Type: txt Example_R2.fastq.txt (1.3 KB, 2 views)
File Type: txt Example_R1.fastq.txt (1.2 KB, 9 views)
File Type: txt TrimmomaticAdapters.fa.txt (386 Bytes, 11 views)

Last edited by Eric Fournier; 01-14-2013 at 12:40 PM. Reason: Fixed in-line alignment spacing
Eric Fournier is offline   Reply With Quote
Old 01-14-2013, 04:57 PM   #4
MeganS
Member
 
Location: US

Join Date: Sep 2010
Posts: 14
Default

I have been using trimmomatic for exactly the same situation and I have been happy with the results. Here is the command I am using:

Code:
java -classpath <path_trimmomatic> org.usadellab.trimmomatic.TrimmomaticPE -phred64 file1.fq file2.fq p1.fastq u1.fastq p2.fastq u2.fastq ILLUMINACLIP:./adapter.fasta:2:30:12 SLIDINGWINDOW:4:20 LEADING:10 TRAILING:10
I think the palindrome clipping is really good, but it took a while to figure out how to get the fasta formatted properly. Here are my entries for palindrome clipping, (note they must "start with 'Prefix', and end in '/1' for the forward adapter and '/2' for the reverse adapter"):

Code:
>Prefix/1
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
>Prefix/2
CAAGCAGAAGACGGCATACGAGATNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
Also, check your quality score encoding. Additionally, I modified the ILLUMINACLIP to not require a minimum prefix for palindrome clipping (public final static int MIN_PREFIX = 1; in IlluminaClippingTrimmer.java)
MeganS is offline   Reply With Quote
Old 01-14-2013, 05:00 PM   #5
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

An alternative to trimmomatic would be cutadapt.
chadn737 is offline   Reply With Quote
Old 01-17-2013, 01:54 PM   #6
Eric Fournier
Member
 
Location: Quebec, Canada

Join Date: Jul 2011
Posts: 21
Default

I've moved on from using Trimmomatic to cutadapt, and I've been able to clean up my sequences pretty well. However, I'm running into a new snag: I just can't seem to reliably align spliced reads.

I've built a test subset of spliced reads (attached), which all align to my reference genome (Bos taurus, UMD3.1) when I use Blast or BLAT. I believe that bowtie2 cannot align spliced reads, so I've moved on to TopHat2. However, I've been unable to find a set of parameters which aligns at least a majority of the spliced reads. My best result so far only aligns 6 of the 18. I am using the following command line:

Code:
./tophat-2.0.6.Linux_x86_64/tophat2 -N 6 --read-gap-length 5 --read-edit-dist 10 --splice-mismatches 2 --library-type fr-unstranded --num-threads 6 --b2-very-sensitive ~/bowtie2-indices/Bos_Taurus/Bos_Taurus Splice_R1_cut.fastq Splice_R2_cut.fastq
Could anyone suggest additional/different parameters to increase the amount of successful alignments?

Thanks for any help!
Attached Files
File Type: txt Splice_R1_cut.fastq.txt (1.9 KB, 2 views)
File Type: txt Splice_R2_cut.fastq.txt (2.0 KB, 0 views)
File Type: txt ForBLAST.fa.txt (3.1 KB, 2 views)
Eric Fournier is offline   Reply With Quote
Old 01-17-2013, 02:03 PM   #7
pbluescript
Senior Member
 
Location: Boston

Join Date: Nov 2009
Posts: 224
Default

Since you have shorter reads, you might need to set --segment-length and --min-anchor-length lower.

You could also try STAR. It finds spliced alignments by looking for the largest portion of a read that aligns, softclipping bases that don't align. That will avoid you even having to use something like cutadapt.
pbluescript is offline   Reply With Quote
Old 01-18-2013, 06:20 AM   #8
Eric Fournier
Member
 
Location: Quebec, Canada

Join Date: Jul 2011
Posts: 21
Default

Alright, I'll try those.

I've tried using STAR, but unfortunately I don't have enough RAM to run it, even in sparse mode. I've started looking into using Amazon Web Services or getting some time on a supercalculator in case I can't get TopHat2 to a satisfactory point.
Eric Fournier is offline   Reply With Quote
Old 01-23-2013, 11:30 AM   #9
Eric Fournier
Member
 
Location: Quebec, Canada

Join Date: Jul 2011
Posts: 21
Default

I've finally managed to get reasonable results using Tophat2 by quality trimming my reads. Even though the low-quality read ends were accurate enough for BLAT/BLAST to align them properly, they contained just too many errors for Tophat2, even with parameters allowing for high flexibility.
Eric Fournier is offline   Reply With Quote
Old 01-23-2013, 11:38 AM   #10
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

In the past, when having similar issues, I find that quality trimming followed by adapter trimming gave the optimal results. Just as poor quality base pairs at the end of the read affect alignment, it also can affect adapter trimming. I also tried doing it in reverse order as well as doing simultaneously (cutadapt can both quality trim and adapter trim) but found that doing the quality trimming first resulted in the best overall alignment and the most reads aligned.
chadn737 is offline   Reply With Quote
Reply

Tags
illumina, low quality, rna-seq

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 05:34 AM.


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