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 09-13-2011, 01:24 AM   #1
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default Exome sequencing analysis manual

Hi Folks,

As I was writing a short guide of Exome analysis in our Institute, I thought it might be of some use to others especially for newbies, who need some kind of starting point to get to analysis of exome data (pretty much like the RNA-seq manual I once read in an older thread...). Instead of explaining everything in 100 new threads one could then point to that manual...

It is the way we do exome analysis at our Institute, but I would be happy if people help improve the manual, add their knowledge and expand it, like a common knowledge base for exome-level analysis.

I attached the pdf version and a .doc version within a zip folder, as the filesize was too large for uploading the doc file alone.

The most updated version can be found in the SeqWiki (http://seqanswers.com/wiki/How-to/exome_analysis)
(just to make it clear, it is not regularly updated and it's only goal is to get people started on the use of tools often used in exome sequencing)

Any comments highly appreciated!

P.S. I added a (very) short visualization chapter
Attached Files
File Type: pdf Exome-analysis.pdf (334.8 KB, 3002 views)

Last edited by ulz_peter; 04-12-2012 at 10:08 PM. Reason: updated manual
ulz_peter is offline   Reply With Quote
Old 09-22-2011, 07:59 AM   #2
liu_xt005
Member
 
Location: Iowa City, IA

Join Date: Jun 2011
Posts: 24
Default Thanks a lot.

I have learned a lot.
liu_xt005 is offline   Reply With Quote
Old 09-22-2011, 09:26 AM   #3
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,542
Default

You have a weird typo "foolwong" just above the FASTQ example.

Also your introduction about the different FASTQ encodings is out of date now. Illumina now follow the Sanger convention. They also changed the read naming convention, in particular the old /1 and /2 suffixes are gone

See this thread for details:
http://seqanswers.com/forums/showthread.php?t=8895

Also I've never heard FASTQ called a "FastAlignment and Quality file" (glossary on last page).
maubp is offline   Reply With Quote
Old 09-22-2011, 11:50 PM   #4
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

thanks for the hints. As we do not produce Illumina data in ourlab (yet) I haven't heard of those changes, although they seem to have been implemented a while ago...

the typo should mean following, I will rewrite that part and repost it...
ulz_peter is offline   Reply With Quote
Old 10-08-2011, 05:23 PM   #5
pc2009open
Junior Member
 
Location: Taiwan

Join Date: Jun 2009
Posts: 3
Default This is a great document.

Thanks a lot. This is a great document. I wish I had read this document earlier.
pc2009open is offline   Reply With Quote
Old 10-08-2011, 05:50 PM   #6
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

As the GATK local realignment around indels portion of the website does not explicitly state to "FixMateInformation", I am curious if that will affect downstream analysis in anyway?

Great document, by the way.
Heisman is offline   Reply With Quote
Old 10-09-2011, 09:50 AM   #7
Jon_Keats
Senior Member
 
Location: Phoenix, AZ

Join Date: Mar 2010
Posts: 279
Default

Very nice document, thanks for sharing.
Jon_Keats is offline   Reply With Quote
Old 10-09-2011, 01:03 PM   #8
NGSfan
Senior Member
 
Location: Austria

Join Date: Apr 2009
Posts: 181
Default

Great work ulz_peter ! That is exactly what my pipeline is like and I'm glad to see that my choices of tools are also someone else's favorites.

I have seen others use different tools but I think the BWA + GATK + ANNOVAR is the best combination of tools so far...
NGSfan is offline   Reply With Quote
Old 10-09-2011, 04:50 PM   #9
pc2009open
Junior Member
 
Location: Taiwan

Join Date: Jun 2009
Posts: 3
Default

Hi ulz_peter,

I have gone through almost the whole process according to your suggestions. However, at the "3.2. Variant quality score recalibration", I encountered some problems. (I used the TATK-1.0.5506 version.)

I got the error message: "Argument with name '--cluster_file' is missing." However, I did not put "--cluster_file" at all.

I looked at some help documents, and found that this kind of "cluster_file" is supposed to be generated by "GenerateVariantClusters". Have you used GenerateVariantClusters before? Is it necessary?

Thanks again for the wonderful manual.
pc2009open is offline   Reply With Quote
Old 10-09-2011, 07:35 PM   #10
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Quote:
Originally Posted by NGSfan View Post
Great work ulz_peter ! That is exactly what my pipeline is like and I'm glad to see that my choices of tools are also someone else's favorites.

I have seen others use different tools but I think the BWA + GATK + ANNOVAR is the best combination of tools so far...
We use Novoalign + SAMtools... I'm curious if there are any papers out there comparing the methods?
Heisman is offline   Reply With Quote
Old 10-09-2011, 09:52 PM   #11
ulz_peter
Senior Member
 
Location: Graz, Austria

Join Date: Feb 2010
Posts: 219
Default

Hi guys ,
Thanks for all your responses. I must admit that the GATK parts are a little outdated (already). I'm gonna switch to the new version this week and will update the manual accordingly...

@pc2009open: I can't find any hint for the use of a cluster_file argument in variant quality score recalibration... Anyone else had seen that?
ulz_peter is offline   Reply With Quote
Old 10-10-2011, 01:49 PM   #12
NGSfan
Senior Member
 
Location: Austria

Join Date: Apr 2009
Posts: 181
Default

Quote:
Originally Posted by Heisman View Post
We use Novoalign + SAMtools... I'm curious if there are any papers out there comparing the methods?
Papers covering all variations and combinations have been hard to find. I did find one under review (Nature proceedings?) where they claim CASAVA 1.8 comes pretty close to GATK.

I think Novoalign is an excellent aligner, although it requires some tweaking to increase sensitivity on indels that are missed with default settings.

We have done a comparison in our lab with BWA , Stampy, Novoalign, and BFAST. Stampy is the best aligner in our hands (detected more of our SNV and INDEL training set), but Novoalign alignments looked a lot cleaner. I think perhaps with tweaking the gap open penalty for indels, Novoalign might have performed better - just takes some effort to test the parameters more to see if can handle all cases.

GATK is definitely ahead of the game for SNV and indel calling (sensitivity and specificity wise). SAMtools is sufficient - probably you can lean on it if you set the parameters to emphasize specificity instead of sensitivity.
NGSfan is offline   Reply With Quote
Old 10-10-2011, 05:25 PM   #13
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

Thank you so much for posting this pipeline, I've been doing the same for some time. Tomorrow I will post some comments about my results so far.

I think you could sum this pipeline to yours:

https://www.vlsci.org.au/sites/defau...e_20Sept11.pdf

Let's make from this thread a big reference for who is doing exome sequencing ... Please !!!
raonyguimaraes is offline   Reply With Quote
Old 10-10-2011, 05:37 PM   #14
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

One question. How many raw snps you are getting after running Unifier Genotyper for the first time ?

Here I'm getting about 300 000 snps and I think there is something wrong with this numbers ... Shouldn't it be around 20 000 snps?

I'm running my analysis again using a BED file from SeqCap EZ Human Exome Library v2.0 (http://www.nimblegen.com/products/se...tml#annotation) but still ... 300 thousands snps are a lot ...

Last edited by raonyguimaraes; 10-10-2011 at 05:41 PM. Reason: english mistaken ... :)
raonyguimaraes is offline   Reply With Quote
Old 10-10-2011, 06:10 PM   #15
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Quote:
Originally Posted by NGSfan View Post
Papers covering all variations and combinations have been hard to find. I did find one under review (Nature proceedings?) where they claim CASAVA 1.8 comes pretty close to GATK.

I think Novoalign is an excellent aligner, although it requires some tweaking to increase sensitivity on indels that are missed with default settings.

We have done a comparison in our lab with BWA , Stampy, Novoalign, and BFAST. Stampy is the best aligner in our hands (detected more of our SNV and INDEL training set), but Novoalign alignments looked a lot cleaner. I think perhaps with tweaking the gap open penalty for indels, Novoalign might have performed better - just takes some effort to test the parameters more to see if can handle all cases.

GATK is definitely ahead of the game for SNV and indel calling (sensitivity and specificity wise). SAMtools is sufficient - probably you can lean on it if you set the parameters to emphasize specificity instead of sensitivity.
Interesting. Thank you for your post. We do a pretty good job (I think) using the latest SAMtools mpileup command with the -A and -B options and setting a minimum mapping quality per read at 50, but I haven't done anything rigorous to determine what our sensitivity/specificity is. I may go ahead an look at comparing it with GATK.
Heisman is offline   Reply With Quote
Old 10-11-2011, 05:30 AM   #16
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

All right,

Here is my pipeline, some metrics and logs ...

Even after using a bed file from "SeqCap EZ Human Exome Library v2.0" i'm still getting 350191 variants ... Hope someone can help me to get to the promised 20k variants.

I'm planning to switch to GATK 1.2 on the next few weeks ...
Attached Files
File Type: gz gatk_logs.tar.gz (172.5 KB, 97 views)
File Type: txt run.txt (7.5 KB, 274 views)
File Type: pdf GATK Metrics.pdf (438.3 KB, 391 views)
raonyguimaraes is offline   Reply With Quote
Old 10-11-2011, 05:45 AM   #17
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Quote:
Originally Posted by raonyguimaraes View Post
All right,

Here is my pipeline, some metrics and logs ...

Even after using a bed file from "SeqCap EZ Human Exome Library v2.0" i'm still getting 350191 variants ... Hope someone can help me to get to the promised 20k variants.

I'm planning to switch to GATK 1.2 on the next few weeks ...
Are you sure you are filtering out mutations that are not called in the exome? If you call mutations in the whole genome you will get many more than 20k if you are doing a capture protocol due to non-specific hybridization.
Heisman is offline   Reply With Quote
Old 10-11-2011, 05:52 AM   #18
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

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 !
raonyguimaraes is offline   Reply With Quote
Old 10-11-2011, 05:54 AM   #19
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Yeah, I haven't used GATK before so I can't really say but that seemed like the most logical thing to me (we see hundreds of thousands of mutations prior to filtering just the CCDS regions).
Heisman is offline   Reply With Quote
Old 10-11-2011, 05:55 AM   #20
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

You are right, without using this bed file I was getting something like 2 million variants ...
raonyguimaraes 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 02:22 AM.


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