Association tests - Solutions
- Logistic regression:
fitNullModel can use a binary phenotype as the outcome variable by specifying the argument family=binomial. Use the status column in the sample annotation to fit a null model for simulated case/control status, with sex and Population as covariates. Then run a single-variant test using this model.
nullmod.status <- fitNullModel(annot, outcome="status", covars=c("sex", "Population"),
family=binomial, verbose=FALSE)
resetIterator(iterator, verbose=FALSE)
assoc <- assocTestSingle(iterator, nullmod.status, test="Score")
## # of selected samples: 1,126
## variant.id chr pos allele.index n.obs freq MAC Score
## 1 1 1 970546 1 1126 0.0039964476 9 0.20256722
## 2 2 1 985900 1 1126 0.0492895204 111 -2.64169956
## 3 3 1 1025045 1 1126 0.0004440497 1 -0.09916904
## 4 4 1 1265550 1 1126 0.0008880995 2 0.81717324
## 5 5 1 1472676 1 1126 0.0071047957 16 0.64418361
## 6 6 1 1735725 1 1126 0.0022202487 5 -0.46319177
## Score.SE Score.Stat Score.pval Est Est.SE PVE
## 1 0.8351783 0.2425437 0.80835892 0.2904095 1.1973491 5.242118e-05
## 2 2.6522412 -0.9960254 0.31923781 -0.3755410 0.3770396 8.840314e-04
## 3 0.2972472 -0.3336248 0.73866267 -1.1223819 3.3642035 9.918446e-05
## 4 0.4033577 2.0259271 0.04277226 5.0226565 2.4791892 3.657417e-03
## 5 1.0778277 0.5976685 0.55006117 0.5545121 0.9277921 3.183080e-04
## 6 0.6396675 -0.7241134 0.46899613 -1.1320153 1.5633123 4.672400e-04
- Inverse normal transform: use the function
nullModelInvNorm to perform an inverse normal transform on the height variable. Compare these residuals with the residuals from the original null model.
nullmod.norm <- nullModelInvNorm(nullmod, verbose=FALSE)
summary(nullmod$resid.marginal)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -38.9394 -7.6295 -0.1089 0.0000 8.0824 44.6166
summary(nullmod.norm$resid.marginal)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.0075251 -0.6717166 -0.0000291 0.0000000 0.6751482 3.0032466
- Repeat the previous exercise on logistic regression, this time running a sliding-window test.
nullmod.status <- fitNullModel(annot, outcome="status", covars=c("sex", "Population"),
family=binomial, verbose=FALSE)
seqResetFilter(seqData, verbose=FALSE)
iterator <- SeqVarWindowIterator(seqData, windowSize=5000, windowShift=2000, verbose=FALSE)
assoc <- assocTestAggregate(iterator, nullmod, test="SKAT", AF.max=0.1, weight.beta=c(1,25))
## # of selected samples: 1,126
## chr start end n.site n.alt n.sample.alt Q pval err
## 1 1 966001 971000 1 9 9 7.318095 0.643988693 0
## 2 1 982001 987000 1 111 107 154.178289 0.058811535 0
## 3 1 1022001 1027000 1 1 1 47.823918 0.005504472 0
## 4 1 1262001 1267000 1 2 2 7.319239 0.308472744 0
## 5 1 1468001 1473000 1 16 16 58.518665 0.293504065 0
## 6 1 1732001 1737000 1 5 5 9.499539 0.509868790 0
## pval.method
## 1 integration
## 2 integration
## 3 integration
## 4 integration
## 5 integration
## 6 integration