An ExpressionSet Object Storing Simulated Genotype Data

Description

An ExpressionSet object storing simulated genotype data. The minor allele frequency (MAF) of cases has the same prior as that of controls.

Usage

data("esSim")

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\}

We set the number of cases as 100, the number of controls as 100, and the number of SNPs as 1000.

The hyperparameters are \alpha_{+}=2, \beta_{+}=5, \pi_{+}=0.1, \alpha_{0}=2, \beta_{0}=5, \pi_{0}=0.8, \alpha_{-}=2, \beta_{-}=5, \pi_{-}=0.1.

Note that when we generate MAFs from the half-flat shape bivariate priors, we might get very small MAFs or get MAFs >0.5. In these cased, we then delete this SNP.

So the final number of SNPs generated might be less than the initially-set number 1000 of SNPs.

For the dataset stored in esSim, there are 872 SNPs. 83 SNPs are in cluster -, 714 SNPs are in cluster 0, and 75 SNPs are in cluster +.

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(esSim)
print(esSim)

pDat=pData(esSim)
print(pDat[1:2,])
print(table(pDat$memSubjs))

fDat=fData(esSim)
print(fDat[1:2,])
print(table(fDat$memGenes))
print(table(fDat$memGenes2))

An ExpressionSet Object Storing Simulated Genotype Data

Description

An ExpressionSet object storing simulated genotype data. The minor allele frequency (MAF) of cases has different prior than that of controls.

Usage

data("esSimDiffPriors")

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_{x+}\left(\theta_{x+}\right)h_{y+}\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_{x+} is the probability density function of the beta distribution Beta\left(\alpha_{x+}, \beta_{x+}\right). The function h_{y+} is the probability density function of the beta distribution Beta\left(\alpha_{y+}, \beta_{y+}\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_{x-}\left(\theta_{x-}\right)h_{y-}\left(\theta_{y-}\right) I\left(\theta_{x-}>\theta_{y-}\right).

The function h_{x-} is the probability density function of the beta distribution Beta\left(\alpha_{x-}, \beta_{x-}\right). The function h_{y-} is the probability density function of the beta distribution Beta\left(\alpha_{y-}, \beta_{y-}\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\}

We set the number of cases as 100, the number of controls as 100, and the number of SNPs as 1000.

The hyperparameters are \alpha_{x+}=2, \beta_{x+}=3, \alpha_{y+}=2, \beta_{y+}=8, \pi_{+}=0.1,

\alpha_{0}=2, \beta_{0}=5, \pi_{0}=0.8,

\alpha_{x-}=2, \beta_{x-}=8, \alpha_{y-}=2, \beta_{y-}=3, \pi_{-}=0.1.

Note that when we generate MAFs from the half-flat shape bivariate priors, we might get very small MAFs or get MAFs >0.5. In these cased, we then delete this SNP.

So the final number of SNPs generated might be less than the initially-set number 1000 of SNPs.

For the dataset stored in esSim, there are 838 SNPs. 64 SNPs are in cluster -, 708 SNPs are in cluster 0, and 66 SNPs are in cluster +.

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)

pDat=pData(esSimDiffPriors)
print(pDat[1:2,])
print(table(pDat$memSubjs))

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

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

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

Author(s)

Yan Xu <yanxu@uvic.ca>, Li Xing <sfulxing@gmail.com>, Jessica Su <rejas@channing.harvard.edu>, Xuekui Zhang <xuekui@uvic.ca>, Weiliang Qiu <Weiliang.Qiu@gmail.com>

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))

Estimate SNP cluster membership

Description

Estimate SNP cluster membership. Only update cluster mixture proportions. Assume all 3 clusters have the same set of hyperparameters.

Usage

estMemSNPs.oneSetHyperPara(es, 
           var.memSubjs = "memSubjs", 
           eps = 1.0e-3,
           MaxIter = 50, 
           bVec = rep(3, 3), 
           pvalAdjMethod = "none", 
           method = "FDR",
           fdr = 0.05,
           verbose = FALSE)

Arguments

Details

We characterize the distribution of genotypes of SNPs by a mixture of 3 Bayesian hierarchical models. The 3 Bayeisan hierarchical models correspond to 3 clusters of SNPs.

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, \beta).

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

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 10 elements

Author(s)

