Question: What Is The Best Pipeline For Human Whole Exome Sequencing?
41
gravatar for Biomed
3.9 years ago by
Biomed3.0k
Bethesda, MD, USA
Biomed3.0k wrote:

We will be getting data from human whole exome sequencing done on Illumina GAIIx. I am planning to put together a pipeline to align (BWA?)(optional as we may get the sequence already aligned)|variant detection (Samtools?)|variant classification-deleteriousness prediction(something like SIFT/PolyPhen)|. I know that there are various options out there but I would like to know based on your experiences what the best collection of these tools is? If you had the chance to start from scratch what kind of a pipeline would you put together? I also know that big sequencing centers have their in house tool sets for these for squeezing the last drop but what is available to the general public is what I am looking for.

Thanks

ADD COMMENTlink modified 6 months ago by vd4mmind40 • written 3.9 years ago by Biomed3.0k
58
gravatar for David Quigley
3.8 years ago by
David Quigley9.1k
San Francisco
David Quigley9.1k wrote:

Here is our basic pipeline. We are novices, but this seems to work for us. Please investigate all of the tools carefully, as (I'll repeat) I am not an expert in this analysis. Note that I am using mostly default values for the analysis; advice I gotten from experienced hands is it's routine to use the intersection of multiple SNP calling methods, and that indel calling is still an art. Note that this is set up for single reads, not paired reads, but the same basic pipeline should apply.

PREPROCESS:

  • Index human genome (Picard), we used HG18 from UCSC.
  • Convert Illumina reads to Fastq format
  • Convert Illumina 1.6 read quality scores to standard Sanger scores

FOR EACH SAMPLE:

  1. Align samples to genome (BWA), generates SAI files.
  2. Convert SAI to SAM (BWA)
  3. Convert SAM to BAM binary format (SAM Tools)
  4. Sort BAM (SAM Tools)
  5. Index BAM (SAM Tools)
  6. Identify target regions for realignment (Genome Analysis Toolkit)
  7. Realign BAM to get better Indel calling (Genome Analysis Toolkit)
  8. Reindex the realigned BAM (SAM Tools)
  9. Call Indels (Genome Analysis Toolkit)
  10. Call SNPs (Genome Analysis Toolkit)
  11. View aligned reads in BAM/BAI (Integrated Genome Viewer)

SAMPLE SCRIPT

/bin/BWA/bwa aln -f /output/FOO.sai -t 3 /seq/REFERENCE/human_18.fasta /seq/FQ/FOO.sanger.fq
/bin/BWA/bwa samse -f /output/FOO.sam /seq/REFERENCE/human_18.fasta /output/FOO.sai /seq/FQ/FOO.sanger.fq
/bin/SAMTOOLS/samtools import /seq/REFERENCE/human_18.fasta /output/FOO.sam /output/FOO.bam
/bin/SAMTOOLS/samtools sort /output/FOO.bam /output/FOO.sorted
/bin/SAMTOOLS/samtools index /output/FOO.sorted.bam /output/FOO.sorted.bam.bai
java -jar /bin/GTK/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /seq/REFERENCE/human_18.fasta -I /output/FOO.sorted.bam  -o /output/FOO.intervals
java -jar /bin/GTK/GenomeAnalysisTK.jar -T IndelRealigner -R /seq/REFERENCE/human_18.fasta -I /output/FOO.sorted.bam -targetIntervals /output/FOO.intervals --output /output/FOO.sorted.realigned.bam
/bin/SAMTOOLS/samtools index /output/FOO.sorted.realigned.bam /output/FOO.sorted.realigned.bam.bai
java -jar /bin/GTK/GenomeAnalysisTK.jar -T IndelGenotyperV2 -R /seq/REFERENCE/human_18.fasta -I /output/FOO.sorted.realigned.bam -O /output/FOO_indel.txt --verbose -o /output/FOO_indel_statistics.txt
java -jar /bin/GTK/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /seq/REFERENCE/human_18.fasta -I /output/FOO.sorted.realigned.bam -varout /output/FOO.geli.calls -vf GELI -stand_call_conf 30.0 -stand_emit_conf 10.0 -pl SOLEXA
ADD COMMENTlink written 3.8 years ago by David Quigley9.1k
1

Very interesting and detailed suggestion. Will try this out once our exome data is ready.

ADD REPLYlink written 3.8 years ago by Khader Shameer14k
1

I haven't yet; our region of interest was quite small, and we only used single reads for that experiment. I think the discussion at http://biostar.stackexchange.com/questions/3925 is a good place to start; the respondents know a lot more than I do about this question.

ADD REPLYlink written 3.0 years ago by David Quigley9.1k
1

Do you do base quality score recalibration? If not, why so? Thanks for sharing!

ADD REPLYlink written 2.6 years ago by Adrian Cortes400

great start. thanks this was what I was looking for when I asked this question. I hope we can build on this and share our experiences as we are not experts either but together I believe we can become experts much faster. I am still hung up on the preprocessing step. how did you Convert Illumina reads to Fastq format and Convert Illumina 1.6 read quality scores to standard Sanger scores? Do you have to convert pipeline 1.6 quality scores to sanger scores?

ADD REPLYlink written 3.8 years ago by Biomed3.0k

To convert Illumina to fastq I think I used MAQ (maq sol2sanger in.sol.fastq out.sanger.fastq). See http://maq.sourceforge.net/fastq.shtml, http://seqanswers.com/forums/showthread.php?t=1801 . It may be that the quality scores are not relevant for BWA, please chime in if you know this. My first alignment used MAQ; I switched to BWA because it can do gapped alignment, and it's faster. seqanswers.com is highly recommended for getting started.

ADD REPLYlink written 3.8 years ago by David Quigley9.1k

I have a question about the workflow provided by you. In the SNP calling script by GATK there is a parameter -vf GELI. what is the meaning of this as i was trying to search this one on the Broad website but there was no description. It will be very helpful if you can describe all the parameters.

Thanks

Saurabh

ADD REPLYlink written 3.7 years ago by Saurabh Baheti0

The version of the UnifiedGenotyper that I used a few weeks ago had a flag -vf which let you output GELI format for the genotype calls. The default was VCF. Apparently they have removed that flag in the latest version and only output VCF, see http://www.broadinstitute.org/gsa/wiki/index.php/Unified_genotyper. To get help with all flags, type java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper --help.

ADD REPLYlink written 3.7 years ago by David Quigley9.1k

Thanks David for the clarification. I have one more question regarding the Indel Calling using GATK. Default parameters for GATK Indel calling are 6X coverage and min of 2 indel-suported read then it will call that variant as indel. But when i compare the samtools indel calls with same parameters, samtools gives double the number of indels than GATK. And on uisng IGV to visualize the exclusive indels by samtools they are not false positive. Just a curious question do you recommend using GATK for indel calls??

ADD REPLYlink written 3.7 years ago by Saurabh Baheti0

Thank you David for precise illustration of whole tools :)

