Quantitative Genomics and Genetics BTRY 4830/6830; PBSB.5201.01 - - PowerPoint PPT Presentation
Quantitative Genomics and Genetics BTRY 4830/6830; PBSB.5201.01 Lecture19: Alternative Tests, Haplotype Testing, and Minimal GWAS Steps Jason Mezey jgm45@cornell.edu April 16, 2019 (T) 10:10-11:25 Summary of lecture 19 Today we will
Jason Mezey jgm45@cornell.edu April 16, 2019 (T) 10:10-11:25
(GLMs) the class that includes linear and logistic regression
approaches in GWAS (i.e., other than regressions)
how haplotypes are defined
when performing a GWAS
BTRY 4830/6830 & PBSB.5201.01 Quantitative Genomics and Genetics Spring 2019
Project (Version 1) Posted April 16; Due 11:59PM May 7
1 Introduction and instructions
The goal of the class project is for you to demonstrate what you have learned by performing a GWAS analysis on real data. To accomplish this, assume that you have been provided data by a collaborator who wants to identify positions of causal polymorphisms (loci). You will perform an in-depth analysis and write a report for your collaborator that explains your methods and results. Instructions: While we provide some general guidelines for how to proceed below, the techniques you use to analyze the data and how you construct your report will be up to you. Do however note the following instructions (PLEASE READ THESE CAREFULLY!!):
(1) Your project must be uploaded by 11:59PM, May 7 - if it is late for any reason, standard grading policies apply. (2) You are allowed to work together with other students in the class to analyze these data. However, note that turning in a report that describes exactly the same analyses as a fellow student is not a good strategy for getting a good grade. Also note that you must write your
(3) This is an ‘open book’ assignment, such that you are allowed to use any resources online, in books, etc. You may also ask third-party (i.e. people not in the class) for suggestions on what analyses to perform but you cannot have a third-party do any of the analyses (or write any code for you!). (4) You are also allowed to use any software or programming language that you would like as part of your analysis. However, we expect that some of the tasks will be performed in R (also note that you are welcome to use any packages, functions, etc. in R). (5) Your final project will include a SINGLE report file and a SINGLE file including all of your R code (ideally an .rmd file!) and / or commands or scripts you used to run other software
commented code that we can run from the command line. If you use other software for some
that provides details on how you ran the software, e.g. what parameters did you use, etc.
(6) The report file must be no more than 8 pages (single-sided), with NO MORE than 5 pages
(7) For your report, you must describe what you did in detail (a good guide is have you provided enough detail such that someone reading your report could replicate what you have done?). You also need to describe the results you have obtained from your analysis. You may also wish to include some text to describe interpretations and conclusions that may be of interest to your collaborator, including statistical and possibly, biological interpretations. For your Figures and Tables, note that clarity and clear labels is a strategy for maximizing your grade. (8) We will grade on two broad criteria: 1. the overall quality of the analyses / report, 2. the amount of effort put into your project. Note that ‘effort’ does not mean run many analyses without thinking carefully about why you are running them or how they fit together to provide a clear picture of results. A guide maximizing your grade on effort is to think carefully about how to produce the best possible report that you can and then put in as many hours as you wish to devote to the project accomplishing this objective (your effort level will be clear to us).
2 The experiment and data
The experiment: Among the recent large scale human genomics resources is Genetic European Variation in Health and Disease (gEUVADIS): http://www.geuvadis.org/ with a samples from 4 different European populations. Each of these individuals were part of the 1000 Genomes project and their genomes were sequenced and analyzed to identify SNP geno-
and mRNA levels were quantified through RNA sequencing. Each of these gene expression measurements may be thought of as a phenotype and one can do a GWAS analysis on each individually, which is called an ?expression Quantitative Trait locus? or eQTL analysis, an unnecessarily fancy name for a GWAS when the phenotype is gene expression. What you have been provided is a small subset of these data that are publicly available. Specifi- cally, you have been provided 50,000 of the SNP genotypes for 344 samples from the CEU (Utah residents with European ancestry), FIN (Finns), GBR (British) and, TSI (Toscani) population. For these same individuals, you have also been provided the expression levels of five genes. You have also been provided information on the population and gender of each of these individuals, and information regarding the position of each gene and SNP in the genome. A description of the broader data set from which these data were extracted can be found in:: http://www.geuvadis.org/web/geuvadis/RNAseq-project and in other papers relating to analysis of the GEUVADIS data.
The data: These have been provided to you in five total files: ‘phenotypes.csv’,‘genotypes.csv’,
‘covars.csv’, ‘gene info.csv’,‘SNP info.csv’. ‘phenotypes.csv’ contains the phenotype data for 344 samples and 5 genes. ‘genotypes.csv’ contains the SNP data for 344 samples and 50000 genotypes. ‘gene info.csv’ contains information about each gene that was measured. The ‘chromosome’ column indicates the chromosome where the gene is located, ‘start’ marks the position in the chromosome where the region of the gene begins and ‘end’ marks the position where the region ends, ‘symbol’ contains the common gene name of the measured transcript and ‘probe’ contains the ids of the transcripts that match with the column names of the phenotype data. ‘SNP info.csv’ contains the additional information on the genotypes and has four columns. The 1st column contains the chromosome number of each SNP, the 2nd column contains the physical position of the SNP on the chromosome, the 3rd column contains the abbreviation used to the ‘rsID’ = the name of each SNP in order.
3 Your assignment and hints for getting started
Your GWAS assignment is to find the position of as many causal polymorphisms as possible for the five expressed genes using the data (note that each ‘hit’ will potentially indicate an eQTL). You may / should use any and as many analysis approaches as you think that are useful to accomplish this goal. In your report, you will need to describe in detail what you did, why you did it, and describe results in a manner that your ‘non-statistical’ collaborator will be able to understand, e.g. explain your terms, provide interpretations, etc. A few hints:
browser, dbSNP, many others) to incorporate biological information into your interpretation and hypotheses concerning what you may have found.
Good luck!
with a case / control phenotype!
individuals where for each we have measured a case / control phenotype and N genotypes, we perform N hypothesis tests
algorithm (twice!) for EACH marker to get the MLE of the parameters under the alternative (= no restrictions on the beta’s!) and under the null (= null hypothesis parameters set to zero!) and use these to calculate our LRT test statistic for each marker
a chi-square distribution (how do we do this is R?)
H0 : a = 0 ∩ d = 0 HA : βa 6= 0 [ βd 6= 0
LRT = 2lnΛ = 2lnL(ˆ θ0|y) L(ˆ θ1|y)
θ1|y) 2l(ˆ θ0|y)
⌅1|y) = l( ˆ µ, ˆ a, ˆ d|y) l( ˆ ⌅0|y) = l( ˆ µ, 0, 0|y)
algorithm to find the parameters that correspond to the maximum
Iterative Re-weighted Least Squares (IRLS) algorithm, which has the following structure:
we assign these numbers and call the resulting vector β[0].
an appropriate function. If the value is below a defined threshold, stop. If not, repeat steps 2,3.
l(β) =
n
⇤
i=1
⇥
equations:
making use of the IRLS algorithm:
2 degree of freedom (why?):
LRT = −2lnΛ = 2l(ˆ θ1|y) − 2l(ˆ θ0|y)
d f=2
l(ˆ θ1|y) =
n
X
i=1
[yiln(γ−1(βµ + xi,aβa + xi,dβd))+(1−yi)ln(γ−1(βµ + xi,aβa + xi,dβd))] l(ˆ θ0|y) =
n
X
i=1
[yiln(γ−1(βµ)) + (1 − yi)ln(γ−1(βµ))]
| | ˆ θ0 = {ˆ βµ, ˆ βa = 0, ˆ βd = 0}
{ ˆ θ1 = {ˆ βµ, ˆ βa, ˆ βd}
distribution of our LRT statistic under the null hypothesis
n (contrast with F-statistics!!) but we know that as n goes to infinite, we know the distribution is i.e. ( ):
θ1|y) 2l(ˆ θ0|y)
ve an exact dis n LRT ! χ2
d f,
df) that depen
equations:
making use of the IRLS algorithm (just add an additional parameter!):
2 degree of freedom (why?):
LRT = −2lnΛ = 2l(ˆ θ1|y) − 2l(ˆ θ0|y) l(ˆ θ1|y) =
n
X
i=1
h yiln(γ−1(ˆ βµ + xi,a ˆ βa + xi,d ˆ βd + xi,z ˆ βz)) + (1 − yi)ln(1 − γ−1(ˆ βµ + xi,a ˆ βa + xi,d ˆ βd + xi,z ˆ βz)) i
d f=2
ˆ θ0 = {ˆ βµ, ˆ βa = 0, ˆ βd = 0, ˆ βz}
ˆ θ1 = {ˆ βµ, ˆ βa, ˆ βd, ˆ βz}
l(ˆ ✓0|y) =
n
X
i=1
yiln(−1(ˆ µ + xi,z ˆ z)) + (1 yi)ln(1 −1(ˆ µ + xi,z ˆ z))
more than one) we start with the same hypotheses:
(run the algorithm twice as before!)
H0 : a = 0 ∩ d = 0 HA : βa 6= 0 [ βd 6= 0
Yi = eβµ+xi,aβa+xi,dβd+xi,zβz 1 + eβµ+xi,aβa+xi,dβd+xi,zβz + ✏i
MLE(ˆ ) = MLE(ˆ µ, ˆ a, ˆ d, ˆ z)
LRT = 2lnΛ = 2l(ˆ θ1|y) 2l(ˆ θ0|y)
ve an exact dis n LRT ! χ2
d f,
df) that depen
l() =
n
X
i=1
[yiln(1(µ + xi,aa + xi,dd + xi,zz))+(1yi)ln(1(µ + xi,aa + xi,dd + xi,zz))] (27)
l(ˆ θ1|y) =
n
X
i=1
h yiln(γ−1(ˆ βµ + xi,a ˆ βa + xi,d ˆ βd + xi,z ˆ βz)) + (1 − yi)ln(1 − γ−1(ˆ βµ + xi,a ˆ βa + xi,d ˆ βd + xi,z ˆ βz)) i
l(ˆ ✓0|y) =
n
X
i=1
yiln(−1(ˆ µ + xi,z ˆ z)) + (1 yi)ln(1 −1(ˆ µ + xi,z ˆ z))
} Yi = −1(X) + ✏i
describe how linear and logistic models fit into this framework
them using three (models that have these properties are considered GLMs):
Y conditional on the independent variable X is in the exponential family of distributions
expected value of the response variable (where we often use the inverse!!)
. Pr(Y |X) ∼ expfamily.
: : E(Y|X) → X, (E(Y|X)) = X
E(Y|X) = −1(X)
= X
same approach, i.e. MLE of the beta parameters using an IRLS algorithm (just substitute the appropriate link function in the equations, etc.)
(where the sampling distribution as the sample size goes to infinite is chi-square)
types of regression modeling you will likely need to apply (!!)
GWAS covers a large number of possible relationships between genotypes and phenotypes, there are a large number
based tests, Kruskal-Wallis tests, other rank based tests, chi- square, Fisher’s exact, Cochran-Armitage, etc. (see PLINK for a somewhat comprehensive list of tests used in GWAS)
GWAS
Yes we should - the reason is different tests have different performance depending on the (unknown) conditions of the system and experiment, i.e. some may perform better than others
will be best suited) we should run a number of tests and compare results
conditional rules) but a reasonable approach is to treat each test as a distinct GWAS analysis and compare the hits across analyses using the following rules:
good evidence that there is a causal polymorphism
are not) we should attempt to determine the reason for the discrepancy (this requires that we understand the tests and experience)
review of possible tests (keep in mind, every time you learn a new test in a statistics class, there is a good chance you could apply it in a GWAS!)
applied
perform N alternative tests for each marker-phenotype combinations, where for each case, we are testing the following hypotheses with different (implicit) codings of X (!!):
H0 : Cov(Y, X) = 0 HA : Cov(Y, X) 6= 0
test (which has deep connections to our logistic regression test under certain assumptions but it has slightly different properties!)
phenotype combinations (left) and calculate the expected numbers in each cell (right):
in this c is χ2
d.f.=2.
ize tends to infinite, i.e. when the sam is d.f. = (#columns-1)(#rows-1) = 2, can therefore calculate the statistic in
LRT = −2lnΛ = −2
3
X
i=1 2
X
j=1
nijln ni n.inj. !
Case Control A1A1 n11 n12 n1. A1A2 n21 n22 n2. A2A2 n31 n32 n3. n.1 n.2 n
Case Control A1A1 (n.1n1.)/n (n.2n1.)/n n1. A1A2 (n.1n2.)/n (n.2n2.)/n n2. A2A2 (n.1n3.)/n (n.2n3.)/n n3. n.1 n.2 n
but it is not exact for smaller sample sizes (although we hope it is close!)
Yes, we can calculate a Fisher’s test statistic for our sample, where the distribution under the null hypothesis is exact for any sample size (I will let you look up how to calculate this statistic and the distribution under the null on your own):
for when we should use one versus the other?
Case Control A1A1 n11 n21 A1A2 n21 n22 A2A2 n31 n32 hi-square test) is also often
square or a Fisher’s exact test
(right):
PLINK):
Case Control A1A1 n11 n12 A1A2 ∪ A2A2 n21 n22 Case Control A1A1 ∪ A1A2 n11 n12 A2A2 n21 n22 Case Control A1 n11 n12 A2 n21 n22
following (now what!?):
correction regardless of the testing approach used, i.e. the result is robust to testing approach.
understand why test results could be different:
true? Are they handling missing data in different ways?
this be?
handle additional phenotypes
genotypes defined using a different approach
testing framework
information but in some cases, testing using haplotypes is more effective than testing one genetic marker at a time
are inherited together
“function” that takes a set of alleles at several loci A, B, C, D, etc. and outputs a haplotype allele:
following alleles (A,G), (A,T),(G,C),(G,C) we could define the following haplotype alleles:
h = f(Ai, Bj, ...)
s h1 = (A, A, C, C)
s (A, G), (A, T), (G, d h2 = (G, T, G, G).
general, we define a haplotype for a set of genetic markers (loci) that are physically linked that are frequently occur in a population
define haplotype alleles that have appreciable frequenecy in the population
many haplotype alleles could we define?
relatively “high” allele frequencies (say >0.05 = arbitrary!)
s h1 = (A, A, C, C)
s (A, G), (A, T), (G, d h2 = (G, T, G, G).
loci (A,G), (A,T),(G,C),(G,C), we have lots of LD (!!) such that there are only 4 alleles in the population (i.e. all other combinations have frequency of zero!):
the population are < 0.01
the other alleles as follows (where * means “any” genetic marker allele):
s h∗
1 = (A, A, C, C),h∗ 2 = (G, T, G, G),h∗ 3 = (A, A, G, C),h∗ 4 = (G, T, C, G)
s h1 = (A, A, ∗, C)
∗ ∗
d h2 = (G, T, ∗, G) d h = h∗ ∪ h∗. ∗ at h1 = h∗
1 ∪h∗ 3
d h2 = (G, T, ∗ d h2 = h∗
2 ∪ h∗ 4.
proceed with a GWAS using our framework (just substitute haplotype alleles and genotypes for genetic marker alleles and genotypes!)
haplotype alleles, we can code our independent variables for our regression model as follows:
effect on our interpretation of where the causal polymorphism is located?)
Xa(h1h1) = −1, Xa(h1h2) = 0, Xa(h2h2) = 1 Xd(h1h1) = −1, Xd(h1h2) = 1, Xd(h2h2) = −1
haplotype testing approach in a GWAS, why might we want to use this approach?
haplotypes in a population:
A1 B1 (C1)⇤ D2 E1 A1 B2 (C1)⇤ D1 E1 A2 B1 (C1)⇤ D1 E1 A1 B1 (C1)⇤ D1 E2
A2 B2 (C2)⇤ D1 E2 A2 B1 (C2)⇤ D2 E2 A1 B2 (C2)⇤ D2 E2 A2 B2 (C2)⇤ D2 E1
haplotype is a better “tag” of the causal polymorphism than any of the surrounding markers
therefore has a higher probability of correctly rejecting the null hypothesis
are performing less total tests in our GWAS (in what sense is this an advantage!?)
it also may not (!!), i.e. sometimes (in fact often!) individual genetic markers are better tags of the causal polymorphism
cannot resolve the position of the causal polymorphism to a position smaller than the range of the haplotype alleles, i.e. large haplotypes can have smaller resolution
data, should we use haplotype testing (i.e. in the future, the importance of haplotype testing may decrease)
(always!) as well as a haplotype test (optional)
GWAS (as with any statistical analysis!) so it doesn’t hurt us to explore our dataset with as many techniques as we want to apply
analyzing GWAS with as many methods as you find useful
if we get conflicting results, which one do we interpret as “correct”!?
fuzzy definition!) is another subject that is in the realm of population genetics and therefore we cannot discuss this in detail in this class (again: I encourage you to take a class on population genetics!)
come from by remembering that the origin of new haplotype alleles are mutations and that new haplotype alleles can be produced by recombination
in the population and therefore what blocks of alleles are inherited as a haplotype (and we therefore use them to define haplotypes using system specific criteria)
haplotypes for given systems and the algorithms used for this purpose (so we will just briefly mention the main concepts here)
genotype markers, decide on the number of genotype markers to put together into a haplotype block, and decide how many haplotype alleles to consider
(system dependent!)
where we have two markers that are right next to each
together in a haplotype block
we are considering a diploid individual who is a heterozygote for both of these markers, which of the marker alleles are physically linked in this individual?
problem and there are many algorithms for accomplishing this goal (note that in the future, technology may make this a non-issue...)
genotypes to include in a haplotype block depends on LD
each other but low LD with other markers, we use this as a guide for defining the haplotype block
0.4 0.6 0.8 5 10 15 0.2 0.4 0.6 0.8 1
RCOR3 TRAF5 C1orf97 RD3 SLC30A1 NEK2 LPGAT1
genes 30 60 1 2
211.5 211.6 211.7 211.8 211.9 212.0
* *
Single marker test Conditional test VBAY Lasso Adaptive Lasso 2D-MCP LOG NEG Independent hit Non−independent hit
10 15
genes
to define, but in general, we define a set where the frequency in a population is above some MAF threshold (which depends on the system)
haplotype alleles (e.g. in humans!)
rarer haplotypes (what are some of these?)
genetic systems work (just like LD!)
alleles, we can use them as we would genetic markers in our GWAS analysis framework
analysis of your GWAS in addition to your single marker analysis (ALWAYS do a single marker analysis)