seqMeta (version 1.6.7)

skatMeta: Combine SKAT analyses from one or more studies

Description

Takes as input `seqMeta` objects (from the prepScores function), and meta-analyzes them.

Usage

skatMeta(..., SNPInfo = NULL, wts = function(maf) {     stats::dbeta(maf, 1, 25) }, method = "saddlepoint", snpNames = "Name", aggregateBy = "gene", mafRange = c(0, 0.5), verbose = FALSE)

Arguments

...

seqMeta 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.

wts

Either a function to calculate testing weights, or a character specifying a vector of weights in the SNPInfo file. For skatMeta the default are the `beta' weights.

method

p-value calculation method. Default is 'saddlepoint', 'integration' is the Davies method used in the SKAT package. See pchisqsum() for more details.

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 seqMeta objects.

mafRange

Range of MAF's to include in the analysis (endpoints included). Default is all SNPs (0

verbose

logical. Whether progress bars should be printed.

Value

a data frame with the following columns:

Details

skatMeta implements an efficient SKAT meta analysis by meta-analyzing scores statistics and their variances. Note: all studies must use coordinated SNP Info files - that is, the SNP names and gene definitions must be the same. Please see the package vignette for more details.

References

Wu, M.C., Lee, S., Cai, T., Li, Y., Boehnke, M., and Lin, X. (2011) Rare Variant Association Testing for Sequencing Data Using the Sequence Kernel Association Test (SKAT). American Journal of Human Genetics.

See Also

prepScores burdenMeta singlesnpMeta skatOMeta

Examples

###load example data for two studies:
### see ?seqMetaExample	
data(seqMetaExample)

####run on each study:
cohort1 <- prepScores(Z=Z1, y~sex+bmi, SNPInfo=SNPInfo, data=pheno1)
cohort2 <- prepScores(Z=Z2, y~sex+bmi, SNPInfo=SNPInfo, kins=kins, data=pheno2)

#### combine results:
##skat
out <- skatMeta(cohort1, cohort2, SNPInfo = SNPInfo)
head(out)

## Not run: 
# ##T1 test
# out.t1 <- burdenMeta(cohort1,cohort2, SNPInfo = SNPInfo, mafRange = c(0,0.01))
# head(out.t1)
# 
# ##single snp tests:
# out.ss <- singlesnpMeta(cohort1,cohort2, SNPInfo = SNPInfo)
# head(out.ss)
# 
# ########################
# ####binary data
# 
# cohort1 <- prepScores(Z=Z1, ybin~1, family=binomial(), SNPInfo=SNPInfo, data=pheno1)
# out.bin <- skatMeta(cohort1, SNPInfo=SNPInfo)
# head(out.bin)
# 
# ####################
# ####survival data
# cohort1 <- prepCox(Z=Z1, Surv(time,status)~strata(sex)+bmi, SNPInfo=SNPInfo, data=pheno1)
# out.surv <- skatMeta(cohort1, SNPInfo=SNPInfo)
# head(out.surv)
# 
# ##### Compare with SKAT on full data set
# require(SKAT)
# n <- nrow(pheno1)
# 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")])
# 
# obj <- SKAT_Null_Model(y~sex+bmi+gl(2,nrow(pheno1)), data=pheno)
# skat.pkg.p <- c(by(SNPInfo$Name, SNPInfo$gene, function(snp.names) {
#            inds <- match(snp.names,colnames(bigZ))
#            if(sum(!is.na(inds)) ==0 ) return(1)
#            SKAT(bigZ[,na.omit(inds)],obj, is_check=TRUE, missing=1)$p.value
#            }))
# 
# head(cbind(out$p,skat.pkg.p))
# 
# #Note: SKAT ignores family strucutre, resulting in p-values that are systematically too small: 
# plot(y=out$p,x=skat.pkg.p, ylab = "SKAT meta p-values", xlab = "SKAT p-values")
# abline(0,1)
# 
# ignore family structure:
# cohort1 <- prepScores(Z=Z1, y~sex+bmi, SNPInfo=SNPInfo, data=pheno1)
# cohort2 <- prepScores(Z=Z2, y~sex+bmi, SNPInfo=SNPInfo, data=pheno2)
# 
# out.nofam <- skatMeta(cohort1,cohort2,SNPInfo=SNPInfo)
# plot(y=out.nofam$p,x=skat.pkg.p, ylab = "SKAT meta p-values", xlab = "SKAT p-values")
# abline(0,1)
# ## End(Not run)

Run the code above in your browser using DataCamp Workspace