ADD REPLYlink written 3.4 years ago by Thaman2.9k

Hello Saurabh, sorry I missed your comment before. In my limited experience, the GATK and MAQ genotype algorithms give very different results. It's not my specialty, so I can't comment on which is preferable; the only thing we've done is Sanger sequencing to work out a gold standard and adjust backwards from there. Hopefully someone who know more than I can point to published work that clears this up.

ADD REPLYlink written 3.4 years ago by David Quigley9.1k

Hi David, thanks for the pipeline example. I notice that although you indexed the reference genome you don't seem to be using an indexed file on any of the parameters. What am I missing? Cheers

ADD REPLYlink written 3.2 years ago by Chris Beck0

Hi David: Have you used any tools to remove/mark PCR duplicates ?

ADD REPLYlink written 3.0 years ago by Khader Shameer14k

Thanks David. That link is very useful.

ADD REPLYlink written 3.0 years ago by Khader Shameer14k

David: I am able to run the entire steps on my paired end data with minor modifications. I would like to know how did you integrated dbSNP / rsIDs to your .VCF files in the UnifiedGenotyper step.

ADD REPLYlink written 3.0 years ago by Khader Shameer14k

In theory you can get this automatically using the -D parameter to pass in a ROD file, available in the GATK bundle at http://www.broadinstitute.org/gsa/wiki/index.php/GATK_resource_bundle . In practice this bundle was literally not present on the broad FTP server when I went looking for it a few weeks ago, so I just pulled down rsIDs from UCSC using the same genome build as my mapping and appended them to my output file with a python script using a hash of location to rsID.

