slide1 n.
Skip this Video
Loading SlideShow in 5 Seconds..
Lisa Strug The Hospital for Sick Children University of Toronto PowerPoint Presentation
Download Presentation
Lisa Strug The Hospital for Sick Children University of Toronto

Loading in 2 Seconds...

play fullscreen
1 / 27

Lisa Strug The Hospital for Sick Children University of Toronto - PowerPoint PPT Presentation

  • Uploaded on

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

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

Lisa Strug The Hospital for Sick Children University of Toronto

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
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

  • 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
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
  • 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
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
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
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
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
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 l ogistic regression
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 data by skotte et al 2012
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:


calculation of e g ij d ij
Calculation of E(Gij|Dij)
  • The calculation of the expected genotype given the sequencing data is given by


  • 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 g ij d ij depends on read depth
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
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
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
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
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
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
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
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
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 groups1
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
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
Rolandic Epilepsy NGS Association Analysis

* Build 37; approximately 450 variants analyzed

  • 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 and Publication
  • Code is available at:
  • 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
  • AndriyDerkach, graduate student in statistics at University of Toronto
  • Deb Pal, Richard Houlston and Ian Tomlinson