SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Initial QC and grooming for Illumina HiSeq2000 paired end RNAseq on Galaxy lindseykelly RNA Sequencing 5 07-30-2014 01:09 PM
Tool for obtaining counts for paired-end strand specific RNAseq data? JenBarb RNA Sequencing 2 03-26-2013 10:59 AM
Problems with TopHat and Paired End RNASeq quique_vzquez Bioinformatics 2 12-14-2012 08:30 AM
Variant calling for high-coverage Illumina data dgmacarthur Bioinformatics 19 04-08-2011 06:34 AM
How do variant callers deal with overlapping paired end reads? krobison Bioinformatics 1 04-30-2010 11:58 AM

Reply
 
Thread Tools
Old 03-29-2013, 12:33 PM   #1
ron128
Member
 
Location: Mumbai

Join Date: Sep 2011
Posts: 38
Default Variant Calling from paired end RNAseq data

Hello everybody, since everybody is getting on the variant calling from RNAseq bandwagon, my P.I wants to get in on the party as well. Too bad for overworked ppl like me So I have this huge dataset of illumina paired end rnaseq which I am trying to run through GATK. I am taking the accepted_hits.bam file from tophat/cufflinks as the input file for GATK. My reference for the tophat/cufflinks pipeline was UCSC annotated genes (knownGenes.gtf) and UCSC hg19 reference.

My pipeline for gatk is:
convert bam>sam, sort sam> insert read groups> fixmates using picard> sam to bam> remove duplicates > reindex bam > realign indels

When i Run the indel realignment for GATK, I get the following error: contig chr 1 missing from reference. So i go ahead and look it up, and use the reorder sam option in picards tools. I modify my pipeline as follows:

convert bam>sam, sort sam> insert read groups> fixmates using picard> createdictionary.jar for hg 19 reference using picard > reorder sam >
sam to bam> remove duplicates > reindex bam > realign indels

My error is not solved even after reordering the sam file and i still get the same error of "chr1 contig not found in your reference"

Is it something to do with the references I have been using? I have used hg19 reference both for my tophat/cufflinks as well as the GATK pipeline.

Thanks a ton in advance!
ron128 is offline   Reply With Quote
Old 03-29-2013, 01:18 PM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

How are you making your .sam? I'm pretty sure that bwa sampe will add read group in for with the -r option.

And I think Picard would add them to a .bam. I'm pretty sure you do NOT have to expand your .bam to a .sam in order to do that.

And you double-cheked to make sure that the name of Chr 1 is exactly the same between your .bam and your reference genome?
swbarnes2 is offline   Reply With Quote
Old 03-30-2013, 12:33 AM   #3
ron128
Member
 
Location: Mumbai

Join Date: Sep 2011
Posts: 38
Default

Hey thanks for a quickie reply Like i said, I am not redoing any alignments here. I have already run the tophat cufflinks pipeline on my data to assemble it into known transcripts using UCSC genes and hg 19 as a reference. SO this way I already have access to an "accepted_hits.bam" as an output from the tophat runs. I am using this accepted_hits.bam as my alignment file for GATK and converting this to a sam file using Picard tols. I am skimming down on my analysis time by not redoing the alignments.

to your query about whether I checked the reference and my sam file for the chr1, yea i did. And this is what is troubling me. I am using the same reference for both my rnaseq and variant calling. So logically there shouldnt be any discrepancies. I was wondering if anybody has faced the same issues when using gatk for variant calling in RNAseq? I am pretty new to this so I might me messing up somewhere in the pipeline..
ron128 is offline   Reply With Quote
Reply

Tags
gatk, rnaseq, rnaseq variant calling, variant 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 03:57 PM.


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