ADD REPLYlink written 3.0 years ago by David Quigley9.1k

Thanks David, that's really helpful. Also I noted that ROD file is not required for dbSNP132 with latest version of GATK. I can directly use VCF. http://getsatisfaction.com/gsa/topics/dbsnp132_rod

ADD REPLYlink written 3.0 years ago by Khader Shameer14k

Thanks for pointing that out; I think that step should be added to this pipeline. See this NG technical report, particularly figure 3: http://www.nature.com/ng/journal/v43/n5/full/ng.806.html

ADD REPLYlink written 2.6 years ago by David Quigley9.1k
21
gravatar for lh3
3.4 years ago by
lh320k
lh320k wrote:

My this answer better serves as a comment to David Quigley's, which is great, but to emphasize the recent improvements, I still give it as separate answer.

Between step 7 and 8, it is recommended to add an additional step:

samtools calmd -Abr FOO.sorted.bam human_18.fasta > FOO.baq.bam

This will substantially improve SNP specificity.

A few other comments to David Quigley's pipeline are:

  1. It is recommended to use Dindel for indel calling. It is the few indel callers that properly model indels (most not, including IndelGenotyperV2). Broad is reimplementing the Dindel model in GATK. The latest samtools also models indels properly, but it is not evaluated as widely as Dindel.

  2. For the human genome build 36, it is recommended to use the version used by the 1000 genomes project. You will lose genes in the pseudoautosomal regions if you use hg18. Hg19/b37 is fine.

  3. The samtools and GATK SNP callers are converging in algorithms and in results. GATK is arguably better than MAQ, which has been confirmed by several independent studies.

ADD COMMENTlink written 3.4 years ago by lh320k
1

BAQ has been added to the GATK pipeline. See more detailed information about how to use it yourself here:

http://www.broadinstitute.org/gsa/wiki/index.php/Per-base_alignment_qualities_(BAQ)_in_the_GATK

ADD REPLYlink modified 11 months ago by brentp17k • written 3.4 years ago by Michael.James.Clark330

Thanks for the update Heng.

ADD REPLYlink written 3.4 years ago by Biomed3.0k

Hi Hemp,

I am confused as to what set of analysis/commands once needs to do after running calmd? Is the .baq.bam file going to input to what's next? Sorry I am newbie to this analysis and am confused what I am supposed to do next with calmd output.

ADD REPLYlink written 2.3 years ago by Angel130

What does samtools calmd do? I do not understand the word "cap" in man-page definition. Thanks.

ADD REPLYlink written 17 months ago by Noolean1.5k
8
gravatar for Ryan Thompson
3.9 years ago by
Ryan Thompson2.4k
TSRI, La Jolla, CA
Ryan Thompson2.4k wrote:

You seem to have listed all the essential components of an exome sequencing pipeline. Exome sequencing is not yet sufficiently well-established to have a single "best-practice" pipeline available. It's still in the roll-your-own stage. You're going to have to experiment with the options for each component (aligner, SNP-caller, functional annotator, etc) to see which give the best results. You'll probably have to write a lot of glue to make the components fit together.

ADD COMMENTlink written 3.9 years ago by Ryan Thompson2.4k

Thanks for the answer. I do realize that there is no single best practice pipeline but I was hoping people might chip in with their own combinations of "this works best for me" collections. I think I 'll just have to get the answer myself at this point.

ADD REPLYlink written 3.9 years ago by Biomed3.0k
7
gravatar for Darked89
3.9 years ago by
Darked893.5k
Barcelona, Spain
Darked893.5k wrote:

