pcs.model — pcs.model • coloc

Functions to prepare principle component models for colocalisation testing

pcs.model(
  object,
  group,
  Y,
  stratum = NULL,
  threshold = 0.8,
  family = if (all(Y %in% c(0, 1))) {     "binomial" } else {     "gaussian" }
)

Arguments

object	A colocPCs object, result of `pcs.prepare()`.
group	1 or 2, indicating which group of samples to extract from principal components matrix
Y	Numeric phenotype vector, length equal to the number of samples from the requested group
stratum	optional vector that gives stratum information
threshold	The minimum number of principal components which captures at least threshold proportion of the variance will be selected. Simulations suggest `threshold=0.8` is a good default value.
family	Passed to `glm()` function. `pcs.model` attempts to guess, either "binomial" if `Y` contains only 0s and 1s, "gaussian" otherwise.

Value

pcs.prepare returns a colocPCs object, pcs.model returns a glm object.

Details

Prepares models of response based on principal components of two datasets for colocalisation testing.

References

Wallace et al (2012). Statistical colocalisation of monocyte gene expression and genetic risk variants for type 1 diabetes. Hum Mol Genet 21:2815-2824. http://europepmc.org/abstract/MED/22403184

Plagnol et al (2009). Statistical independence of the colocalized association signals for type 1 diabetes and RPS26 gene expression on chromosome 12q13. Biostatistics 10:327-34. http://www.ncbi.nlm.nih.gov/pubmed/19039033

Examples


 ## simulate covariate matrix (X) and continuous response vector (Y)
 ## for two populations/triats Y1 and Y2 depend equally on f1 and f2
 ## within each population, although their distributions differ between
 ## populations.  They are compatible with a null hypothesis that they
 ## share a common causal variant, with the effect twice as strong for
 ## Y2 as Y1
set.seed(1)
 X1 <- matrix(rbinom(5000,1,0.4),ncol=10)
 Y1 <- rnorm(500,apply(X1[,1:2],1,sum),2)
 X2 <- matrix(rbinom(5000,1,0.6),ncol=10)
 Y2 <- rnorm(500,2*apply(X2[,1:2],1,sum),5)

 ## generate principal components object
 colnames(X1) <- colnames(X2) <- make.names(1:ncol(X1))
 pcs <- pcs.prepare(X1,X2)

 ## generate glm objects
 m1 <- pcs.model(pcs, group=1, Y=Y1)
#> selecting 8 components out of 10 to capture 0.8299584 of total variance.
 m2 <- pcs.model(pcs, group=2, Y=Y2)
#> selecting 8 components out of 10 to capture 0.8299584 of total variance.

 ## Alternatively, if one (or both) datasets have a known stratification, here simulated as
 S <- rbinom(500,1,0.5)
 ## specify this in pcs.model as
 m1 <- pcs.model(pcs, group=1, Y=Y1, stratum=S)
#> selecting 8 components out of 10 to capture 0.8299584 of total variance.

 ## test colocalisation using PCs
 coloc.test(m1,m2,plot.coeff=FALSE,bayes=FALSE)
#> Warning: 1 variables dropped from regression X:
#> 	 S 
#>           eta.hat         chisquare                 n p.value.chisquare 
#>         2.8941037         8.6521699         8.0000000         0.2786026