Description Usage Arguments Details Value Author(s) References See Also Examples
Takes as input 'skatCohort' objects (from the skatCohort
function), and meta analyzes them.
1 2 |
... |
skatCohort objects |
SNPInfo |
The SNP Info file. This should contain the fields listed in snpNames and aggregateBy. Only SNPs in this table will be meta analyzed, so this may be used to restrict the analysis. |
snpNames |
The field of SNPInfo where the SNP identifiers are found. Default is 'Name' |
aggregateBy |
The field of SNPInfo on which the skat results were aggregated. Default is 'gene'. Though gene groupings are not explicitely required for single snp analysis, it is required to find where single snp information is stored in the skatCohort objects. |
cohortBetas |
Whether or not to include cohort-level effects in the output. |
verbose |
logical. Whether progress bars should be printed. |
This function meta analyzes score tests for single variant effects. Though the test is formally a score test, coefficients and standard errors are also returned, which can be interpreted as a 'one-step' approximation to the maximum likelihood estimates.
a data frame with the gene, snp name, meta analysis.
Arie Voorman, Jennifer Brody
Lin, DY and Zeng, D. On the relative efficiency of using summary statistics versus individual-level data in meta-analysis. Biometrika. 2010.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 | data(skatExample)
####run on each cohort:
cohort1 <- skatCohort(Z=Z1, y~sex+bmi, SNPInfo = SNPInfo, data =pheno1)
cohort2 <- skatCohort(Z=Z2, y~sex+bmi, SNPInfo = SNPInfo, data =pheno2)
#### combine results:
out <- singlesnpMeta(cohort1, cohort2, SNPInfo = SNPInfo)
head(out)
## Not run:
##compare
bigZ <- matrix(NA,2*n,nrow(SNPInfo))
colnames(bigZ) <- SNPInfo$Name
for(gene in unique(SNPInfo$gene)){
snp.names <- SNPInfo$Name[SNPInfo$gene == gene]
bigZ[1:n,SNPInfo$gene == gene][,
snp.names %in% colnames(Z1)] <- Z1[,
na.omit(match(snp.names,colnames(Z1)))]
bigZ[(n+1):(2*n),SNPInfo$gene == gene][,
snp.names %in% colnames(Z2)] <- Z2[,
na.omit(match(snp.names,colnames(Z2)))]
}
pheno <- rbind(pheno1[,c("y","sex","bmi")],pheno2[,c("y","sex","bmi")])
out3 <- apply(bigZ,2,function(z){
if(any(!is.na(z))){
z[is.na(z)] <- mean(z,na.rm=TRUE)
mod <- lm(y ~ sex+bmi+c(rep(1,n),rep(0,n))+z,data=pheno)
beta <- mod$coef["z"]
se <- sqrt(vcov(mod)["z","z"])
p <- pchisq( (beta/se)^2,df=1,lower=F)
return(c(beta,se,p))
} else {
return(c(0,Inf,1))
}
})
out3 <- t(out3[,match(out$Name,colnames(out3))])
##plot
par(mfrow=c(2,2))
plot(x=out3[,1],y=out$beta, xlab = "complete data (Wald)",
ylab = "meta-analysis (Score)", main = "coefficients");abline(0,1)
plot(x=out3[,2],y=out$se, xlab = "complete data (Wald)", y
lab = "meta-analysis (Score)", main = "standard errors");abline(0,1)
plot(x=out3[,3],y=out$p, xlab = "complete data (Wald)",
ylab = "meta-analysis (Score)", main = "p-values");abline(0,1)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.