Not exome sequencing specific, but you may start with improving your sequence first:

  • Illumina's pipeline output may differ between versions. Newer version (1.6) should handle better "crowded" lanes

  • you may try to use an alternative basecaller (Ibis http://bioinf.eva.mpg.de/Ibis/)

  • this one use with caution: there are few programs for data correction, i.e. Correction tool for SOAPdenovo.

This should improve numbers of mapped reads. I have to admit I did not test the whole pipeline (I got hijacked by something else for time being), and my goal was not SNP but transcript (any transcript) detection.

As for mapping itself, it depends if you also look for novel splicing / exon skipping etc. If yes, then you need a mapper capable of spliced mapping.

An incomplete list of mappers: http://openwetware.org/wiki/Wikiomics:RNA-Seq

ADD COMMENTlink written 3.9 years ago by Darked893.5k
6
gravatar for Erik Garrison
3.4 years ago by
Erik Garrison1.1k
Somerville, MA
Erik Garrison1.1k wrote:

For simultaneous SNP, MNP, and INDEL calling you can also use FreeBayes, which has the added benefit of optionally using other samples from the same population to increase callset quality.

It generates VCF from BAM files. It benefits greatly from realignment preprocessing steps which David Quigley describes.

ADD COMMENTlink written 3.4 years ago by Erik Garrison1.1k
3
gravatar for Laura
3.4 years ago by
Laura1.2k
Cambridge UK
Laura1.2k wrote:

For functional annotation of your snps/indels once you have them you might want to look at the ensembl variation effect predictor which can be used both on the website for small number of snps or using the api for large numbers

http://www.ensembl.org/tools.html

ADD COMMENTlink written 3.4 years ago by Laura1.2k
3
gravatar for jvijai
3.4 years ago by
jvijai960
jvijai960 wrote:

Thank you for this wonderful thread which is full of useful information from my perspective.

If the platform was SOLiD, what variation of this excellent post and (kindly provided) code from David Quigley would be required?

1) Does use of human_g1k_v37 provide advantage over HG18? (sorry if this has been addressed elsewhere. could not find the specific answer).

2) When using BWA with Colorspace data. There is this thread in Seqanswers mentioning possible problems (longer read lengths while indexing). http://bit.ly/hB4IO8 and there is mention of colorspace support on BWA manual http://bit.ly/dgl1XD

3) Any issues noted with GATK. I have not used GATK, but have overheard conversations between people who have tried about having trouble with atleast some aspects of it.

Thanks

JVJ

ADD COMMENTlink written 3.4 years ago by jvijai960
3
gravatar for Roman Valls Guimerà
2.8 years ago by
Stockholm
Roman Valls Guimerà350 wrote:

We've been using Brad Chapman's one for a while now in production ourselves here in Sweden:

https://github.com/chapmanb/bcbb/blob/master/nextgen/README.md

ADD COMMENTlink written 2.8 years ago by Roman Valls Guimerà350
1
gravatar for Kevin
3.9 years ago by
Kevin440
Kevin440 wrote:

You might wanna check out usegalaxy.org for NGS pipeline. There are others as well. offhand I cannot remember them.

ADD COMMENTlink written 3.9 years ago by Kevin440
6

I know galaxy but I can not see how I can use it to 1. Align 2. Variant call 3. Variant deleteriousness predict. It is nice for manipulating common file formats, convert between them and matching your data against common public data sets but this was not my question.

ADD REPLYlink written 3.9 years ago by Biomed3.0k
5

Hi @Kevin. Thanks for contributing to BioStar. Your reply is rather shortish and vague. Please consider putting more time to write answers that are clear, detailed, and useful. You will thus be more helpful to the members of this community. Cheers

ADD REPLYlink written 3.9 years ago by Eric Normandeau7.1k
1
gravatar for Gard Thomassen
3.2 years ago by
Gard Thomassen10 wrote:

Interesting thread.

Does GATK work well with exome sequencing when it comes to the realigning part? When using it on RNA-seq I run into truble as it was originally designed for complete genome sequencing and it kind of freaks out when it finds that you only have transcript data and no data for the introns... This messes up the sliding window algorithm used in the realignment part. Have you fixed the code to take care of this? A fix was suggested on the Broad-help-site but I have not gotten around to pursue it. Or has the GATK been upgradeded to cope with exome and RNA-seq NGS data?

Best regards Gard

