View source for Journal:Bioinformatics workflow for clinical whole genome sequencing at Partners HealthCare Personalized Medicine

{{Infobox journal article
|name         = 
|image        = 
|alt          = <!-- Alternative text for images -->
|caption      = 
|title_full   = Bioinformatics workflow for clinical whole genome sequencing at Partners HealthCare Personalized Medicine
|journal      = ''Journal of Personalized Medicine''
|authors      = Tsai, E.A.; Shakbatyan, R.; Evan, J.; Rossetti, P.; Graham, C.; Sharma, H.; Lin, C.-F., Lebo, M.S.
|affiliations = Partners HealthCare, Brigham and Women’s Hospital, Harvard Medical School
|contact      = Tel.: +1-617-768-8292
|editors      = Weiss, S.T.; Liggett, S.B.
|pub_year     = 2016
|vol_iss      = '''6'''(1)
|pages        = 12
|doi          = [http://doi.org/10.3390/jpm6010012 10.3390/jpm6010012]
|issn         = 2075-4426
|license      = [http://creativecommons.org/licenses/by/4.0/ Creative Commons Attribution 4.0 International]
|website      = [http://www.mdpi.com/2075-4426/6/1/12/htm http://www.mdpi.com/2075-4426/6/1/12/htm]
|download     = [http://www.mdpi.com/2075-4426/6/1/12/pdf http://www.mdpi.com/2075-4426/6/1/12/pdf] (PDF)
}}
{{ombox	
| type      = content	
| style    = width: 500px;	
| text      = This article should not be considered complete until this message box has been removed. This is a work in progress.	
}}
==Abstract==
Effective implementation of precision medicine will be enhanced by a thorough understanding of each patient’s genetic composition to better treat his or her presenting symptoms or mitigate the onset of disease. This ideally includes the sequence [[information]] of a complete genome for each individual. At Partners HealthCare Personalized Medicine, we have developed a clinical process for whole genome sequencing (WGS) with application in both healthy individuals and those with disease. In this manuscript, we will describe our [[bioinformatics]] strategy to efficiently process and deliver [[Genomics|genomic]] data to geneticists for clinical interpretation. We describe the handling of data from FASTQ to the final variant list for clinical review for the final report. We will also discuss our methodology for validating this workflow and the cost implications of running WGS.

'''Keywords''': clinical sequencing, WGS, NGS, next generation sequencing, bioinformatics, validation, precision medicine

==Introduction==
Precision medicine is becoming an increasing focus in medical research.<ref name="CollinsANew15">{{cite journal |title=A new initiative on precision medicine |journal=New England Journal of Medicine |author=Collins, F.S.; Varmus, H. |volume=372 |issue=9 |pages=793–5 |year=2015 |doi=10.1056/NEJMp1500523 |pmid=25635347}}</ref> To achieve the resolution necessary to personalize clinical care, greater attention has been drawn towards higher resolution of the patient genome. [[wikipedia::Next generation sequencing|Next generation sequencing]] (NGS) provided a cost-effective method for targeted sequencing of known disease genes at base pair resolution.<ref name="AlfaresResults15">{{cite journal |title=Results of clinical genetic testing of 2,912 probands with hypertrophic cardiomyopathy: Expanded panels offer limited additional sensitivity |journal=Genetics in Medicine |author=Alfares, A.A.; Kelly, M.A.; McDermott, G. et al. |volume=17 |issue=11 |pages=880–8 |year=2015 |doi=10.1038/gim.2014.205 |pmid=25611685}}</ref> Moreover, the advent of exome sequencing enabled rapid discovery of genes causing Mendelian disorders. While gene panels and exome sequencing have proved fast and cost-effective for delivering genomic results back to the patient, these technologies are limited by our current knowledge of the exome, which changes over time. Additionally, the use of targeted capture may introduce biases to the data, including PCR duplicates, depth of coverage disparities, and failures at difficult to amplify target regions.<ref name="HarismendyEval09">{{cite journal |title=Evaluation of next generation sequencing platforms for population targeted sequencing studies |journal=Genome Biology |author=Harismendy, O.; Ng, P.C.; Strausberg, R.L. et al. |volume=10 |issue=3 |pages=R32 |year=2009 |doi=10.1186/gb-2009-10-3-r32 |pmid=19327155 |pmc=PMC2691003}}</ref>

Practical considerations such as sequencing costs, data processing and maintenance, and data analysis complexities are important considerations when a [[laboratory]] is considering a new NGS program. These issues are amplified in whole genome sequencing (WGS) due to the volume of the data and have long been barriers to entry for [[Clinical laboratory|clinical laboratories]] looking to adopt WGS. Despite the ability of WGS to interrogate the entirety of the genome, clinical interpretation still often focuses on only 3% of the genome (i.e., exome data, pharmacogenomics risk variants, and single nucleotide variants associated with complex disease risk).<ref name="BieseckerTheClinSeq09">{{cite journal |title=The ClinSeq Project: Piloting large-scale genome sequencing for research in genomic medicine |journal=Genome Research |author=Biesecker, L.G.; Mullikin, J.C.; Facio, F.M. et al. |volume=19 |issue=9 |pages=1665-74 |year=2009 |doi=10.1101/gr.092841.109 |pmid=19602640 |pmc=PMC2752125}}</ref><ref name="DeweyClinical14">{{cite journal |title=Clinical interpretation and implications of whole-genome sequencing |journal=JAMA |author=Dewey, F.E.; Grove, M.E.; Pan, C. et al. |volume=311 |issue=10 |pages=1035-45 |year=2014 |doi=10.1001/jama.2014.1717 |pmid=24618965 |pmc=PMC4119063}}</ref><ref name="JiangDetect13">{{cite journal |title=Detection of clinically relevant genetic variants in autism spectrum disorder by whole-genome sequencing |journal=American Journal of Human Genetics |author=Jiang, Y.H.; Yuen, R.K.; Jin, X. et al. |volume=93 |issue=2 |pages=249-63 |year=2013 |doi=10.1016/j.ajhg.2013.06.012 |pmid=23849776 |pmc=PMC3738824}}</ref><ref name="LupskiWhole10">{{cite journal |title=Whole-genome sequencing in a patient with Charcot-Marie-Tooth neuropathy |journal=New England Journal of Medicine |author=Lupski, J.R.; Reid, J.G.; Gonzaga-Jauregui, C. et al. |volume=362 |issue=13 |pages=1181-91 |year=2010 |doi=10.1056/NEJMoa0908094 |pmid=20220177 |pmc=PMC4036802}}</ref> Therefore, WGS services may be overlooked for clinical applications as they trend towards increased costs and longer turnaround times due to a heavier computational load, increased number of variants for analysis, and larger data archives. However, the steadily decreasing cost of sequencing and storage now allow laboratories to consider genome sequencing. WGS, and more specifically PCR-free WGS, also decreases the need to re-sequence each time the coding sequences of targeted regions change, novel genes are discovered, or a new reference genome is released. The balance between cost, turnaround time, accuracy, and completeness has to be addressed when launching a WGS program. Here, we describe the workflow we adopted and the challenges we met supporting the bioinformatics of WGS in a clinical setting.

==Results==
===Bioinformatics validation===
Our pipeline performs robustly as different entry points and different runs of the same sample returned exactly the same variants. Using HapMap sample NA12878 and previous Sanger confirmed regions, we identified thresholds for variants of Quality by Depth (QD) ≥ 4 and Fisher Strand Bias (FS) ≤ 30 as providing the optimal balance of sensitivity and specificity. This sample has been well characterized by our laboratory and also by other groups.<ref name="MacArthurASys12">{{cite journal |title=A systematic survey of loss-of-function variants in human protein-coding genes |journal=Science |author=MacArthur, D.G.; Balasubramanian, S.; Frankish, A. et al. |volume=335 |issue=6070 |pages=823-8 |year=2012 |doi=10.1126/science.1215040 |pmid=22344438 |pmc=PMC3299548}}</ref><ref name="ZookInteg14">{{cite journal |title=Integrating human sequence data sets provides a resource of benchmark SNP and indel genotype calls |journal=Nature Biotechnology |author=Zook, J.M., Chapman, B.; Wang, J. et al. |volume=32 |issue=3 |pages=246–51 |year=2014 |doi=10.1038/nbt.2835 |pmid=24531798}}</ref> Of 425 total confirmed variants, all 425 variants were detected by genome sequencing for a sensitivity of 100% (95% CI: 99.1%–100% for SNVs and 79.6%–100% for indels; Table 1a). In addition, four likely reference sequence errors, positions where only the alternative allele has ever been identified, were correctly genotyped by genome sequencing as homozygous for the alternative allele, three of which had incorrect genotype calls with previous orthogonal assays. In addition to these true positive variants, calls were also made for 21 false positive (FP) variants, including 20 substitutions and one indel (Table 1b). After filtration of variants based on optimal QD and FS thresholds discussed, only one FP remained.

{| 
 | STYLE="vertical-align:top;"|
{| class="wikitable" border="1" cellpadding="5" cellspacing="0" width="70%"
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|'''Table 1.''' Sensitivity and Specificity of WGS. Utilizing HapMap sample NA12878, we previously identified and confirmed via Sanger sequencing 425 variants in 195 genes across ~700 kb of sequence. (a) Analysis of raw variants identified 21 false positive (FP) variants within this region of NA12878. After application of laboratory-defined thresholds (FS ≤ 30 and QD ≥ 4), this limited the number of FP to only one with no increase in false negatives (FN). (b) All 425 Sanger confirmed variants were identified in the WGS data above our quality thresholds. (c) Concordance with 1000 Genomes data for variants across the genome was high, particularly for SNVs. Concordance remained high when matching on predicted genotype.
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|(a) Specificity
 |-
  ! style="padding-left:10px; padding-right:10px;" colspan="2"|Variant type
  ! style="padding-left:10px; padding-right:10px;" colspan="2"|FP (before Thresholds)
  ! style="padding-left:10px; padding-right:10px;" colspan="2"|FP (after Thresholds)
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|SNVs
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|20
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|1
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|Indels
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|1
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|0
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|(b) Sensitivity
 |-
  ! style="padding-left:10px; padding-right:10px;" colspan="2"|Variant type
  ! style="padding-left:10px; padding-right:10px;"|#
  ! style="padding-left:10px; padding-right:10px;"|FN
  ! style="padding-left:10px; padding-right:10px;"|Sensitivity
  ! style="padding-left:10px; padding-right:10px;"|95% Cl
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|SNVs
  | style="background-color:white; padding-left:10px; padding-right:10px;"|410
  | style="background-color:white; padding-left:10px; padding-right:10px;"|0
  | style="background-color:white; padding-left:10px; padding-right:10px;"|100%
  | style="background-color:white; padding-left:10px; padding-right:10px;"|99.1%–100%
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="2"|Indels
  | style="background-color:white; padding-left:10px; padding-right:10px;"|15
  | style="background-color:white; padding-left:10px; padding-right:10px;"|0
  | style="background-color:white; padding-left:10px; padding-right:10px;"|100%
  | style="background-color:white; padding-left:10px; padding-right:10px;"|79.6%–100%
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="6"|(c) Concordance with 1000 Genomes data
 |-
  ! style="padding-left:10px; padding-right:10px;"|Variant type
  ! style="padding-left:10px; padding-right:10px;"|1K Genomes Variants
  ! style="padding-left:10px; padding-right:10px;"|Present in NGS Calls
  ! style="padding-left:10px; padding-right:10px;"|% Present in NGS Calls
  ! style="padding-left:10px; padding-right:10px;"|Present in NGS Calls with Matched Genotypes
  ! style="padding-left:10px; padding-right:10px;"|% Present in NGS Calls with Matched Genotypes
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;"|SNVs
  | style="background-color:white; padding-left:10px; padding-right:10px;"|2,762,933
  | style="background-color:white; padding-left:10px; padding-right:10px;"|2,735,592
  | style="background-color:white; padding-left:10px; padding-right:10px;"|99.01%
  | style="background-color:white; padding-left:10px; padding-right:10px;"|2,730,826
  | style="background-color:white; padding-left:10px; padding-right:10px;"|98.84%
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Indels
  | style="background-color:white; padding-left:10px; padding-right:10px;"|327,474
  | style="background-color:white; padding-left:10px; padding-right:10px;"|299,300
  | style="background-color:white; padding-left:10px; padding-right:10px;"|91.39%
  | style="background-color:white; padding-left:10px; padding-right:10px;"|285,401
  | style="background-color:white; padding-left:10px; padding-right:10px;"|87.15%
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|'''Total'''
  | style="background-color:white; padding-left:10px; padding-right:10px;"|3,090,407
  | style="background-color:white; padding-left:10px; padding-right:10px;"|3,034,892
  | style="background-color:white; padding-left:10px; padding-right:10px;"|98.20%
  | style="background-color:white; padding-left:10px; padding-right:10px;"|3,016,227
  | style="background-color:white; padding-left:10px; padding-right:10px;"|97.60%
 |- 
|}
|}

Our genome-wide comparison of the variants detected through this pipeline compared to those detected at a similar coverage through the 1000 Genomes Project for NA12878 confirmed similar results (Table 1c). Out of 2.7 million SNVs and 285,000 indels detected through our pipeline, we found that 98.8% of these SNVs were also called in the 1000 Genomes dataset. In addition, 97.6% of the total variants called were concordant for genotype. Annotation and variant filtration were also performed on NA12878. Greater than 50 variants were randomly selected for manual inspection to ensure that they were properly annotated and filtered.

===Known regions of poor coverage===
All of the genomes delivered to us from Illumina’s [[Clinical Laboratory Improvement Amendments|CLIA]]-laboratory have at least 30X coverage, with an average coverage of 43X. However, this coverage varies across the coding sequences of the genome and affects both clinically relevant and clinically unknown regions of the genome. A gene-level list of the percentage of callable bases is provided in Supplementary Table S1. In total, of the 1381 genes with at least five asserted Pathogenic or Likely pathogenic variants by clinical laboratories in ClinVar, 94 had <90% coverage across their coding region. The 20 genes with the poorest coverage metrics included many with high prevalence and clinical relevance (Table 2). In our clinical process, for indication gene-list driven analyses, we report back the coverage of analyzed genes, highlighting those with coverage issues. This provides a useful guide to determining potential false negative findings, including those due to regions that are difficult to sequence with NGS technologies.

{| 
 | STYLE="vertical-align:top;"|
{| class="wikitable" border="1" cellpadding="5" cellspacing="0" width="70%"
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="5"|'''Table 2.''' Top 20 Poorly Covered Genes with Clinical Relevance. Clinical relevance is defined as having at least five Pathogenic or Likely pathogenic variants in ClinVar reported in the gene by submitting laboratories or working groups.
 |-
  ! style="padding-left:10px; padding-right:10px;"|Gene
  ! style="padding-left:10px; padding-right:10px;"|# Clinically Significant Variants
  ! style="padding-left:10px; padding-right:10px;"|% Callable
  ! style="padding-left:10px; padding-right:10px;"|Disease
  ! style="padding-left:10px; padding-right:10px;"|Disease Prevalence
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;"|STRC
  | style="background-color:white; padding-left:10px; padding-right:10px;"|8
  | style="background-color:white; padding-left:10px; padding-right:10px;"|20
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Sensorineural hearing loss
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Common
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|ADAMTSL2
  | style="background-color:white; padding-left:10px; padding-right:10px;"|5
  | style="background-color:white; padding-left:10px; padding-right:10px;"|32
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Geleophysic dysplasia
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|CYP21A2
  | style="background-color:white; padding-left:10px; padding-right:10px;"|13
  | style="background-color:white; padding-left:10px; padding-right:10px;"|44
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Congenital adrenal hyperplasia
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Common
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|ARX
  | style="background-color:white; padding-left:10px; padding-right:10px;"|19
  | style="background-color:white; padding-left:10px; padding-right:10px;"|45
  | style="background-color:white; padding-left:10px; padding-right:10px;"|X-linked infantile spasm syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|MECP2
  | style="background-color:white; padding-left:10px; padding-right:10px;"|250
  | style="background-color:white; padding-left:10px; padding-right:10px;"|53
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rett syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Common
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|GJB1
  | style="background-color:white; padding-left:10px; padding-right:10px;"|16
  | style="background-color:white; padding-left:10px; padding-right:10px;"|53
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Charcot-Marie-Tooth disease
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Common
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|ABCD1
  | style="background-color:white; padding-left:10px; padding-right:10px;"|33
  | style="background-color:white; padding-left:10px; padding-right:10px;"|57
  | style="background-color:white; padding-left:10px; padding-right:10px;"|X-linked adrenoleukodystrophy
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Moderate
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|EMD
  | style="background-color:white; padding-left:10px; padding-right:10px;"|11
  | style="background-color:white; padding-left:10px; padding-right:10px;"|57
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Emery-Dreifuss muscular dystrophy
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Moderate
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|G6PD
  | style="background-color:white; padding-left:10px; padding-right:10px;"|16
  | style="background-color:white; padding-left:10px; padding-right:10px;"|58
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Glucose-6-phosphate dehydrogenase deficiency
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Common
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|GATA1
  | style="background-color:white; padding-left:10px; padding-right:10px;"|12
  | style="background-color:white; padding-left:10px; padding-right:10px;"|60
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Dyserythropoietic anemia and thrombocytopenia
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|AVPR2
  | style="background-color:white; padding-left:10px; padding-right:10px;"|15
  | style="background-color:white; padding-left:10px; padding-right:10px;"|62
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Nephrogenic diabetes insipidus
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|EDA
  | style="background-color:white; padding-left:10px; padding-right:10px;"|37
  | style="background-color:white; padding-left:10px; padding-right:10px;"|63
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Hypohidrotic ectodermal dysplasia
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Moderate
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|SLC16A2
  | style="background-color:white; padding-left:10px; padding-right:10px;"|11
  | style="background-color:white; padding-left:10px; padding-right:10px;"|63
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Allan-Herndon-Dudley syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|FLNA
  | style="background-color:white; padding-left:10px; padding-right:10px;"|42
  | style="background-color:white; padding-left:10px; padding-right:10px;"|64
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Otopalatodigital syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|EBP
  | style="background-color:white; padding-left:10px; padding-right:10px;"|24
  | style="background-color:white; padding-left:10px; padding-right:10px;"|64
  | style="background-color:white; padding-left:10px; padding-right:10px;"|X-linked chondrodysplasia punctata
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|RPGR
  | style="background-color:white; padding-left:10px; padding-right:10px;"|17
  | style="background-color:white; padding-left:10px; padding-right:10px;"|64
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Retinitis pigmentosa
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Common
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|TAZ
  | style="background-color:white; padding-left:10px; padding-right:10px;"|17
  | style="background-color:white; padding-left:10px; padding-right:10px;"|64
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Barth syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|IDS
  | style="background-color:white; padding-left:10px; padding-right:10px;"|16
  | style="background-color:white; padding-left:10px; padding-right:10px;"|64
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Hunter syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Moderate
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|FGD1
  | style="background-color:white; padding-left:10px; padding-right:10px;"|8
  | style="background-color:white; padding-left:10px; padding-right:10px;"|64
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Aarskog-Scott syndrome
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Rare
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|GPR143
  | style="background-color:white; padding-left:10px; padding-right:10px;"|6
  | style="background-color:white; padding-left:10px; padding-right:10px;"|65
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Ocular albinism
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Moderate
 |- 
|}
|}

===Cost analysis and scalability===
The majority of our bioinformatics costs of WGS lies in data storage and not in computational processing. As part of our clinical workflow, the FASTQ, BAM, and VCF files are periodically archived on a replicated, secondary storage site. We also store MD5 checksums to ensure data consistency during transfer and storage. These disks are primarily used to store data, and any access to these files must be done after they are copied back onto a faster storage location (e.g., primary storage). Based on the calculations of CPU time used and the lifespan of our cluster as well as the hourly costs of a bioinformatics analyst, we estimated that it costs $77 for the computational time and $128 for the hands-on time necessary to process a genome. However, approximately 300 GB per genome of storage is required to keep files containing unaligned and aligned reads as well as the final variant calls (respective FASTQ, BAM, and VCF files). 

In Table 3, we estimate the annual cost of storage per genome to be around $40–$55, depending on whether the genome is stored on a primary storage device or a secondary (or archival) disk; the bioinformatics process and storage of a genome for one year for our laboratory is approximately $245. Disk storage of sequencing data incurs significant cost in addition to the sequencing and interpretive components, and depending upon laboratory needs and policies, the availability of the data on either active or deep storage may differ. This pipeline is built on a high performance computing (HPC) cluster, so the ability to process many genomes at a given time is heavily dependent on the size of the cluster. With our current infrastructure, we do not run more than 10 genomes simultaneously. The pipeline is scalable with respect to the cluster size, but we have not explored scaling this process to simultaneously handle hundreds or thousands of genomes.

{| 
 | STYLE="vertical-align:top;"|
{| class="wikitable" border="1" cellpadding="5" cellspacing="0" width="70%"
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;" colspan="4"|'''Table 3.''' Cost Analysis for Storage of WGS data. Primary storage assumes unreplicated, active storage with high input/output (I/O) capacity. Secondary storage assumes replicated, deep storage with low I/O capacity. The cost of processing a genome and data retention on the primary and secondary storage for one year is ~$245.
 |-
  ! style="padding-left:10px; padding-right:10px;"|Storage Type
  ! style="padding-left:10px; padding-right:10px;"|Genome/Month ($)
  ! style="padding-left:10px; padding-right:10px;"|Genome/Year ($)
  ! style="padding-left:10px; padding-right:10px;"|Genome/5 Years ($)
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Primary
  | style="background-color:white; padding-left:10px; padding-right:10px;"|4.42
  | style="background-color:white; padding-left:10px; padding-right:10px;"|53.04
  | style="background-color:white; padding-left:10px; padding-right:10px;"|265.20
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|Secondary
  | style="background-color:white; padding-left:10px; padding-right:10px;"|3.48
  | style="background-color:white; padding-left:10px; padding-right:10px;"|41.76
  | style="background-color:white; padding-left:10px; padding-right:10px;"|208.80
 |- 
  | style="background-color:white; padding-left:10px; padding-right:10px;"|'''Total'''
  | style="background-color:white; padding-left:10px; padding-right:10px;"|7.90
  | style="background-color:white; padding-left:10px; padding-right:10px;"|94.80
  | style="background-color:white; padding-left:10px; padding-right:10px;"|474.00
 |- 
|}
|}

==Discussion==
The validated bioinformatics workflow we described above has been used to process hundreds of genomes. We have been updating its functionality and re-validating the components of this process multiple times since initially implementing the pipeline. The modularity of our process allows us to repeat only the annotation, upload, and filtration processes when new or updated annotation resources are available. Similarly, filtrations can be easily re-queried from the persistent variant data in the Oracle database.

We have also compared the results of our NA12878 with the Genome in a Bottle Consortium (GIAB) dataset driven by the National Institute of Standards and Technology (NIST).<ref name="ZookInteg14" /> A similar tool known as the Genetic Testing Reference Materials Coordination Program (GeT-RM) allows users to query across different alignments of the same dataset to see if the variant of interest found by your pipeline also is found by other pipelines. This dataset includes variant data generated by our laboratory. Both of these resources are tremendously valuable in validating clinical genome sequencing.

The limitation of genome sequencing is an important consideration. Although it is often times termed “whole genome,” it has many regions that cannot be accurately determined by current technologies. These regions are inaccessible to current, standard clinical WGS analysis and include regions of high homology for which reads cannot be uniquely mapped, regions where the reference genome contains errors, regions with multiple reference haplotypes, and tandem repeats that extend beyond the sequenced read length. These regions can overlap with clinically relevant areas and can potentially lead to false negative results if not dealt with carefully. In Table 2, we have indicated a subset of these that consistently have the lowest coverage in our WGS assay. Current techniques for managing these regions include repeat-primed PCR for repeat expansion disorders, long-range PCR for genes with homology issues, and allele-based barcoding techniques to anchor reads to unique neighboring regions. Aligning to the GRCh38 reference genome may help for regions with multiple reference alleles, but it is still unclear how those regions will be implemented in a clinical workflow. Eventually the use of longer read technologies, and possibly ''de novo'' assembly, will enable more accurate aligning and detection of variants in these difficult regions.

As genome sequencing becomes more routine, additional enhancements will continuously be made to the process. In our laboratory, these enhancements focus on a user interface for variant filtration to enable live filtering for our clinical staff, identification of copy number and structural variation from the WGS data, and upgrading the pipeline to align to the GRCh38 reference genome. This work requires a team of bioinformaticians with a broad skillset that can interact directly with both the clinical staff and the information technology teams. Finding this dynamic can prove to be one of the most challenging aspects of creating a clinical genome sequencing program.

==Methods and materials==
===Bioinformatics pipeline===
The bioinformatics pipeline processes data from the hard drives delivered by Illumina Clinical Services Laboratory (San Diego, CA, USA) to the variant report files generated for clinical interpretation in a fairly automated fashion. As these hard drives are received, they are accessioned in the laboratory and uploaded to our high performance computing (HPC) environment. The BAM files delivered on the hard drives contain reads aligned by the pipeline developed at Illumina using CASAVA (or ISAAC in recent data deliveries). These files contained both aligned and unaligned reads, so the conversion from BAM back to FASTQ format to decouple mapping information from the read sequence does not incur loss of sequencing information.

Data processing from the FASTQ file to the VCF file was performed primarily using the parameters recommended by each respective software package.<ref name="Picard">{{cite web |url=http://broadinstitute.github.io/picard/ |title=Picard |work=GitHub |author=Li, H.; Broad Institute |accessdate=22 October 2015}}</ref><ref name="LiFast09">{{cite journal |title=Fast and accurate short read alignment with Burrows-Wheeler transform |journal=Bioinformatics |author=Li, H.; Durbin, R. |volume=25 |issue=14 |pages=1754-60 |year=2009 |doi=10.1093/bioinformatics/btp324 |pmid=19451168 |pmc=PMC2705234}}</ref><ref name="LiTheSeq09">{{cite journal |title=The Sequence Alignment/Map format and SAMtools |journal=Bioinformatics |author=Li, H.; Handsaker, B.; Wysoker, A. et al. |volume=25 |issue=16 |pages=2078-9 |year=2009 |doi=10.1093/bioinformatics/btp352 |pmid=19505943 |pmc=PMC2723002}}</ref><ref name="VanDerAuweraFromFast13">{{cite journal |title=From FastQ data to high confidence variant calls: The Genome Analysis Toolkit best practices pipeline |journal=Current Protocols in Bioinformatics |author=Van der Auwera, G.A.; Carneiro, M.O.; Hartl, C. et al. |volume=43 |pages=11.10.1-11.10.33 |year=2013 |doi=10.1002/0471250953.bi1110s43 |pmid=25431634 |pmc=PMC4243306}}</ref> Paired-end alignment of the sequencing reads to the hg19 reference genome was performed using bwa 0.6.1-r104. The aligned reads were sorted and PCR duplicates removed using samtools 0.1.18. Local indel realignment, base quality recalibration, variant calling by UnifiedGenotyper, and variant recalibration were performed using Genome Analysis ToolKit (GATK) 2.2.5 and the recommendations in the Best Practices Workflow by the GATK development team at the Broad Institute (Cambridge, MA, USA). The entire workflow is summarized in Figure 1. For more efficient computing, our validated process includes parallelization by dividing files into smaller pieces for even distribution across the cluster. This method was described in our application of this pipeline to the MedSeq Project, a randomized clinical trial assessing the impact of genome sequencing in clinical practice<ref name="McLaughlinASys14">{{cite journal |title=A systematic approach to the reporting of medically relevant findings from whole genome sequencing |journal=BMC Medical Genetics |author=McLaughlin, H.M.; Ceyhan-Birsoy, O.; Christensen, K.D. et al. |volume=15 |pages=134 |year=2014 |doi=10.1186/s12881-014-0134-1 |pmid=25714468 |pmc=PMC4342199}}</ref> and to the discovery of a putative locus causing a novel recessive syndrome presenting with skeletal malformation and malignant lymphoproliferative disease.<ref name="ChopraInherit15">{{cite journal |title=Inherited CHST11/MIR3922 deletion is associated with a novel recessive syndrome presenting with skeletal malformation and malignant lymphoproliferative disease |journal=Molecular Genetics & Genomic Medicine |author=Chopra, S.S.; Leshchiner, I.; Duzkale, H. et al. |volume=3 |issue=5 |pages=413–23 |year=2015 |doi=10.1002/mgg3.152 |pmid=26436107 |pmc=PMC4585449}}</ref> Similar pipelines have been used successfully in other clinical labs using NGS previously.<ref name="BieseckerTheClinSeq09" /><ref name="JiangDetect13" /><ref name="ChoudharyEval14">{{cite journal |title=Evaluation of an integrated clinical workflow for targeted next-generation sequencing of low-quality tumor DNA using a 51-gene enrichment panel |journal=BMC Medical Genomics |author=Choudhary, A.; Mambo, E.; Sanford, T. et al. |volume=7 |pages=62 |year=2014 |doi=10.1186/s12920-014-0062-0 |pmid=25395014 |pmc=PMC4241214}}</ref><ref name="MillerA26_15">{{cite journal |title=A 26-hour system of highly sensitive whole genome sequencing for emergency management of genetic diseases |journal=Genome Medicine |author=Miller, N.A.; Farrow, E.G.; Gibson, M. et al. |volume=7 |pages=100 |year=2015 |doi=10.1186/s13073-015-0221-8 |pmid=26419432 |pmc=PMC4588251}}</ref><ref name="ReidLaunch14">{{cite journal |title=Launching genomics into the cloud: deployment of Mercury, a next generation sequence analysis pipeline |journal=BMC Bioinformatics |author=Reid, J.G.; Carroll, A.; Veeraraghavan, N. et al. |volume=15 |pages=30 |year=2014 |doi=10.1186/1471-2105-15-30 |pmid=24475911 |pmc=PMC3922167}}</ref>


[[File:Fig1 Tsai JofPersMed2016 6-1.png|700px]]
{{clear}}
{| 
 | STYLE="vertical-align:top;"|
{| border="0" cellpadding="5" cellspacing="0" width="700px"
 |-
  | style="background-color:white; padding-left:10px; padding-right:10px;"| <blockquote>'''Figure 1.''' WGS Alignment and Variant Calling Pipeline. There are multiple entry points to our pipeline where it can be re-triggered due to system failures or outside datasets. Standard processing of genome data from the Illumina Clinical Services Laboratory starts at the top entry point. From there, FASTQ sequences are aligned to the reference hg19 genome using the burrows-wheeler aligner (bwa). Since the alignment is computationally intensive, we divided the sequence files into smaller files. The alignments, known as “raw” BAM files in our pipeline, are processed through a series of steps prior to variant calling. The “final” BAM file is the resulting file after removing PCR artifacts, local indel realignment, and base quality recalibration. This is used as the input file to the variant caller portion of the pipeline. Variant calling happens in two phases, where the variants are identified and then their quality scores are recalibrated in the final VCF file.</blockquote>
 |- 
|}
|}

Successive variant annotation is accomplished through a set of scripts that independently annotate the dataset with a wide collection of information, including: (1) transcript and gene annotations from Alamut (Interactive Biosoftware, Rouen, France); (2) gene annotations from variant effect predictor (VEP); (3) variant annotations from 1000 Genomes Project, ClinVar, and Exome Sequencing Project (ESP); and (4) clinical interpretation from our laboratory maintained in GeneInsight (Cambridge, MA, USA).<ref name="1000GenomesAnInt12">{{cite journal |title=An integrated map of genetic variation from 1,092 human genomes |journal=BMC Bioinformatics |author=1000 Genomes Project Consortium; Abecasis, G.R.; Auton, A. et al. |volume=491 |issue=7422 |pages=56–65 |year=2012 |doi=10.1038/nature11632 |pmid=23128226 |pmc=PMC3498066}}</ref><ref name="AronsonTheGene11">{{cite journal |title=The GeneInsight Suite: a platform to support laboratory and provider use of DNA-based genetic testing |journal=Human Mutation |author=Aronson, S.J.; Clark, E.H.; Babb, L.J.; et al. |volume=32 |issue=5 |pages=532–6 |year=2011 |doi=10.1002/humu.21470 |pmid=21432942 |pmc=PMC3082613}}</ref><ref name="McLarenDeriv10">{{cite journal |title=Deriving the consequences of genomic variants with the Ensembl API and SNP Effect Predictor |journal=Bioinformatics |author=McLaren, W.; Pritchard, B.; Rios, D.; et al. |volume=26 |issue=16 |pages=2069-70 |year=2010 |doi=10.1093/bioinformatics/btq330 |pmid=20562413 |pmc=PMC2916720}}</ref><ref name="NHLBIExome">{{cite web |url=http://evs.gs.washington.edu/EVS/ |title=Exome Variant Server |author=NHLBI Exome Sequencing Project |publisher=University of Washington |accessdate=22 October 2015}}</ref> As a clinical laboratory, we are focused on genomic regions that are known or likely to be associated with disease. Thus, we have limited the annotation of variants called on the genome to only those that lie in our target regions. The target region contains the RefSeq coding sequence for all genes mapped to hg19, pharmacogenomic variants, and association or risk variants taken from the NHGRI GWAS catalog.<ref name="WelterTheNHGRI14">{{cite journal |title=The NHGRI GWAS Catalog, a curated resource of SNP-trait associations |journal=Nucleic Acids Research |author=Welter, D.; MacArthur, J.; Morales, J. et al. |volume=42 |issue=D1 |pages=D1001-D1006 |year=2014 |doi=10.1093/nar/gkt1229 |pmid=24316577 |pmc=PMC3965119}}</ref> These regions are buffered by 50 base pairs (bp). The VCF file is limited to only this target region prior to variant annotation. Additionally, during filtration and coverage analysis, we focus in on the regions most likely to contain clinically reportable disease causing variants. Thus, we created clinical region of interest (ROI) files containing all coding regions of all exons for each gene ±15 bp or ±2 bp for filtration and coverage, respectively.

The sample data and variant annotations are uploaded to an Oracle SQL database. The variant filtration process was designed to query this database and return the variants of interest in an excel spreadsheet for careful clinical review. Variant filtration has been an evolving process, and filtration specifications vary depending on the prevalence, penetrance, and predicted mode of inheritance of the genetic disorder. Table 4 lists some of the common filtrations used to generate the final report. Often, the final variant report is generated using Boolean logic to combine more than one of these filters. The overview of the entire bioinformatics process, from alignment and variant calling to final report generation is outlined in Figure 2.

==References==
{{Reflist|colwidth=30em}}

==Notes==
This presentation is faithful to the original, with only a few minor changes to presentation. In some cases important information was missing from the references, and that information was added. In one case, a project URL had changed from Sourceforge to GitHub, and the updated URL was used.

<!--Place all category tags here-->
[[Category:LIMSwiki journal articles (added in 2016)‎]]
[[Category:LIMSwiki journal articles (all)‎]]
[[Category:LIMSwiki journal articles on bioinformatics‎‎]]