Yan Xu <yanxu@uvic.ca>, Li Xing <sfulxing@gmail.com>, Jessica Su <rejas@channing.harvard.edu>, Xuekui Zhang <xuekui@uvic.ca>, Weiliang Qiu <Weiliang.Qiu@gmail.com>

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)
fDat = fData(esSimDiffPriors)
print(fDat[1:2,])
print(table(fDat$memGenes))

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

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

Simulate Genotype Data from a Mixture of 3 Bayesian Hierarchical Models

Description

Simulate Genotype Data from a Mixture of 3 Bayesian Hierarchical Models. The minor allele frequency (MAF) of cases has the same prior as that of controls.

Usage

simGenoFunc(nCases = 100, 
            nControls = 100, 
            nSNPs = 1000, 
            alpha.p = 2, 
            beta.p = 5, 
            pi.p = 0.1, 
            alpha0 = 2, 
            beta0 = 5, 
            pi0 = 0.8, 
            alpha.n = 2, 
            beta.n = 5, 
            pi.n = 0.1, 
            low = 0.02, 
            upp = 0.5, 
            verbose = FALSE)

Arguments

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\}

Note that when we generate MAFs from the half-flat shape bivariate priors, we might get very small MAFs or get MAFs >0.5. In these cased, we then delete this SNP.

So the final number of SNPs generated might be less than the initially-set number of SNPs.

Value

An ExpressionSet object stores genotype data.

Author(s)

Yan Xu <yanxu@uvic.ca>, Li Xing <sfulxing@gmail.com>, Jessica Su <rejas@channing.harvard.edu>, Xuekui Zhang <xuekui@uvic.ca>, Weiliang Qiu <Weiliang.Qiu@gmail.com>

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

set.seed(2) 

esSim = simGenoFunc(
  nCases = 100,
  nControls = 100,
  nSNPs = 500,
  alpha.p = 2, beta.p = 5, pi.p = 0.1,
  alpha0 = 2, beta0 = 5, pi0 = 0.8,
  alpha.n = 2, beta.n = 5, pi.n = 0.1,
  low = 0.02, upp = 0.5, verbose = FALSE
)

print(esSim)
pDat = pData(esSim)
print(pDat[1:2,])
print(table(pDat$memSubjs))

fDat = fData(esSim)
print(fDat[1:2,])
print(table(fDat$memGenes))
print(table(fDat$memGenes2))

Simulate Genotype Data from a Mixture of 3 Bayesian Hierarchical Models

Description

Simulate Genotype Data from a Mixture of 3 Bayesian Hierarchical Models. The minor allele frequency (MAF) of cases has different priors than that of controls.

Usage

simGenoFuncDiffPriors(
    nCases = 100, 
    nControls = 100, 
    nSNPs = 1000, 
    alpha.p.ca = 2, 
    beta.p.ca = 3, 
    alpha.p.co = 2, 
    beta.p.co = 8, 
    pi.p = 0.1, 
    alpha0 = 2, 
    beta0 = 5, 
    pi0 = 0.8, 
    alpha.n.ca = 2, 
    beta.n.ca = 8, 
    alpha.n.co = 2, 
    beta.n.co = 3, 
    pi.n = 0.1, 
    low = 0.02, 
    upp = 0.5, 
    verbose = FALSE)

Arguments

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\}

Note that when we generate MAFs from the half-flat shape bivariate priors, we might get very small MAFs or get MAFs >0.5. In these cased, we then delete this SNP.

So the final number of SNPs generated might be less than the initially-set number of SNPs.

Value

An ExpressionSet object stores genotype data.

Author(s)

Yan Xu <yanxu@uvic.ca>, Li Xing <sfulxing@gmail.com>, Jessica Su <rejas@channing.harvard.edu>, Xuekui Zhang <xuekui@uvic.ca>, Weiliang Qiu <Weiliang.Qiu@gmail.com>

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

set.seed(2)

esSimDiffPriors = simGenoFuncDiffPriors(
  nCases = 100,
  nControls = 100,
  nSNPs = 500,
  alpha.p.ca = 2, beta.p.ca = 3,
  alpha.p.co = 2, beta.p.co = 8, pi.p = 0.1,
  alpha0 = 2, beta0 = 5, pi0 = 0.8,
  alpha.n.ca = 2, beta.n.ca = 8,
  alpha.n.co = 2, beta.n.co = 3, pi.n = 0.1,
  low = 0.02, upp = 0.5, verbose = FALSE
)

print(esSimDiffPriors)

pDat = pData(esSimDiffPriors)
print(pDat[1:2,])
print(table(pDat$memSubjs))

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