ADD COMMENTlink written 3.2 years ago by Gard Thomassen10

Truth be told I don't know, Gard. When I posted our pipeline we had just done a genome sequencing project; we were "about to do" an exome project "very soon" but I haven't got the data yet. At least I did say I was a novice at this kind of analysis. If you find a good answer please post, as my answer is deficient without it.

ADD REPLYlink written 3.2 years ago by David Quigley9.1k
1
gravatar for User 9126
2.8 years ago by
User 912630
User 912630 wrote:

This is really wonderful thread.

I am little bit confused in using this pipeline if i have say 10 samples.

Should I first align all files to the genome independently,combine as single bam file and proceed further.? or should I need to process every file independently through all these steps?

Thanks Santhosh

ADD COMMENTlink written 2.8 years ago by User 912630

This is not a discussion forum. Questions like yours usually go as a separate item on this site. Link to this question when asking your question.

ADD REPLYlink written 2.8 years ago by Hranjeev1.2k
0
gravatar for Michael Gombert
2.6 years ago by
Michael Gombert0 wrote:

Since i was very unsatisfied with confirming variants produced by CASAVA, I tried to setup this pipeline described abouve by Quigley and L for our HiSeq/SeqCapEZ exome data.

After a long weekend and some hassle with incompatible Java versions and a reference concatenated (hGRC37, but copied together in a non-karyotypic way) in the wrong way, I now make it to the 'samtools calmd...' part.

But all I get back is the error message 'Floating point exception', which I have no idea how to deal with. I feel short of going beserk. Anyone seen this before that can tell me how to deal with it?

Best regards, Michael Gombert

ADD COMMENTlink written 2.6 years ago by Michael Gombert0
1

Michael, this is not a discussion thread, although it may appear that way, if you have a question to ask please post it as a new question so it can be answered and archived to help people in the future

ADD REPLYlink written 2.6 years ago by Daniel Swan10k
0
gravatar for vd4mmind
6 months ago by
vd4mmind40
vd4mmind40 wrote:

Thank you for the valuable thread. I have some more query for which I need some suggestions, I am new to GATK and want to use it for my exome sequencing data analysis. I have been a bit lost reading all the blogs , comments and the technical forums. So here is something I want to say and please correct and guide me through the procedure. I have downloaded the hg19 files from the UCSC browser and created the reference genome but do I need to again use the one which is there in GATK repository and then align my samples for downstream analysis? Also I want to run the GATK in my institute cluster. So if am not wrong I should create the directory of the latest GATK version and transfer all the necessary files via Filezilla in the cluster directory with the same name. Now this I have already done. So next thing is to download the bundle from the repository where I see 2 versions , so which one should I download? 2.5 or 2.3? Also once I download the bundle do I have to download anything else? So here it is which I should be downloading right in my cluster. The jar file and the resource folder with the .java files and then in the main directory of the GATK version folder in my cluster I should download the bundle version (2.5 or 2.3) and then unzip all the files that are there in the bundle directory. Right? Please let me know. Then I should be ready to use the GATK for the different downstream processes listed below:

Identify target regions for realignment (Genome Analysis Toolkit) ->Realign BAM to get better Indel calling (Genome Analysis Toolkit) ->Reindex the realigned BAM (SAM Tools) ->Call Indels (Genome Analysis Toolkit) ->Call SNPs (Genome Analysis Toolkit)->View aligned reads in BAM/BAI (Integrated Genome Viewer)

Please let me know if this looks correct or not. The VCF files from the 1kG and the DBSNP are already there in compressed form in the bundle repository of the GATK website which I am currently downloading and I can use them directly after unzipping them.

ADD COMMENTlink written 6 months ago by vd4mmind40
1

You're likely to get the same comment as the two above questions did, namely that this is not a discussion thread and that you should ask your question separately.

ADD REPLYlink written 6 months ago by jobinv890

I have created a new thread, suggestions are welcome Need help with the GATK usage

ADD REPLYlink written 6 months ago by vd4mmind40
Please log in to add an answer.

Help
Access
  • RSS
  • Stats
  • API

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.0.0
Traffic: 575 users visited in the last hour