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

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

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.

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.

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

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.

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

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

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

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

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

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

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