Transcript Strug_Banff2014

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:
1. Study design considerations e.g. case-control differences in
ethnicity or other confounders
1. 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
Illumina HiSeq 2000
(lower error rate *)
Map and local re-map
to reference genome
Sequencing on Ion Torrent
PGM (lower coverage for
AT-rich genome and higher
error rate*)
Map and local re-map
to reference genome
Raw reads with base
quality (FastQ file)
Alligned reads with
map qualities (BAM
file)
SNP detection and genotype calling
• Genotype probability P(genotype|Data) for each
sample.
• SNPs and Indels by joining info across samples
* Quail et al, BMC Genomics 2012, 13:341
VCF files with genotype
calls and genotype
probabilities etc
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.
æ P(G | D ) ö
(1)
ij
÷÷ > R
log10 çç
è P(G(2) | Dij ) ø
NGS Case-Control Association Designs
1. 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 be overly conservative
2. 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
3. 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
 Yi - phenotype value for sample i (1-case, 0-control).
 Gij- Genotype information for sample i at jth variant.
• To find out whether the genotype is associated with the
phenotype, we use
logitP(Yi =1| Gij = gi ) = b0 + b1gi
• The score statistic used to test the null hypothesis H 0 : b1 = 0
_ö
æ
S = åçYi -Y ÷gi
è
ø
i=1
n
• And the corresponding score test statistic
S 2j
Tj =
var(S j )
Association for Sequencing Data
by Skotte et al. (2012)
• Notation:
 Dij - sequence data (reads and errors, BAM file) for jth variant of ith subject.
 Gij – unobserved genotype (coded as 0, 1, 2).
• Joint probability of the observed data (Yi, Dij) for i=1,…,n
n
n
2
 P(Yi , Dij )    P(Yi | Gij  g ) P( Dij , Gij  g )
i 1
• With
i 1 g 0
logitP(Yi = 1| Gij = g) = b0 + b1g
The score used to test the null hypothesis H0: 1  0
n
S 2j
(1)
S j = å (Yi -Y )E(Gij | Dij ),T j =
i=1
var(S j )
Calculation of E(Gij|Dij)
• The calculation of the expected genotype given the sequencing data is
2
given by
E(Gij | Dij ) = å gP(Gij = g | Dij )
g=0
Where
 The genotype likelihood P(Dij | Gij ) 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 P(Gij = g) 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: Var (Gij ) = Var ( E (Gij | Dij )) + E (Var (Gij | Dij ))
 At high read depth, Var ( E (Gij | Dij )) ® Var (Gij ) because:

P(Gij  g | Dij )
goes to 1 for true genotype in
 And E (Var (G
ij
)
| Dij ) ® 0
 When read depth is not
high enough
E (Var (G | D )) > 0
ij
(
ij
)
Var E (Gij | Dij ) < Var (Gij )
Therefore, Var ( E (Gij | Dij )) is read depth dependent
E  Gij |Dij    gP(Gij  g | Dij )
2
g 0
Robust Variance Estimation
2 ˆ
2 ˆ
ˆ
Var(S
)
=
(1-Y
)
Var
(E(G
|
D
))
+
(Y
)
Varcontrols (E(Gij | Dij ))
å
å
j
case
ij
ij
cases
controls
S 2j
Tj =
~ c12
Var(S j )
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 Ncases where, for Ncontrols >> Ncases ,
Var(Sj) is underestimated
• Thus variance estimation must distinguish between the two groups
ˆ
é
ù
ˆ
• We estimate Var
Cases ë E ( Gij | Dij )û by Varcase (Gij )with estimated
genotype frequencies P(Gij = g)
ˆ
é
ù
• And Var
Controls ë E ( Gij | Dij )û 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: T j 
S 2j
ˆ (S )
Var
j
 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 assess Power 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
Method
CAST
Type of
analysis
Sample
Level of the test
Size
0.05 10-2 10-3 10-4
RVS
500 case 0.96 0.89 0.72 0.51
True
Genotypes
1500
0.97 0.92 0.79 0.61
controls
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
SNP
Position*
P-value (rank)
using Genotype
Calls:
RE versus Complete
Genomics controls
rs6484504
31424823
0.00008 (1)
0.0010 (3)
rs578666
31404484
0.0002 (2)
0.0003 (2)
rs674035
31399014
0.0007 (3)
0.0002 (1)
rs11031375
31428184
0.003 (4)
0.009 (7)
rs662702
31809070
0.012 (5)
0.052 (191)
rs11031330
31275073
0.011 (7)
0.0018(4)
rs603202
31866585
NA
0.0024(5)
* Build 37; approximately 450 variants analyzed
P-value (rank)
Using RVS:
RE versus 1000
Genomes controls
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
Derkach A, 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
 Andriy Derkach, graduate student in statistics at
University of Toronto
 Deb Pal, Richard Houlston and Ian Tomlinson