going from RNA seq TopHat output to variant calls

[efoss]

I have RNA seq data aligned to a reference genome using TopHat. I would now like to take these SAM/BAM files as an input and get as an output information for where sequence variants are - chromosome, base pair coordinates, where the variants are (gene names, introns, non-genic regions, etc.), what type of mutations they are (SNPs, substitutions, deletions), what effect they have on amino acid sequences (frame shift, nonsense, missense, silent, etc.) and ideally also whether the variant has been reported as a SNP. DNAnexus will do these things but it's pretty expensive and (I believe) not ideal for RNA seq with splice junctions. Does anyone have suggestions for useful tools?

Thank you.

Eric

[liux]

Any progress on this? would love to see how other people does mRNAseq variants calling with TopHat output.

[Dameon]

I use GATK to call variants from TopHat aligned BAM files. First, you'll need to add @RG information and sort using PICARD tools so as to configure the BAM files for GATK; otherwise, it will fail. Depending on what species you are interrogating, you can then realign around indels and recalibrate the quality scores, or go straight to the Unified Genotyper. Because you would be expecting differentially expressed genes with very low and variable coverage across exons, set the --stand_emit_conf and --stand_call_conf to something really low, like 2, and then use the variant annotater option (-A I think) in the Unified Genotyper to add the ReadPosRankSumTest quality score. Take the VCF file generated by GATK and run it through SNPeff (if human, submit the GATK vcf file to SeattleSNP)and then take that vcf file as raw input to GATK's VariantAnnotator to annotate the raw GATK vcf file. Now filter for what you are interested in. Enjoy.

[efoss]

I use GATK to call variants from TopHat aligned BAM files. First, you'll need to add @RG information and sort using PICARD tools so as to configure the BAM files for GATK; otherwise, it will fail. Depending on what species you are interrogating, you can then realign around indels and recalibrate the quality scores, or go straight to the Unified Genotyper. Because you would be expecting differentially expressed genes with very low and variable coverage across exons, set the --stand_emit_conf and --stand_call_conf to something really low, like 2, and then use the variant annotater option (-A I think) in the Unified Genotyper to add the ReadPosRankSumTest quality score. Take the VCF file generated by GATK and run it through SNPeff (if human, submit the GATK vcf file to SeattleSNP)and then take that vcf file as raw input to GATK's VariantAnnotator to annotate the raw GATK vcf file. Now filter for what you are interested in. Enjoy.

Hi Dameon,

Thanks very much. I've run GATK with DNA but not RNA. Do you see any problem with using GATK with RNA seq? The Broad Institute people are kind of ambiguous about whether it works with RNA seq. Anyway, I'll give it a try. Thanks for the detailed instructions.

Best,

Eric

[Dameon]

Hi Dameon,

Thanks very much. I've run GATK with DNA but not RNA. Do you see any problem with using GATK with RNA seq? The Broad Institute people are kind of ambiguous about whether it works with RNA seq. Anyway, I'll give it a try. Thanks for the detailed instructions.

Best,

Eric

The only problems I forsee of using GATK to call variants from RNA-seq data is the filtering. You want to set the Unified Genotyper as sensitive as possible, don't worry about this as GATK is very aggressive in calling SNPs by default, and then use as many options as possible from VariantAnnotator to whittle down the variants to what you believe to be true SNP calls. It would probably help to use --glm SNP so that you only have to worry about filtering for false positive SNP calls for now. Let me know how everything turns out.

[bharati]

normalization of the aligned data

do we not need to go for any normalization method before calling variations on mRNA Seq data?

[pbluescript]

This is a tricky problem and simply using Tophat with GATK will give you an incredible amount of false positives.
Read the comments on this paper to get an idea of the issues as well as some methods to deal with it:

http://www.sciencemag.org/content/333/6038/53.abstract

Here are several other papers that deal with this issue:

http://genomebiology.com/2012/13/4/r26

http://www.nature.com/nmeth/journal/v9/n6/full/nmeth.1982.html

http://rnajournal.cshlp.org/content/early/2012/07/25/rna.033233.112.abstract

There are more out there too, but the basic idea is that if you want to call variants from RNA Seq data, you have to be very careful.

[sindrle]

One question, how will the difference between single-end and paired-end seq effect SNPs call i mRNAseq?

[crazyhottommy]

you may have a look at this http://allaboutbioinfo.blogspot.com/...53107057687822

[sindrle]

Thats fantastic!
Do you have any more good things to read like this one?

Thanks a lot!

[crazyhottommy]

well, several more here

http://www.rna-seqblog.com/technology/methods/data-analysis/snp-detection/snpir-identification-of-genomic-variants-from-rna-seq-data/

http://www.rna-seqblog.com/technology/methods/data-analysis/unspliced-mapping-tools/blackops-increasing-confidence-in-variant-detection-through-mappability-filtering/

http://www.rna-seqblog.com/technology/publications/analysis-and-design-of-rna-sequencing-experiments-for-identifying-rna-editing-and-other-single-nucleotide-variants/

http://www.rna-seqblog.com/news/commentary/seqing-snps-in-rna-data/

Thats fantastic!
Do you have any more good things to read like this one?

Thanks a lot!

[sindrle]

Its so bad Im so tired of courses today, this was really inspiring! Will read it all tomorrow.

Thank you!!

[sindrle]

I use GATK to call variants from TopHat aligned BAM files. First, you'll need to add @RG information and sort using PICARD tools so as to configure the BAM files for GATK; otherwise, it will fail. Depending on what species you are interrogating, you can then realign around indels and recalibrate the quality scores, or go straight to the Unified Genotyper. Because you would be expecting differentially expressed genes with very low and variable coverage across exons, set the --stand_emit_conf and --stand_call_conf to something really low, like 2, and then use the variant annotater option (-A I think) in the Unified Genotyper to add the ReadPosRankSumTest quality score. Take the VCF file generated by GATK and run it through SNPeff (if human, submit the GATK vcf file to SeattleSNP)and then take that vcf file as raw input to GATK's VariantAnnotator to annotate the raw GATK vcf file. Now filter for what you are interested in. Enjoy.

Why SeattleSNP instead of SNPeff for humans? And which software form the SeattleSNP are you referreing to?