SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
PCR duplicate removal for whole genome sequencing vs. whole exome sequencing cliff Bioinformatics 1 09-27-2011 07:29 AM
Qs in exome sequencing data analysis Maone Genomic Resequencing 4 06-17-2011 07:32 AM
Maone, newbie in exome sequencing and data analysis Maone Introductions 0 06-15-2011 07:11 AM
Hands-on ngs workshop - human exome sequencing and microbial whole genome sequencing vikram Events / Conferences 0 12-08-2010 08:36 PM
GS FLX data analysis software manual drgoettel 454 Pyrosequencing 3 07-14-2009 02:47 AM

Reply
 
Thread Tools
Old 10-12-2011, 11:39 PM   #21
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

I just adapted the manual to fit it in the Wiki How to section:

Any changes, recommendations, complaints, etc. welcome:

http://seqanswers.com/wiki/How-to/exome_analysis
ulz_peter is offline   Reply With Quote
Old 10-12-2011, 11:51 PM   #22
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

Quote:
Originally Posted by raonyguimaraes View Post
On the Unifier Genotyper I'm using the following parameters:

# # #Standard Raw VCF
java -Xmx15g -jar $GATK_DIR/GenomeAnalysisTK.jar -T UnifiedGenotyper \
-l INFO \
-I $OUT_DIR/exome.real.dedup.recal.bam \
-R $REFERENCE \
-B:intervals,BED $EXON_CAPTURE_FILE \
-B:dbsnp,VCF $DBSNP \
-glm BOTH \
-stand_call_conf 50.0 \
-stand_emit_conf 20.0 \
-dcov 300 \
-A AlleleBalance \
-A DepthOfCoverage \
-A FisherStrand \
-o $OUT_DIR/exome.raw.vcf \
-log $LOG_DIR/UnifiedGenotyper.log \
-nt 4

