| Title: | GPA (Genetic analysis incorporating Pleiotropy and Annotation) |
|---|---|
| Description: | This package provides functions for fitting GPA, a statistical framework to prioritize GWAS results by integrating pleiotropy information and annotation data. In addition, it also includes ShinyGPA, an interactive visualization toolkit to investigate pleiotropic architecture. |
| Authors: | Dongjun Chung, Emma Kortemeier, Carter Allen |
| Maintainer: | Dongjun Chung <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.19.0 |
| Built: | 2024-11-29 05:56:36 UTC |
| Source: | https://github.com/bioc/GPA |
This package provides functions for fitting GPA, a statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data, along with ShinyGPA, a visualization toolkit to investigate the pleiotropic architecture using GWAS results.
| Package: | GPA |
| Type: | Package |
| Version: | 0.99.9 |
| Date: | 2020-04-14 |
| License: | GPL (>= 2) |
| LazyLoad: | yes |
This package contains a main class, GPA, which represents GPA model fit.
This package contains four main methods for the GPA framework (Chung et al., 2014),
GPA, assoc, pTest, and aTest.
GPA method fits the GPA model
and assoc method implements association mapping.
pTest and aTest methods implement hypothesis testing
for pleiotropy and annotation enrichment, respectively.
This package contains two main methods for the ShinyGPA visualization toolkit (Kortemeier et al., 2017), fitAll and shinyGPA. fitAll function generates all the intermediate results needed to run shinyGPA. shinyGPA opens the ShinyGPA interface, which takes the results generated from fitAll as input.
Dongjun Chung, Emma Kortemeier
Maintainer: Dongjun Chung <[email protected]>
Chung D*, Yang C*, Li C, Gelernter J, and Zhao H (2014), "GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data," PLoS Genetics, 10: e1004787. (* joint first authors)
Kortemeier E, Ramos PS, Hunt KJ, Kim HJ, Hardiman G, and Chung D (2018), "ShinyGPA: An interactive and dynamic visualization toolkit for genetic studies," PLOS One, 13(1): e0190949.
GPA, assoc, pTest, aTest,
GPA, fitAll, shinyGPA.
## Not run: # simulation setting nBin <- 1000 pi1 <- 0.2 common <- 0.5 betaAlpha <- c( 0.6, 0.6 ) annP <- c( 0.2, 0.4, 0.4, 0.4 ) seed <- 12345 # simulation setting nCommon <- round( pi1 * common * nBin ) nUniq <- round( pi1 * ( 1 - common ) * nBin ) nBg <- nBin - 2 * nUniq - nCommon # M * K matrix of GWAS p-value set.seed( seed ) pvec1 <- c( rbeta( nCommon, betaAlpha[1], 1 ), rbeta( nUniq, betaAlpha[1], 1 ), runif( nUniq ), runif( nBg ) ) pvec2 <- c( rbeta( nCommon, betaAlpha[2], 1 ), runif( nUniq ), rbeta( nUniq, betaAlpha[2], 1 ), runif( nBg ) ) pmat <- cbind( pvec1, pvec2 ) # M * D matrix of annotation ann <- c( sample( c(1,0), nCommon, replace=TRUE, prob = c( annP[4], 1 - annP[4] ) ), sample( c(1,0), nUniq, replace=TRUE, prob = c( annP[2], 1 - annP[2] ) ), sample( c(1,0), nUniq, replace=TRUE, prob = c( annP[3], 1 - annP[3] ) ), sample( c(1,0), nBg, replace=TRUE, prob = c( annP[1], 1 - annP[1] ) ) ) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL ) cov.GPA.noAnn <- cov( fit.GPA.noAnn ) # GPA with annotation data fit.GPA.wAnn <- GPA( pmat, ann ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) # GPA under pleiotropy H0 fit.GPA.pleiotropy.H0 <- GPA( pmat, NULL, pleiotropyH0=TRUE ) # association mapping assoc.GPA.wAnn <- assoc( fit.GPA.wAnn, FDR=0.05, fdrControl="global" ) # hypothesis testing for pleiotropy test.pleiotropy <- pTest( fit.GPA.noAnn, fit.GPA.pleiotropy.H0 ) # hypothesis testing for annotation enrichment test.annotation <- aTest( fit.GPA.noAnn, fit.GPA.wAnn ) # simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.4 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # Fit GPA for all possible pairs of GWAS datasets out <- fitAll( pmat ) # Run the ShinyGPA app using the ouput from fitAll() shinyGPA(out) ## End(Not run)## Not run: # simulation setting nBin <- 1000 pi1 <- 0.2 common <- 0.5 betaAlpha <- c( 0.6, 0.6 ) annP <- c( 0.2, 0.4, 0.4, 0.4 ) seed <- 12345 # simulation setting nCommon <- round( pi1 * common * nBin ) nUniq <- round( pi1 * ( 1 - common ) * nBin ) nBg <- nBin - 2 * nUniq - nCommon # M * K matrix of GWAS p-value set.seed( seed ) pvec1 <- c( rbeta( nCommon, betaAlpha[1], 1 ), rbeta( nUniq, betaAlpha[1], 1 ), runif( nUniq ), runif( nBg ) ) pvec2 <- c( rbeta( nCommon, betaAlpha[2], 1 ), runif( nUniq ), rbeta( nUniq, betaAlpha[2], 1 ), runif( nBg ) ) pmat <- cbind( pvec1, pvec2 ) # M * D matrix of annotation ann <- c( sample( c(1,0), nCommon, replace=TRUE, prob = c( annP[4], 1 - annP[4] ) ), sample( c(1,0), nUniq, replace=TRUE, prob = c( annP[2], 1 - annP[2] ) ), sample( c(1,0), nUniq, replace=TRUE, prob = c( annP[3], 1 - annP[3] ) ), sample( c(1,0), nBg, replace=TRUE, prob = c( annP[1], 1 - annP[1] ) ) ) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL ) cov.GPA.noAnn <- cov( fit.GPA.noAnn ) # GPA with annotation data fit.GPA.wAnn <- GPA( pmat, ann ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) # GPA under pleiotropy H0 fit.GPA.pleiotropy.H0 <- GPA( pmat, NULL, pleiotropyH0=TRUE ) # association mapping assoc.GPA.wAnn <- assoc( fit.GPA.wAnn, FDR=0.05, fdrControl="global" ) # hypothesis testing for pleiotropy test.pleiotropy <- pTest( fit.GPA.noAnn, fit.GPA.pleiotropy.H0 ) # hypothesis testing for annotation enrichment test.annotation <- aTest( fit.GPA.noAnn, fit.GPA.wAnn ) # simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.4 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # Fit GPA for all possible pairs of GWAS datasets out <- fitAll( pmat ) # Run the ShinyGPA app using the ouput from fitAll() shinyGPA(out) ## End(Not run)
Association mapping.
assoc( object, ... ) ## S4 method for signature 'GPA' assoc( object, FDR=0.05, fdrControl="global", pattern=NULL )assoc( object, ... ) ## S4 method for signature 'GPA' assoc( object, FDR=0.05, fdrControl="global", pattern=NULL )
object |
GPA model fit. |
FDR |
FDR level. |
fdrControl |
Method to control FDR. Possible values are "global" (global FDR control) and "local" (local FDR control). Default is "global". |
pattern |
Pattern for association mapping.
By default (i.e., |
... |
Other parameters to be passed through to generic |
assoc uses the direct posterior probability approach of Newton et al. (2004)
to control global FDR in association mapping.
Users can specify the pattern using 1 and * in pattern argument,
where 1 and * indicate phenotypes of interest and phenotypes that are not of interest, respectively.
For example, when there are three phenotypes,
pattern="111" means a SNP associated with all of three phenotypes,
while pattern="11*" means a SNP associated with the first two phenotypes
(i.e., association with the third phenotype is ignored (averaged out)).
If pattern=NULL, returns a binary matrix indicating association of SNPs for each phenotype,
where its rows and columns match those of input p-value matrix for function GPA.
Otherwise, returns a binary vector indicating association of SNPs for the phenotype combination of interest.
Dongjun Chung
Chung D*, Yang C*, Li C, Gelernter J, and Zhao H (2014), "GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data," PLoS Genetics, 10: e1004787. (* joint first authors)
Newton MA, Noueiry A, Sarkar D, and Ahlquist P (2004), "Detecting differential gene expression with a semiparametric hierarchical mixture method," Biostatistics, Vol. 5, pp. 155-176.
# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) fit.GPA.wAnn <- GPA( pmat, ann , maxIter = 100 ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) assoc.GPA.wAnn <- assoc( fit.GPA.wAnn, FDR=0.05, fdrControl="global" )# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) fit.GPA.wAnn <- GPA( pmat, ann , maxIter = 100 ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) assoc.GPA.wAnn <- assoc( fit.GPA.wAnn, FDR=0.05, fdrControl="global" )
Hypothesis testing for annotation enrichment.
aTest( fitWithoutAnn, fitWithAnn, vDigit=1000 )aTest( fitWithoutAnn, fitWithAnn, vDigit=1000 )
fitWithoutAnn |
GPA model fit without using annotation data. |
fitWithAnn |
GPA model fit with using annotation data. |
vDigit |
Number of digits for reporting parameter estimates and standard errors. For example, setting it to 1000 means printing out values up to three digits below zero. |
aTest implements the hypothesis testing for annotation enrichment.
It requires two GPA model fits,
one fitted with using annotation data and one fitted without using annotation data,
and evaluates annotation enrichment for risk-associated SNPs using the likelihood ratio test.
Returns a list with components:
q |
q estimates. |
statistics |
Statistics of the test for annotation enrichment. |
pvalue |
p-value of the test for annotation enrichment. |
Dongjun Chung
Chung D*, Yang C*, Li C, Gelernter J, and Zhao H (2014), "GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data," PLoS Genetics, 10: e1004787. (* joint first authors)
# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL, maxIter = 100 ) # GPA with annotation data fit.GPA.wAnn <- GPA( pmat, ann, maxIter = 100 ) # hypothesis testing for annotation enrichment test.annotation <- aTest( fit.GPA.noAnn, fit.GPA.wAnn )# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL, maxIter = 100 ) # GPA with annotation data fit.GPA.wAnn <- GPA( pmat, ann, maxIter = 100 ) # hypothesis testing for annotation enrichment test.annotation <- aTest( fit.GPA.noAnn, fit.GPA.wAnn )
Fit GPA model and the GPA model under H0 for all possible pairs of GWAS datasets.
fitAll( pmat, maxIter=2000, stopping="relative", epsStopLL=1e-10, parallel=FALSE, nCore=8 )fitAll( pmat, maxIter=2000, stopping="relative", epsStopLL=1e-10, parallel=FALSE, nCore=8 )
pmat |
p-value matrix from GWAS data, where row and column correspond to SNP and phenotype, respectively. |
maxIter |
Maximum number of EM iteration. Default is 2000. |
stopping |
Stopping rule for EM iteration.
Possible values are |
epsStopLL |
Threshold to stop the EM iteration. Default is 1e-100. |
parallel |
Utilize multiple CPUs for parallel computing
using |
nCore |
Number of CPUs when parallel computing is utilized. |
fitAll function fits the GPA model and the GPA model under H0 for all possible pairs of GWAS datasets. Its output can be used as an input for the shinyGPA function.
A list with 6 elements, including
pmat (original GWAS p-value matrix),
combs (a matrix of GWAS pair indices),
combList (a matrix of GWAS pair indices),
pTestPval (a matrix of pleiotropy test p-values),
fitGPA (a list of the GPA fit for each pair), and
fitH0 (a list of the GPA fit under H0 for each pair).
Dongjun Chung, Emma Kortemeier
Kortemeier E, Ramos PS, Hunt KJ, Kim HJ, Hardiman G, and Chung D (2017), "ShinyGPA: An interactive and dynamic visualization toolkit for genetic studies."
# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.4 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # Fit GPA for all possible pairs of GWAS datasets out <- fitAll( pmat, maxIter = 100 )# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.4 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # Fit GPA for all possible pairs of GWAS datasets out <- fitAll( pmat, maxIter = 100 )
Fit GPA model.
GPA( gwasPval, annMat=NULL, pleiotropyH0=FALSE, empiricalNull=FALSE, maxIter=2000, stopping="relative", epsStopLL=1e-10, initBetaAlpha=0.1, initPi=0.1, initQ=0.75, lbPi=NA, lbBetaAlpha=0.001, lbQ=0.001, lbPval=1e-30, vDigit=1000, verbose=1 )GPA( gwasPval, annMat=NULL, pleiotropyH0=FALSE, empiricalNull=FALSE, maxIter=2000, stopping="relative", epsStopLL=1e-10, initBetaAlpha=0.1, initPi=0.1, initQ=0.75, lbPi=NA, lbBetaAlpha=0.001, lbQ=0.001, lbPval=1e-30, vDigit=1000, verbose=1 )
gwasPval |
p-value matrix from GWAS data, where row and column correspond to SNP and phenotype, respectively. |
annMat |
Binary matrix from annotation data, where row and column correspond to SNP and annotation, respectively. |
pleiotropyH0 |
Fit GPA under the null hypothesis of pleiotropy test?
Possible values are |
empiricalNull |
Empirically estimate null distribution for GPA?
Possible values are |
maxIter |
Maximum number of EM iteration. Default is 2000. |
stopping |
Stopping rule for EM iteration.
Possible values are |
epsStopLL |
Threshold to stop the EM iteration. Default is 1e-100. |
initBetaAlpha |
Initial value for alpha estimate. Default is 0.1. |
initPi |
Initial value for pi estimate. Default is 0.1. |
initQ |
Initial value for q estimate. Default is 0.75. |
lbPi |
Lower bound for pi estimate.
If |
lbBetaAlpha |
Lower bound for alpha estimate. Default is 0.001. |
lbQ |
Lower bound for q estimate. Default is 0.001. |
lbPval |
Lower bound for GWAS p-value.
Any GWAS p-values smaller than |
vDigit |
Number of digits for reporting parameter estimates. For example, setting it to 1000 means printing out values up to three digits below zero. |
verbose |
Amount of progress report during the fitting procedure. Possible values are 0 (minimal output), 1, 2, or 3 (maximal output). Default is 1. |
GPA fits the GPA model. It requires to provide GWAS p-value to gwasPval,
while users can also provide annotation data to annMat.
It is assumed that number of rows of matrix provided to gwasPval
equals to that provided to annMat.
pTest implements the hypothesis testing for pleiotropy.
It requires two GPA model fits, one of interest and one under the null hypothesis,
and they can be obtained by setting pleiotropyH0=FALSE and pleiotropyH0=TRUE,
respectively.
aTest implements the hypothesis testing for annotation enrichment.
It requires two GPA model fits,
one fitted with using annotation data and one fitted without using annotation data,
and they can be obtained by providing annotation data to annMat and not, respectively.
Construct GPA class object.
Dongjun Chung
Chung D*, Yang C*, Li C, Gelernter J, and Zhao H (2014), "GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data," PLoS Genetics, 10: e1004787. (* joint first authors)
Kortemeier E, Ramos PS, Hunt KJ, Kim HJ, Hardiman G, and Chung D (2018), "ShinyGPA: An interactive and dynamic visualization toolkit for genetic studies," PLOS One, 13(1): e0190949.
# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL, maxIter = 100 ) cov.GPA.noAnn <- cov( fit.GPA.noAnn ) # GPA with annotation data fit.GPA.wAnn <- GPA( pmat, ann, maxIter = 100 ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) # GPA under the null hypothesis of pleiotropy test fit.GPA.pleiotropy.H0 <- GPA( pmat, NULL, pleiotropyH0=TRUE, maxIter = 100 )# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL, maxIter = 100 ) cov.GPA.noAnn <- cov( fit.GPA.noAnn ) # GPA with annotation data fit.GPA.wAnn <- GPA( pmat, ann, maxIter = 100 ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) # GPA under the null hypothesis of pleiotropy test fit.GPA.pleiotropy.H0 <- GPA( pmat, NULL, pleiotropyH0=TRUE, maxIter = 100 )
This class represents GPA model fit.
When users use fdr method, users can specify the pattern using 1 and * in pattern argument,
where 1 and * indicate phenotypes of interest and phenotypes that are not of interest, respectively.
For example, when there are three phenotypes,
pattern="111" means a SNP associated with all of three phenotypes,
while pattern="11*" means a SNP associated with the first two phenotypes
(i.e., association with the third phenotype is ignored (averaged out)).
Objects can be created by calls of the form new("GPA", ...).
fit:Object of class "list",
representing the fitted GPA model.
setting:Object of class "list",
representing the setting for GPA model fitting.
gwasPval:Object of class "matrix",
representing the p-value matrix from GWAS data.
annMat:Object of class "matrix",
representing the annotation matrix.
signature(object = "GPA"): provide brief summary of the object.
signature(x = "GPA"):
provide the matrix of posterior probability that a SNP belongs to each combination of association status.
signature(object = "GPA", pattern=NULL): provide local FDR.
By default (i.e., pattern=NULL),
it returns a matrix of local FDR that a SNP is not associated with each phenotype (i.e., marginal FDR),
where the order of columns is same as that in input GWAS data.
If a pattern is specified, a vector of corresponding local FDR is provided.
See the details about how users can specify the pattern.
signature(object = "GPA", silent=FALSE, vDigitEst=1000, vDigitSE=1000 ):
provide the covariance matrix for parameter estimates of GPA model.
If silent=TRUE, it suppresses the summary output.
vDigitEst and vDigitSE control number of digits for reporting parameter estimates and standard errors. For example, setting it to 1000 means printing out values up to three digits below zero.
signature(object = "GPA"):
extract parameter estimates from GPA model fit.
signature(object = "GPA"):
extract standard errors for parameter estimates from GPA model fit.
Dongjun Chung
Chung D*, Yang C*, Li C, Gelernter J, and Zhao H (2014), "GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data," PLoS Genetics, 10: e1004787. (* joint first authors)
showClass("GPA") # simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) fit.GPA.wAnn <- GPA( pmat, ann ) fit.GPA.wAnn pp.GPA.wAnn <- print( fit.GPA.wAnn ) fdr.GPA.wAnn <- fdr( fit.GPA.wAnn ) fdr11.GPA.wAnn <- fdr( fit.GPA.wAnn, pattern="11" ) fdr1..GPA.wAnn <- fdr( fit.GPA.wAnn, pattern="1*" ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) est.GPA.wAnn <- estimates( fit.GPA.wAnn ) se.GPA.wAnn <- se( fit.GPA.wAnn )showClass("GPA") # simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) ann <- rbinom(n = nrow(pmat), size = 1, prob = 0.15) ann <- as.matrix(ann,ncol = 1) fit.GPA.wAnn <- GPA( pmat, ann ) fit.GPA.wAnn pp.GPA.wAnn <- print( fit.GPA.wAnn ) fdr.GPA.wAnn <- fdr( fit.GPA.wAnn ) fdr11.GPA.wAnn <- fdr( fit.GPA.wAnn, pattern="11" ) fdr1..GPA.wAnn <- fdr( fit.GPA.wAnn, pattern="1*" ) cov.GPA.wAnn <- cov( fit.GPA.wAnn ) est.GPA.wAnn <- estimates( fit.GPA.wAnn ) se.GPA.wAnn <- se( fit.GPA.wAnn )
Hypothesis testing for pleiotropy.
pTest( fit, fitH0, vDigit=1000 )pTest( fit, fitH0, vDigit=1000 )
fit |
Fit of the GPA model of interest. |
fitH0 |
GPA model fit under the null hypothesis of pleiotropy test. |
vDigit |
Number of digits for reporting parameter estimates and standard errors. For example, setting it to 1000 means printing out values up to three digits below zero. |
pTest implements the hypothesis testing for pleiotropy.
It requires two GPA model fits, one of interest and one under the null hypothesis
(obtained by setting pleiotropyH0=TRUE when running GPA function),
and evaluates genetical correlation among multiple phenotypes using the likelihood ratio test.
Returns a list with components:
pi |
pi estimates. |
piSE |
Standard errors for pi estimates. |
statistics |
Statistics of the pleiotropy test. |
pvalue |
p-value of the pleiotropy test. |
Dongjun Chung
Chung D*, Yang C*, Li C, Gelernter J, and Zhao H (2014), "GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data," PLoS Genetics, 10: e1004787. (* joint first authors)
# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL, maxIter = 100 ) # GPA under the null hypothesis of pleiotropy test fit.GPA.pleiotropy.H0 <- GPA( pmat, NULL, pleiotropyH0=TRUE, maxIter = 100 ) # hypothesis testing for pleiotropy test.pleiotropy <- pTest( fit.GPA.noAnn, fit.GPA.pleiotropy.H0 )# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # GPA without annotation data fit.GPA.noAnn <- GPA( pmat, NULL, maxIter = 100 ) # GPA under the null hypothesis of pleiotropy test fit.GPA.pleiotropy.H0 <- GPA( pmat, NULL, pleiotropyH0=TRUE, maxIter = 100 ) # hypothesis testing for pleiotropy test.pleiotropy <- pTest( fit.GPA.noAnn, fit.GPA.pleiotropy.H0 )
Run ShinyGPA app.
shinyGPA( out=NULL )shinyGPA( out=NULL )
out |
output of |
shinyGPA runs the ShinyGPA app. It takes the output of the fitAll function, which fits the GPA model for all possible pairs of GWAS datasets, as input.
Provides visualization to investigate pleiotropic architecture using GWAS results.
Dongjun Chung, Emma Kortemeier
Kortemeier E, Ramos PS, Hunt KJ, Kim HJ, Hardiman G, and Chung D (2018), "ShinyGPA: An interactive and dynamic visualization toolkit for genetic studies," PLOS One, 13(1): e0190949.
# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # Fit GPA for all possible pairs of GWAS datasets out <- fitAll( pmat, maxIter = 100 ) # Run the ShinyGPA app using the ouput from fitAll() # shinyGPA(out)# simulator function simulator <- function( risk.ind, nsnp=20000, alpha=0.6 ) { m <- length(risk.ind) p.sig <- rbeta( m, alpha, 1 ) pvec <- runif(nsnp) pvec[ risk.ind ] <- p.sig return(pvec) } # run simulation set.seed(12345) nsnp <- 1000 alpha <- 0.3 pmat <- matrix( NA, nsnp, 5 ) pmat[,1] <- simulator( c(1:200), nsnp=nsnp, alpha=alpha ) pmat[,2] <- simulator( c(51:250), nsnp=nsnp, alpha=alpha ) pmat[,3] <- simulator( c(401:600), nsnp=nsnp, alpha=alpha ) pmat[,4] <- simulator( c(451:750), nsnp=nsnp, alpha=alpha ) pmat[,5] <- simulator( c(801:1000), nsnp=nsnp, alpha=alpha ) # Fit GPA for all possible pairs of GWAS datasets out <- fitAll( pmat, maxIter = 100 ) # Run the ShinyGPA app using the ouput from fitAll() # shinyGPA(out)