270 likes | 281 Views
Case-Control Genetic Association Studies with Next-Generation Sequence Data and External Controls: The Robust Variance Score. Lisa Strug The Hospital for Sick Children University of Toronto
E N D
Case-Control Genetic Association Studies with Next-Generation Sequence Data and External Controls: The Robust Variance Score Lisa Strug The Hospital for Sick Children University of Toronto Emerging Statistical Challenges and Methods For Analysis of Massive Genomic Data in Complex Human Disease Studies June, 2014 Banff
Motivation • Association studies using next generation sequencing (NGS) are expensive • Public NGS data exists (e.g. 1000 genomes, Complete Genomics) • Can we use our NGS cases with (publicly available) ‘out of study’ sequenced control groups in genetic association studies? • Supplement our control NGS data with public NGS data? • Or as the only control group, much like the Wellcome Trust Case Control Consortium (2007) did for GWAS with SNPs
Benefits of Using External NGS Control Groups • Using (publicly) available control groups would reduce costs • Allows one to focus on sequencing cases • Could increase statistical power • Would provide a means to prioritize variants for follow-up genotyping or functional studies • Could be used to increase sample size of sequenced controls
Challenges • Several factors could bias association studies when cases and • controls are sequenced as part of different experimental designs: • Study design considerations e.g. case-control differences in ethnicity or other confounders • Factors related to sequencing technology and parameters • Enrichment and base calling procedures (different platforms) • Alignment (algorithm and reference genome) • Read depth (LRD results in biased allele frequency) • SNP detection and genotype calling algorithms (e.g GATK, SAMtools)
Effects of read depth and platforms on genotype calling Sequencing on Ion Torrent PGM (lower coverage for AT-rich genome and higher error rate*) Sequencing on IlluminaHiSeq 2000 (lower error rate *) Raw reads with base quality (FastQ file) Map and local re-map to reference genome Map and local re-map to reference genome Allignedreads with map qualities (BAM file) SNP detection and genotype calling VCF files with genotype calls and genotype probabilities etc • Genotype probability P(genotype|Data) for each sample. • SNPs and Indels by joining info across samples * Quail et al, BMC Genomics 2012, 13:341
Bias When Using Genotype Calls • Systematic differences in read depth can effect Type I error (Kim et al. 2011) due to differential missclassification/screening (rare homozygotes miscalled/screened more in the LRD sample) • Simulation: Genotype calls for 1000 variants • 50 cases at 100x and 150 controls at 4x R, selection threshold – filters out called genotypes that have low confidence, assigns them missing. E.g.
NGS Case-Control Association Designs • Sequence cases for variant discovery and follow-up with genotyping (Liu and Leal, 2012; Longmate et al., 2010; Sanna et al. 2011) • Cost effective and can control Type I error inflation • Cannot detect protective variants present only in discovery sample and may beoverly conservative • Adjusting for read depth or weighting variant calls by quality score (Daye et al. 2012; Garner 2011) • Parameters are not estimable when cases and controls distinguishable by read depth • Randomly down-sampling BAM files (available in GATK toolkit) • Not as powerful as an approach using all data
Our Proposal • Re-purpose and extend an approach by Skotte et al (2012) that incorporates genotype uncertainty into the association score test by using genotype likelihoods • In comparison to using called genotypes • Can improve power • Avoid spurious findings • Limit the need for some subjective quality thresholds e.g. R
Workflow proposed for NGS association using external NGS controls &Different alignment algorithms are implicitly accounted for by the RVS because the unit of analysis is genotype probability rather than the genotype calls in the association analysis.
Recall: Score Statistic for Logistic Regression • - phenotype value for sample i (1-case, 0-control). • - Genotype information for sample i at jth variant. • To find out whether the genotype is associated with the phenotype, we use • The score statistic used to test the null hypothesis • And the corresponding score test statistic
Association for Sequencing Databy Skotte et al. (2012) • Notation: • Dij - sequence data (reads and errors, BAM file) for jth variant of ithsubject. • Gij– unobserved genotype (coded as 0, 1, 2). • Joint probability of the observed data (Yi, Dij) for i=1,…,n • With • The score used to test the null hypothesis H0: (1)
Calculation of E(Gij|Dij) • The calculation of the expected genotype given the sequencing data is given by Where • The genotype likelihood is calculated from all the reads of the sample and is provided in the output of standard genotype calling packages, eg .VCF file, or calculated from the aligned reads using the simple Bayesian genotyper (McKenna et al, 2010). • Genotype frequencies are calculated from the full sample • by the EM algorithm (McKenna et al. 2010; Skotte et al. 2012) • Also need to calculate the var(Sj) and therefore var[E(Gij|Dij)]
Variance of E(Gij|Dij) Depends on Read Depth • The law of the total variance: • At high read depth, because: • goes to 1 for true genotypein • And • When read depth is not • high enough Therefore, is readdepth dependent
Robust Variance Estimation for the true variances • If read depth is the same for cases and controls, the logistic regression estimate of the Var(Sj) could be used, otherwise, this estimate is biased • Bias is a function of Ncontrols and Ncaseswhere, for Ncontrols >> Ncases , Var(Sj) is underestimated • Thus variance estimation must distinguish between the two groups • We estimate by with estimated genotype frequencies • And is estimated by its sample variance in controls
Robust Variance Score (RVS) Test • Is the score test with the proposed variance estimation • Score test for single SNP analysis: • For joint analysis of J rare variants, combine vector S=(S1,…,SJ)into single test statistic by common approaches such as • CAST (Morgenthaler and Thilly 2007), Weighted Sum (Madsen and Browning 2009), SKAT (Wu et al. 2011), SKAT-O (Lee et al. 2012). • Estimate covariance matrices for vector S separately in cases and controls and combine as previously done. • Evaluate test statistic by asymptotic distribution or by bootstrap resampling (Hall and Hart, 1990) .
Evaluating the RVS • I. Using Simulation to assessPower and Type I error • II. RVS applied to 1000 Genomes data: high read depth (HRD) exomes versus low read depth (LRD) whole genomes • III. RVS applied to Rolandic Epilepsy NGS HRD cases versus 1000 Genomes LRD controls, in a previously identified region of association
I. Evaluating RVS via Simulation • Simulation Setting: • single variant analysis under the null, MAF =1%,10%, 20%, 30%, 40% • Rare variant analysis: collapse 5 rare variants, MAF ranging from 0.1% to 5% • 500 cases, 100x average read depth with error N(0.01,0.025) • 1500 controls, 4x average read depth with error N(0.01,0.025) • We present RVS with CAST; similar results for other tests • Under the alternative, all 5 variants have OR=1.5
Single Variant Analysis Type I Error QQ plot of p-values for an analysis of 1000 variants simulated under the null model. 500 high read depth cases and 1500 low read depth controls. MAF equal to A) 0.01, B) 0.1, C) 0.2 and D) 0.4.
RVS Type I Error and Power: Compared to True Genotypes Quantile-quantile plot of p-values Table: Empirical power
II. Evaluating RVS Using 1000 Genomes Data • The 1000 Genomes project samples (CEU + GBR), phase 3 release [20130502] • One sample consists of exome data on 56 individuals (~50x) • Second sample consists of 113 individuals sequenced at low read depth (~6.5x) • Aligned chromosome 11 reads downloaded • Use GATK (DePristo et al. 2011) on the combined sample of aligned reads, generate multi-sample VCF file • Apply common filters, then extract genotype calls and genotype likelihoods from the VCF file • Compared RVS to the score test using genotype calls
RVS Type I Error in 1000 Genomes Control Groups Quantile-quantile plot of p-values • Single SNP analysis, MAF>5%
RVS Type I Error in 1000 Genomes Control Groups Quantile-quantile plots of p-values • Analysis with CAST, MAF<5%, 5 variants
III. RVS in NGS cases with Rolandic Epilepsy (RE) and public controls (1000 genome) • 27 individuals of European decent with RE, sequenced (~197x) in a previously linked and associated 600kb region of chr 11 • Compare them to 113 whole genome controls from 1000 genomes (MAF>0.05) • Generate multi-sample VCF file on the combined set with GATK • Apply common filters, then extract genotype calls and genotype likelihoods • Compare RVS results to an analysis with genotype calls using a sample of 200 Europeans sequenced by Complete Genomics (~35x)
Rolandic Epilepsy NGS Association Analysis * Build 37; approximately 450 variants analyzed
Conclusions • Score test based on genotype calls/likelihoods has inflated type I error when case and control groups have different read depths • Read depth bias can be avoided if both cases and controls are HRD. • We cannot guarantee HRD in both groups at every locus. • RVS allows one to incorporate external control groups into NGS association studies; assuming matching on other factors • RVS can be used for single or joint rare variant analysis. • RVS can be extended to accommodate in-study controls augmented with an out-of-study control group. • RVS will be extended to accommodate covariate adjustment.
Code and Publication • Code is available at: • https://github.com/strug-lab • Strug.research.sickkids.ca • DerkachA, Chiang T, Gong J, Addis L, Dobbins S, Tomlinson I, Houlston R, Pal DK, Strug LJ. 2014. Association analysis using next generation sequence data from publicly available control groups: The robust variance score statistics. Bioinformatics, Epub. PMID: 24733292
Acknowledgments • AndriyDerkach, graduate student in statistics at University of Toronto • Deb Pal, Richard Houlston and Ian Tomlinson