The company where this where done guarantees 30X of coverage ... (http://www.otogenetics.com/human_exome_page.htm)

I know this number should reduce after Variant Recalibrator ... I just want to know how many variants people are getting on this step.

By filtering out mutations you mean using the BED File to call only at the target regions ? If so, yes !
I use stand_emit_conf 10.0 and we usually get ~60k SNPs
ulz_peter is offline   Reply With Quote
Old 10-13-2011, 12:49 AM   #23
mirabilia
Junior Member
 
Location: Italy

Join Date: Nov 2010
Posts: 3
Default

hi folks,
thanks for sharing your expertise...it's a great help for a quite newbie like me.
I'm wondering if this analysis pipeline is suitable also for prokaryotic case or needs some adjustments. In case, could you suggest me some references?

thx!
mirabilia is offline   Reply With Quote
Old 10-13-2011, 12:59 AM   #24
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

I haven't tried it with prokaryiotic samples, but it should work actually (bwa, picard and samtools definitely work with prokaryotic data, not too sure about the GATK though...)

You need to adjust it though, for example index your own reference sequences and analysis depends on what sequence variation you'd expect (this pipeline works for diploid genomes only, though you might use some parts of it for different purposes)

Hope that helps.
ulz_peter is offline   Reply With Quote
Old 10-13-2011, 01:42 AM   #25
hanifk
Member
 
Location: China

Join Date: Oct 2010
Posts: 18
Default

thanks very much
hanifk is offline   Reply With Quote
Old 10-13-2011, 03:53 AM   #26
mirabilia
Junior Member
 
Location: Italy

Join Date: Nov 2010
Posts: 3
Default

Quote:
Originally Posted by ulz_peter View Post
I haven't tried it with prokaryiotic samples, but it should work actually (bwa, picard and samtools definitely work with prokaryotic data, not too sure about the GATK though...)

You need to adjust it though, for example index your own reference sequences and analysis depends on what sequence variation you'd expect (this pipeline works for diploid genomes only, though you might use some parts of it for different purposes)

Hope that helps.
thanks a lot ulz_peter!
Could you please clarify which steps of your pipeline are specifically for diploid genomes in order I can customize for my purposes?
mirabilia is offline   Reply With Quote
Old 10-13-2011, 10:21 AM   #27
Michael.James.Clark
Senior Member
 
Location: Palo Alto

Join Date: Apr 2009
Posts: 213
Default

Very cool! I was planning to put together a little Google Site going through how I analyze exome-seq that's very similar to this. Now I'm not sure I should bother!
__________________
Mendelian Disorder: A blogshare of random useful information for general public consumption. [Blog]
Breakway: A Program to Identify Structural Variations in Genomic Data [Website] [Forum Post]
Projects: U87MG whole genome sequence [Website] [Paper]
Michael.James.Clark is offline   Reply With Quote
Old 10-13-2011, 01:29 PM   #28
Orr Shomroni
Member
 
Location: Netherlands

Join Date: Oct 2011
Posts: 26
Default Thank you guys

Thank you ulz_peter and raonyguimaraes. I'm starting doing NGS quite soon, and being a newbie, this pipeline is quite helpful. Also it seems very similar to what my instructor recommended I should do (BWA, SamTools/Varscan, and Annovar). She also said something about using Sift and Polyphen to predict the effect of the mutation on the gene functionality (continuous score that is benign below a certain threshold, and destructive above it). Anyone familiar with those techniques?
__________________
"Though it may seem that all's been said and done, originality still lives on" - some unoriginal guy who had nothing better to write as his signature
Orr Shomroni is offline   Reply With Quote
Old 10-13-2011, 01:30 PM   #29
Michael.James.Clark
Senior Member
 
Location: Palo Alto

Join Date: Apr 2009
Posts: 213
Default

Quote:
Originally Posted by Orr Shomroni View Post
Thanks you ulz_peter and raonyguimaraes. I'm starting doing NGS quite soon, and being a newbie, this pipeline also seems very similar to what my instructor recommended me to do (BWA, SamTools/Varscan, and Annovar). She also said something about using Sift and Polyphen to predict the effect of the mutation on the gene functionality (continuous score that is benign below a certain threshold, and destructive above it). Anyone knows what I'm talking about?
Annovar can annotate with SIFT and Polyphen now.
__________________
Mendelian Disorder: A blogshare of random useful information for general public consumption. [Blog]
Breakway: A Program to Identify Structural Variations in Genomic Data [Website] [Forum Post]
Projects: U87MG whole genome sequence [Website] [Paper]
Michael.James.Clark is offline   Reply With Quote
Old 10-13-2011, 04:43 PM   #30
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

I think we should all give a try to VAAST as well

http://www.yandell-lab.org/software/vaast.html

A probabilistic disease-gene finder for personal genomes.
Yandell M, Huff C, Hu H, Singleton M, Moore B, Xing J, Jorde LB, Reese MG.
Source

Department of Human Genetics, Eccles Institute of Human Genetics, University of Utah and School of Medicine, Salt Lake City, UT 84112, USA. myandell@genetics.utah.edu

VAAST (the Variant Annotation, Analysis & Search Tool) is a probabilistic search tool for identifying damaged genes and their disease-causing variants in personal genome sequences. VAAST builds on existing amino acid substitution (AAS) and aggregative approaches to variant prioritization, combining elements of both into a single unified likelihood framework that allows users to identify damaged genes and deleterious variants with greater accuracy, and in an easy-to-use fashion. VAAST can score both coding and noncoding variants, evaluating the cumulative impact of both types of variants simultaneously. VAAST can identify rare variants causing rare genetic diseases, and it can also use both rare and common variants to identify genes responsible for common diseases. VAAST thus has a much greater scope of use than any existing methodology. Here we demonstrate its ability to identify damaged genes using small cohorts (n = 3) of unrelated individuals, wherein no two share the same deleterious variants, and for common, multigenic diseases using as few as 150 cases.
raonyguimaraes is offline   Reply With Quote
Old 10-13-2011, 05:09 PM   #31
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

I am not affiliated with VAAST in any way, but I have used it extensively and absolutely love it. I don't want to derail this thread in any way but I can certainly answer questions about it.
Heisman is offline   Reply With Quote
Old 10-13-2011, 11:43 PM   #32
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

Quote:
Originally Posted by mirabilia View Post
thanks a lot ulz_peter!
Could you please clarify which steps of your pipeline are specifically for diploid genomes in order I can customize for my purposes?
I didn't find it yet, but there was a statement on the GATK homepage that the options descirbed there (which are basically pretty mcuh the same as I use) only work for diploid genomes and expected shifts of allele frequency must be adressed. So the question is: what are you planning to do: find rare alleles within some strains, sequence a genetically homogeneous strain...
ulz_peter is offline   Reply With Quote
Old 10-17-2011, 08:37 AM   #33
liu_xt005
Member
 
Location: Iowa City, IA

Join Date: Jun 2011
Posts: 24
Default Problem with SNP-calling

Following ulz_peter's original doc, I have some problem when doing the SNP-calling.

java -Xmx4g -jar /path/GenomeAnalysisTK-1.1-35-ge253f6f/GenomeAnalysisTK.jar \
-glm BOTH \
-R hg18.fa \
-T UnifiedGenotyper \
-I myinput.marked.realigned.fixed.recal.bam \
-D dbsnp132_hg18.txt \
-o myoutput.snps.vcf \
-metrics snps.metrics \
-stand_call_conf 50.0 \
-stand_emit_conf 10.0 \
-dcov 1000 \
-A DepthOfCoverage \
-A AlleleBalance \
-L hg18_exonIntervals.bed

This "-L" option does not work.
I got the hg18_exonIntervals.bed from UCSC as ulz_peter's original doc shows.
I run the SNP-calling without the "-L" line.
Then the variant quality score recalibration step does not work, generating an empty output.tranches file.

Can somebody help me out? Thanks a lot.

Last edited by liu_xt005; 10-17-2011 at 10:19 AM.
liu_xt005 is offline   Reply With Quote
Old 10-17-2011, 09:37 PM   #34
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

What is the error message when you specify the -L argument?

I actually stopped using the Variant quality Score recalibration as it often did not work out for me (I never work on more than 2 exomes at a time).

I out the version withouth the recalibration on the SEQanswers Wiki/How-To section. You may have a look there, as I will update that in the future and stop uploading newer versions of the PDF file...
ulz_peter is offline   Reply With Quote
Old 10-18-2011, 03:07 AM   #35
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Smile

Quote:
Originally Posted by liu_xt005 View Post
Following ulz_peter's original doc, I have some problem when doing the SNP-calling.

java -Xmx4g -jar /path/GenomeAnalysisTK-1.1-35-ge253f6f/GenomeAnalysisTK.jar \
-glm BOTH \
-R hg18.fa \
-T UnifiedGenotyper \
-I myinput.marked.realigned.fixed.recal.bam \
-D dbsnp132_hg18.txt \
-o myoutput.snps.vcf \
-metrics snps.metrics \
-stand_call_conf 50.0 \
-stand_emit_conf 10.0 \
-dcov 1000 \
-A DepthOfCoverage \
-A AlleleBalance \
-L hg18_exonIntervals.bed

This "-L" option does not work.
I got the hg18_exonIntervals.bed from UCSC as ulz_peter's original doc shows.
I run the SNP-calling without the "-L" line.
Then the variant quality score recalibration step does not work, generating an empty output.tranches file.

Can somebody help me out? Thanks a lot.
Did you try with -B:targetIntervals,BED hg18_exonIntervals.bed?

By the way I couldn't figure out how to use this on version 1.2
raonyguimaraes is offline   Reply With Quote
Old 10-18-2011, 03:53 AM   #36
mirabilia
Junior Member
 
Location: Italy

Join Date: Nov 2010
Posts: 3
Default

Quote:
Originally Posted by ulz_peter View Post
I didn't find it yet, but there was a statement on the GATK homepage that the options descirbed there (which are basically pretty mcuh the same as I use) only work for diploid genomes and expected shifts of allele frequency must be adressed. So the question is: what are you planning to do: find rare alleles within some strains, sequence a genetically homogeneous strain...
The first part of the analysis, until the local realignment around indels, is of course suitable for all genomes I have. Quality score recalibration and SNP calling need customization and I'm thinking on... Actually, I'd like to catch rare alleles in sequenced bacterial populations, like repeats, tandem repeats and structural variants, but I suspect it's a really difficult task to discriminate between real rare variations and noise introduced by such data.
Obviously any kind of suggestion, it's really appreciate!
mirabilia is offline   Reply With Quote
Old 10-18-2011, 07:13 AM   #37
liu_xt005
Member
 
Location: Iowa City, IA

Join Date: Jun 2011
Posts: 24
Thumbs up -L option problem solved

Quote:
Originally Posted by ulz_peter View Post
What is the error message when you specify the -L argument?

I actually stopped using the Variant quality Score recalibration as it often did not work out for me (I never work on more than 2 exomes at a time).

I out the version withouth the recalibration on the SEQanswers Wiki/How-To section. You may have a look there, as I will update that in the future and stop uploading newer versions of the PDF file...
Dear ulz_peter and raonyguimaraes,

Thanks VERY MUCH to both of you.
The problem seems to be solved by removing the random and hap intervals from the .bed file.

ulz_peter,
raonyguimaraes posted a similar pipeline by Gayle Philip.
SNPs and Indels are recalibrated/filtered separately and combined after.
I am trying that, and think that it is a good idea to exclude Indels from Gaussian models.
liu_xt005 is offline   Reply With Quote
Old 10-20-2011, 04:16 PM   #38
MQ-BCBB
Member
 
Location: Maryland

Join Date: May 2009
Posts: 25
Default

Great resource!!!
Just a few comments:
In the picard/MarkDuplicates.jar, the option 'CREATE_INDEX=true' should be added.
With respect to adding read group information, instead of using the bwa sampe -r option, picard AddOrReplaceReadGroups.jar is an easier way to go as it tells you which options are required. Thanks for sharing!
MQ-BCBB is offline   Reply With Quote
Old 10-20-2011, 09:16 PM   #39
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

You're absolutely right on the first one.
For the sampe -r option: I'd like to keep it in the bwa part as another step handling the BAM file takes a lot of time than just doing it once... but I could add that as an alternative
ulz_peter is offline   Reply With Quote
Old 10-25-2011, 08:38 AM   #40
emilyjia2000
Member
 
Location: usa

Join Date: May 2011
Posts: 59
Default reference dictionary

When I tried to use GATK to do the local realignment according to ulz_peter's instruction, one error message occurred: Invalid command line: Failed to load reference dictionary. Could anybody let me know where to get this reference dictionary? how to use it in the command line?

Thanks in advance.
emilyjia2000 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 03:22 AM.


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