SEQanswers

Go Back   SEQanswers > Applications Forums > Genomic Resequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
Calling SNPS following Tophat alignment of RNA-SEQ reads cormicp Bioinformatics 4 05-30-2012 10:17 AM
calling SNPs asankaf General 2 02-04-2009 07:45 PM

Reply
 
Thread Tools
Old 08-29-2013, 01:13 PM   #1
metheuse
Member
 
Location: US

Join Date: Jan 2013
Posts: 77
Default Questions in calling SNPs from RNA-seq reads

I'm trying to call SNPs from RNA-seq data.
The workflow I'm trying to follow is: tophat2 --> remove duplicates --> local realignment and base quality recalibration by GATK --> call SNPs by samtools mpileup or GATK

These are some specific questions:

1. Should I filter the alignment from tophat2 to get only the uniquely mapped reads (reads that could only be mapped to single location on the genome)? If yes, how should I filter them?
I know there is a --multihits (-g) parameter in tophat that can be set to 1 to report only one hit per mapped read. But this is not the strict definition of "unique mapping". I saw people saying that in the alignment sam file, if MAPQ=50, that means this read can only be mapped to a single position. Could I, for example, run tophat2 with g=20 (default), and then select only the aligned reads with MAPQ=50?

2. GATK cannot accept N in cigar. How should I deal with the aligned reads with N in the cigar? GATK has two options: the first is to remove these reads; the second is to allow N in cigar. The former strategy may lose a lot of information, while the latter is at risk (I don't know what the risk is exactly).

3. To use samtools mpileup, how should I set the -C (parameter for adjusting mapQ) parameter? I know it recommends 50 for bwa aligned reads, but what it is for tophat2 alignment?

4. If you have experience in this kind of data analysis, can you see any other potential problems in the pipeline or the parameters of each program? Do you have any recommendation?

By the way, I know there are a lot of caveats of using RNA-seq data to call SNPs. The major purpose of this project is to measure expression, but calling SNPs is to try to get more use of the data (or the money).

Thanks a lot for any suggestions or advice!
metheuse is offline   Reply With Quote
Old 10-15-2013, 10:59 AM   #2
bruce01
Senior Member
 
Location: Dublin

Join Date: Mar 2011
Posts: 156
Default

I've done this a few times (call variants in RNAseq). From my experience it was 'easiest' to follow the GATK Best Practices (http://www.broadinstitute.org/gatk/guide/best-practices). I say 'easiest' because GATK is a bit of a hassle to run (the first time, anyway, second time was ok=).

So I aligned with BWA (Tophat gives potential splice junctions; if you think they may be holding many important variants continue with Tophat, but the majority of SNPs in exon of interest (stops, non-synonymous) will not be in that region I think... and if they are you may be better trying the splice variant/differential splicing analysis...). I used UnifiedGenotyper and HaplotypeCaller to get my VCFs. I tried mpileup before, it was fine but everyone kept saying 'GATK is the gold standard' so I used it in the end: peer (review) pressure=D

So I recommend: BWA, remove duplicates (Picard), realign, recalibrate, call variants, <other>.

Hope that is of some help.
bruce01 is offline   Reply With Quote
Old 10-15-2013, 11:36 AM   #3
metheuse
Member
 
Location: US

Join Date: Jan 2013
Posts: 77
Default

Quote:
Originally Posted by bruce01 View Post
I've done this a few times (call variants in RNAseq). From my experience it was 'easiest' to follow the GATK Best Practices (http://www.broadinstitute.org/gatk/guide/best-practices). I say 'easiest' because GATK is a bit of a hassle to run (the first time, anyway, second time was ok=).

So I aligned with BWA (Tophat gives potential splice junctions; if you think they may be holding many important variants continue with Tophat, but the majority of SNPs in exon of interest (stops, non-synonymous) will not be in that region I think... and if they are you may be better trying the splice variant/differential splicing analysis...). I used UnifiedGenotyper and HaplotypeCaller to get my VCFs. I tried mpileup before, it was fine but everyone kept saying 'GATK is the gold standard' so I used it in the end: peer (review) pressure=D

So I recommend: BWA, remove duplicates (Picard), realign, recalibrate, call variants, <other>.

Hope that is of some help.
Thanks for your information! It's helpful.
metheuse is offline   Reply With Quote
Old 10-15-2013, 11:41 AM   #4
bruce01
Senior Member
 
Location: Dublin

Join Date: Mar 2011
Posts: 156
Default

Hopefully it is helpful, if you have ideas on how to use the variant data specifically in terms of DE genes or otherwise maybe post back?! Also have a look at SNPdat to annotate your VCF, works very well for me.
bruce01 is offline   Reply With Quote
Old 10-15-2013, 11:43 AM   #5
metheuse
Member
 
Location: US

Join Date: Jan 2013
Posts: 77
Default

Quote:
Originally Posted by bruce01 View Post
Hopefully it is helpful, if you have ideas on how to use the variant data specifically in terms of DE genes or otherwise maybe post back?! Also have a look at SNPdat to annotate your VCF, works very well for me.
Thanks, I'll take a look at SNPdat. We usually use SnipEff to annotate VCF.
metheuse is offline   Reply With Quote
Old 01-20-2014, 11:32 AM   #6
FrankiB
Member
 
Location: Sherbrooke

Join Date: Dec 2013
Posts: 23
Default

What are the specific settings (options, parameters) to use when mapping RNA seq data on human genome using BWA software?
FrankiB is offline   Reply With Quote
Old 01-21-2014, 06:34 AM   #7
bruce01
Senior Member
 
Location: Dublin

Join Date: Mar 2011
Posts: 156
Default

To quote Homer: "the two greatest words in the English language: de-fault!"

That is, I use default settings. Also depends on your data, as per: http://bio-bwa.sourceforge.net/bwa.shtml
bruce01 is offline   Reply With Quote
Old 01-22-2014, 07:36 AM   #8
FrankiB
Member
 
Location: Sherbrooke

Join Date: Dec 2013
Posts: 23
Default

OK thanks.

But first I have to create the index. I'm using the human_g1k_v37 version. It has been 2 hours now that I've send the command and I have still this line:

[bwa_index] Pack FASTA...

Is everything okay?
FrankiB is offline   Reply With Quote
Old 01-22-2014, 07:52 AM   #9
bruce01
Senior Member
 
Location: Dublin

Join Date: Mar 2011
Posts: 156
Default

Hi FrankiB, yes I would say everything is ok. As a general rule if there is no error given, and the rest of your computer is working ok, then you just need to wait for things to build. Have a look at 'screen' and 'top' if you don't know them, nice way to make sure your command is actually doing something!
bruce01 is offline   Reply With Quote
Old 01-29-2015, 09:58 AM   #10
amodupe
Junior Member
 
Location: newark, DE

Join Date: Apr 2012
Posts: 1
Default

Quote:
Originally Posted by bruce01 View Post
I've done this a few times (call variants in RNAseq). From my experience it was 'easiest' to follow the GATK Best Practices (http://www.broadinstitute.org/gatk/guide/best-practices). I say 'easiest' because GATK is a bit of a hassle to run (the first time, anyway, second time was ok=).

So I aligned with BWA (Tophat gives potential splice junctions; if you think they may be holding many important variants continue with Tophat, but the majority of SNPs in exon of interest (stops, non-synonymous) will not be in that region I think... and if they are you may be better trying the splice variant/differential splicing analysis...). I used UnifiedGenotyper and HaplotypeCaller to get my VCFs. I tried mpileup before, it was fine but everyone kept saying 'GATK is the gold standard' so I used it in the end: peer (review) pressure=D

So I recommend: BWA, remove duplicates (Picard), realign, recalibrate, call variants, <other>.

Hope that is of some help.
Hi bruce01, what did you use for your VCFs, UnifiedGenotyper or HaplotypeCaller? or did you use both?
amodupe is offline   Reply With Quote
Old 01-29-2015, 11:46 AM   #11
bruce01
Senior Member
 
Location: Dublin

Join Date: Mar 2011
Posts: 156
Default

I use(d) HaplotypeCaller which is now recommended by GATK support as the one to use in all circumstances.
bruce01 is offline   Reply With Quote
Old 03-17-2016, 07:50 AM   #12
exome
Junior Member
 
Location: Los Angeles

Join Date: Mar 2016
Posts: 1
Default

Detecting mutations from RNA-Seq is not a typical approach to detect mutations, mainly due to the intrinsic complexity in the transcriptome (e.g., splicing, high/low gene expression level).
exome is offline   Reply With Quote
Reply

Tags
gatk, rna-seq, snp calling

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 02:06 AM.


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