estMemSNPs {GWASbyCluster}

Estimate SNP cluster membership

Description

Estimate SNP cluster membership. Only update cluster mixture proportions. Assume the 3 clusters have different sets of hyperparameters.

Usage

estMemSNPs(es, 
           var.memSubjs = "memSubjs", 
           eps = 0.001, 
           MaxIter = 50, 
           bVec = rep(3, 3), 
           pvalAdjMethod = "fdr", 
           method = "FDR",
           fdr = 0.05,
           verbose = FALSE)

Arguments

es

An ExpressionSet object storing SNP genotype data. It contains 3 matrices. The first matrix, which can be extracted by exprs method (e.g., exprs(es)), stores genotype data, with rows are SNPs and columns are subjects.

The second matrix, which can be extracted by pData method (e.g., pData(es)), stores phenotype data describing subjects. Rows are subjects, and columns are phenotype variables.

The third matrix, which can be extracted by fData method (e.g., fData(es)), stores feature data describing SNPs. Rows are SNPs and columns are feature variables.

var.memSubjs

character. The name of the phenotype variable indicating subject's case-control status. It must take only two values: 1 indicating case and 0 indicating control.

eps

numeric. A small positive number as threshold for convergence of EM algorithm.

MaxIter

integer. A positive integer indicating maximum iteration in EM algorithm.

bVec

numeric. A vector of 2 elements. Indicates the parameters of the symmetric Dirichlet prior for proportion mixtures.

pvalAdjMethod

character. Indicating p-value adjustment method. c.f. p.adjust.

method

method to obtain SNP cluster membership based on the responsibility matrix. The default value is “FDR”. The other possible value is “max”. see details.

fdr

numeric. A small positive FDR threshold used to call SNP cluster membership

verbose

logical. Indicating if intermediate and final results should be output.

Details

In this simulation, we generate additive-coded genotypes for 3 clusters of SNPs based on a mixture of 3 Bayesian hierarchical models.

In cluster +, the minor allele frequency (MAF) \theta_{x+} of cases is greater than the MAF \theta_{y+} of controls.

In cluster 0, the MAF \theta_{0} of cases is equal to the MAF of controls.

In cluster -, the MAF \theta_{x-} of cases is smaller than the MAF \theta_{y-} of controls.

The proportions of the 3 clusters of SNPs are \pi_{+}, \pi_{0}, and \pi_{-}, respectively.

We assume a “half-flat shape” bivariate prior for the MAF in cluster +

2h_{+}\left(\theta_{x+}\right)h_{+}\left(\theta_{y+}\right) I\left(\theta_{x+}>\theta_{y+}\right),

where I(a) is hte indicator function taking value 1 if the event a is true, and value 0 otherwise. The function h_{+} is the probability density function of the beta distribution Beta\left(\alpha_{+}, \beta_{+}\right).

We assume \theta_{0} has the beta prior Beta(\alpha_0, \beta_0).

We also assume a “half-flat shape” bivariate prior for the MAF in cluster -

2h_{-}\left(\theta_{x-}\right)h_{-}\left(\theta_{y-}\right) I\left(\theta_{x-}>\theta_{y-}\right).

The function h_{-} is the probability density function of the beta distribution Beta\left(\alpha_{-}, \beta_{-}\right).

Given a SNP, we assume Hardy-Weinberg equilibrium holds for its genotypes. That is, given MAF \theta, the probabilities of genotypes are

Pr(geno=2) = \theta^2

Pr(geno=1) = 2\theta\left(1-\theta\right)

Pr(geno=0) = \left(1-\theta\right)^2

We also assume the genotypes 0 (wild-type), 1 (heterozygote), and 2 (mutation) follows a multinomial distribution Multinomial\left\{1, \left[ \theta^2, 2\theta\left(1-\theta\right), \left(1-\theta\right)^2 \right]\right\}

For each SNP, we calculat its posterior probabilities that it belongs to cluster k. This forms a matrix with 3 columns. Rows are SNPs. The 1st column is the posterior probability that the SNP belongs to cluster -. The 2nd column is the posterior probability that the SNP belongs to cluster 0. The 3rd column is the posterior probability that the SNP belongs to cluster +. We call this posterior probability matrix as responsibility matrix. To determine which cluster a SNP eventually belongs to, we can use 2 methods. The first method (the default method) is “FDR” method, which will use FDR criterion to determine SNP cluster membership. The 2nd method is use the maximum posterior probability to decide which cluster a SNP belongs to.

Value

A list of 12 elements

wMat

matrix of posterior probabilities. The rows are SNPs. There are 3 columns. The first column is the posterior probability that a SNP belongs to cluster - given genotypes of subjects. The second column is the posterior probability that a SNP belongs to cluster 0 given genotypes of subjects. The third column is the posterior probability that a SNP belongs to cluster + given genotypes of subjects.

memSNPs

a vector of SNP cluster membership for the 3-cluster partitionfrom the mixture of 3 Bayesian hierarchical models.

memSNPs2

a vector of binary SNP cluster membership. 1 indicates the SNP has different MAFs between cases and controls. 0 indicates the SNP has the same MAF in cases as that in controls.

piVec

a vector of cluster mixture proportions.

alpha.p

the first shape parameter of the beta prior for MAF obtaind from initial 3-cluster partitions based on GWAS for cluster +.

beta.p

the second shape parameter of the beta prior for MAF obtaind from initial 3-cluster partitions based on GWAS for cluster +.

alpha0

the first shape parameter of the beta prior for MAF obtaind from initial 3-cluster partitions based on GWAS for cluster 0.

beta0

the second shape parameter of the beta prior for MAF obtaind from initial 3-cluster partitions based on GWAS for cluster 0.

alpha.n

the first shape parameter of the beta prior for MAF obtaind from initial 3-cluster partitions based on GWAS for cluster -.

beta.n

the second shape parameter of the beta prior for MAF obtaind from initial 3-cluster partitions based on GWAS for cluster -.

loop

number of iteration in EM algorithm

diff

sum of the squared difference of cluster mixture proportions between current iteration and previous iteration in EM algorithm. if eps < eps, we claim the EM algorithm converges.

res.limma

object returned by limma

Author(s)

Yan Xu <[e-mail]>, Li Xing <[e-mail]>, Jessica Su <[e-mail]>, Xuekui Zhang <[e-mail]>, Weiliang Qiu <[e-mail]>

References

Yan X, Xing L, Su J, Zhang X, Qiu W. Model-based clustering for identifying disease-associated SNPs in case-control genome-wide association studies. Scientific Reports 9, Article number: 13686 (2019) https://www.nature.com/articles/s41598-019-50229-6.

Examples

data(esSimDiffPriors)
print(esSimDiffPriors)

es=esSimDiffPriors[1:500,]
fDat = fData(es)
print(fDat[1:2,])
print(table(fDat$memGenes))

res = estMemSNPs(
  es = es, 
  var.memSubjs = "memSubjs")

print(table(fDat$memGenes, res$memSNPs))