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 Score
## 1 1 1 970546 1 1126 0.0039964476 0.20256722
## 2 2 1 985900 1 1126 0.0492895204 -2.64169956
## 3 3 1 1025045 1 1126 0.0004440497 -0.09916904
## 4 4 1 1265550 1 1126 0.0008880995 0.81717324
## 5 5 1 1472676 1 1126 0.0071047957 0.64418361
## 6 6 1 1735725 1 1126 0.0022202487 -0.46319177
## Score.SE Score.Stat Score.pval
## 1 0.8351783 0.2425437 0.80835892
## 2 2.6522412 -0.9960254 0.31923781
## 3 0.2972472 -0.3336248 0.73866267
## 4 0.4033577 2.0259271 0.04277226
## 5 1.0778277 0.5976685 0.55006117
## 6 0.6396675 -0.7241134 0.46899613
- 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_0 pval_0
## 1 1 966001 971000 1 9 9 7.318094 0.643988703
## 2 1 982001 987000 1 111 107 154.178280 0.058811539
## 3 1 1022001 1027000 1 1 1 47.823916 0.005504472
## 4 1 1262001 1267000 1 2 2 7.319239 0.308472754
## 5 1 1468001 1473000 1 16 16 58.518662 0.293504072
## 6 1 1732001 1737000 1 5 5 9.499539 0.509868791
## err_0
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0