Despite the recent advances of modern GWAS methods, it is still remains an important problem of addressing calculation an effect size and corresponding p-value for the whole gene rather than for single variant. We developed an R-package rqt, which offers gene-level GWAS meta-analysis. The package can be easily included into bioinformatics pipeline or used as stand-alone. Contact: ilya.zhbannikov@duke.edu for questions of usage the or any other issues.
Below we provide several examples that show GWAS meta-analysis on gene-level layer.
The workflow of gene-level meta analysis consists of the following steps: (i) reducing the number of predictors, thereby alleviating correlation problem in variants (accounting for LD); (ii) then the regression mod-el is fitted on the reduced dataset to obtain corresponding regression coefficient (“effect sizes”); (iii) these coefficients are then to be pooled into a total index representing a total gene-level effect size and corresponding statistics is calculated. P- and q- values are then calculated using this statistics from asymptotic approximation or permutation procedure; (iv) the final step is combining gene-level p-values calculated from each study with Fisher’s combined probability method.
In order to install the package, the user must first install R (). After that, can be installed with:
devtools::install_github("izhbannikov/rqt", buildVignette=TRUE)In requires the following datasets: (i) (a by 1) matrix (i.e. a vector); and (ii) - an object of class containing one assay: (a by ) matrix, where - is the total number of individuals in the study and is the total number of genetic variants. Optionally, can accept covariates, in form of by matrix, where is the total number of covariates used in the study. Phenotype can be dichotomous (0/1, where 1 indicates control and 0 case).
In meta-analysis, requires a list of ( - number of datasets used in meta-analysis) and optionally it accepts covariates in form described above.
library(rqt)
data <- data.matrix(read.table(system.file("extdata/test.bin1.dat",
package="rqt"), header=TRUE))
pheno <- data[,1]
geno <- data[, 2:dim(data)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj <- rqt(phenotype=pheno, genotype=geno.obj)
res <- geneTest(obj, method="pca", out.type = "D")
print(res)## Phenotype:
## [1] 1 1 1 1 1 1
## ...
##
## Genotype:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## [1,] 0 0 0 1 0 2 1 0 0 0 2 1 2 0 0 2 1 1 0 0 1 0 0 0 0 2 1 0
## [2,] 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 2 1 2 1 0 1 0 0 0 2 1 1
## [3,] 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0
## [4,] 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 0 1 0 2
## [5,] 0 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 2 0
## [6,] 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 2 0 0 1 0 1 1 0 0 0 2 0 0
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [1,] 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1
## [2,] 2 0 0 0 1 2 1 2 0 1 1 0 1 0 0 1 0 0 0 2 0 1 0 1 0
## [3,] 0 0 0 0 0 1 1 1 2 2 0 0 0 1 0 1 0 2 1 1 1 0 0 0 1
## [4,] 0 0 0 1 2 2 0 1 1 1 1 0 0 0 1 1 0 0 1 1 2 0 1 0 2
## [5,] 0 0 0 0 1 1 0 0 1 1 0 0 2 2 0 0 0 1 2 1 0 0 0 0 1
## [6,] 0 0 0 0 1 0 1 1 0 0 2 0 1 1 2 0 0 0 1 1 1 0 0 1 1
## 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## [1,] 1 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 2 1 1 1 0 0 2 1
## [2,] 2 1 0 1 0 0 0 1 0 1 2 2 1 1 0 1 0 1 0 0 0 0 1 1 2
## [3,] 2 1 0 0 0 0 0 0 0 1 0 0 1 1 0 1 2 1 0 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0 2 1 0 1 0 0 2 1 0 0 1 1 1 0 0 0 0 0 1
## [5,] 1 1 1 0 0 0 2 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 0 1
## [6,] 0 1 1 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 1 1 1 1 1 1
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## [1,] 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
## [2,] 1 1 2 0 0 0 1 2 0 2 1 0 0 2 1 0 0 0 0 1 1 2
## [3,] 1 1 0 0 0 0 2 0 0 2 0 0 0 0 1 1 0 0 1 2 0 0
## [4,] 0 1 2 0 1 0 1 2 0 0 0 0 0 0 1 0 0 0 1 0 0 0
## [5,] 0 2 0 1 0 0 2 0 0 0 0 1 0 0 1 0 0 1 1 2 1 1
## [6,] 2 1 0 0 1 0 0 0 0 2 0 1 0 0 0 0 0 1 1 0 1 1
## ...
##
## Covariates:
## data frame with 0 columns and 0 rows
##
##
## Results:
##
## $Qstatistic
## Q1 Q2 Q3
## 1 0.3798991 13.83897 0.3798991
##
## $pValue
## pVal1 pVal2 pVal3
## 1 0.5376573 0.5862745 0.7862228
##
## $beta
## [1] -0.005612555
##
## $var.pooled
## [1] 8.291881e-05
##
## $mean.vif
## [1] 1
##
## $model
## [1] "PCA"library(rqt)
data <- data.matrix(read.table(system.file("extdata/test.cont1.dat",
package="rqt"), header = TRUE))
pheno <- data[,1]
geno <- data[, 2:dim(data)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj <- rqt(phenotype=pheno, genotype=geno.obj)
res <- geneTest(obj, method="pca", out.type = "C")
print(res)## Phenotype:
## [1] 3.422452 2.457394 2.708564 4.589394 5.461723 4.438707
## ...
##
## Genotype:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## [1,] 0 0 0 1 0 2 1 0 0 0 2 1 2 0 0 2 1 1 0 0 1 0 0 0 0 2 1 0
## [2,] 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 2 1 2 1 0 1 0 0 0 2 1 1
## [3,] 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0
## [4,] 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 0 1 0 2
## [5,] 0 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 2 0
## [6,] 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 2 0 0 1 0 1 1 0 0 0 2 0 0
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [1,] 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1
## [2,] 2 0 0 0 1 2 1 2 0 1 1 0 1 0 0 1 0 0 0 2 0 1 0 1 0
## [3,] 0 0 0 0 0 1 1 1 2 2 0 0 0 1 0 1 0 2 1 1 1 0 0 0 1
## [4,] 0 0 0 1 2 2 0 1 1 1 1 0 0 0 1 1 0 0 1 1 2 0 1 0 2
## [5,] 0 0 0 0 1 1 0 0 1 1 0 0 2 2 0 0 0 1 2 1 0 0 0 0 1
## [6,] 0 0 0 0 1 0 1 1 0 0 2 0 1 1 2 0 0 0 1 1 1 0 0 1 1
## 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## [1,] 1 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 2 1 1 1 0 0 2 1
## [2,] 2 1 0 1 0 0 0 1 0 1 2 2 1 1 0 1 0 1 0 0 0 0 1 1 2
## [3,] 2 1 0 0 0 0 0 0 0 1 0 0 1 1 0 1 2 1 0 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0 2 1 0 1 0 0 2 1 0 0 1 1 1 0 0 0 0 0 1
## [5,] 1 1 1 0 0 0 2 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 0 1
## [6,] 0 1 1 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 1 1 1 1 1 1
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## [1,] 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
## [2,] 1 1 2 0 0 0 1 2 0 2 1 0 0 2 1 0 0 0 0 1 1 2
## [3,] 1 1 0 0 0 0 2 0 0 2 0 0 0 0 1 1 0 0 1 2 0 0
## [4,] 0 1 2 0 1 0 1 2 0 0 0 0 0 0 1 0 0 0 1 0 0 0
## [5,] 0 2 0 1 0 0 2 0 0 0 0 1 0 0 1 0 0 1 1 2 1 1
## [6,] 2 1 0 0 1 0 0 0 0 2 0 1 0 0 0 0 0 1 1 0 1 1
## ...
##
## Covariates:
## data frame with 0 columns and 0 rows
##
##
## Results:
##
## $Qstatistic
## Q1 Q2 Q3
## 1 0.7901911 46.38915 12.15855
##
## $pValue
## pVal1 pVal2 pVal3
## 1 0.3740423 6.908983e-05 0.001169993
##
## $beta
## [1] -0.0164608
##
## $var.pooled
## [1] 0.0003429017
##
## $mean.vif
## [1] 1
##
## $model
## [1] "PCA"This method is used for continous outcome, i.e. out.type = “C”.
library(rqt)
data <- data.matrix(read.table(system.file("extdata/test.cont1.dat",
package="rqt"), header = TRUE))
pheno <- data[,1]
geno <- data[, 2:dim(data)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj <- rqt(phenotype=pheno, genotype=geno.obj)
res <- geneTest(obj, method="pls", out.type = "C")
print(res)## Phenotype:
## [1] 3.422452 2.457394 2.708564 4.589394 5.461723 4.438707
## ...
##
## Genotype:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## [1,] 0 0 0 1 0 2 1 0 0 0 2 1 2 0 0 2 1 1 0 0 1 0 0 0 0 2 1 0
## [2,] 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 2 1 2 1 0 1 0 0 0 2 1 1
## [3,] 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0
## [4,] 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 0 1 0 2
## [5,] 0 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 2 0
## [6,] 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 2 0 0 1 0 1 1 0 0 0 2 0 0
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [1,] 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1
## [2,] 2 0 0 0 1 2 1 2 0 1 1 0 1 0 0 1 0 0 0 2 0 1 0 1 0
## [3,] 0 0 0 0 0 1 1 1 2 2 0 0 0 1 0 1 0 2 1 1 1 0 0 0 1
## [4,] 0 0 0 1 2 2 0 1 1 1 1 0 0 0 1 1 0 0 1 1 2 0 1 0 2
## [5,] 0 0 0 0 1 1 0 0 1 1 0 0 2 2 0 0 0 1 2 1 0 0 0 0 1
## [6,] 0 0 0 0 1 0 1 1 0 0 2 0 1 1 2 0 0 0 1 1 1 0 0 1 1
## 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## [1,] 1 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 2 1 1 1 0 0 2 1
## [2,] 2 1 0 1 0 0 0 1 0 1 2 2 1 1 0 1 0 1 0 0 0 0 1 1 2
## [3,] 2 1 0 0 0 0 0 0 0 1 0 0 1 1 0 1 2 1 0 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0 2 1 0 1 0 0 2 1 0 0 1 1 1 0 0 0 0 0 1
## [5,] 1 1 1 0 0 0 2 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 0 1
## [6,] 0 1 1 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 1 1 1 1 1 1
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## [1,] 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
## [2,] 1 1 2 0 0 0 1 2 0 2 1 0 0 2 1 0 0 0 0 1 1 2
## [3,] 1 1 0 0 0 0 2 0 0 2 0 0 0 0 1 1 0 0 1 2 0 0
## [4,] 0 1 2 0 1 0 1 2 0 0 0 0 0 0 1 0 0 0 1 0 0 0
## [5,] 0 2 0 1 0 0 2 0 0 0 0 1 0 0 1 0 0 1 1 2 1 1
## [6,] 2 1 0 0 1 0 0 0 0 2 0 1 0 0 0 0 0 1 1 0 1 1
## ...
##
## Covariates:
## data frame with 0 columns and 0 rows
##
##
## Results:
##
## $Qstatistic
## Q1 Q2 Q3
## 1 0.1590242 0.0001860569 0.1590242
##
## $pValue
## pVal1 pVal2 pVal3
## 1 0.6900565 1 0.9035676
##
## $beta
## [1] 0.4623658
##
## $var.pooled
## [1] 1.344337
##
## $mean.vif
## [1] 1
##
## $model
## Partial least squares regression , fitted with the kernel algorithm.
## Call:
## plsr(formula = pheno ~ ., ncomp = numcomp, data = dd, validation = "none")This method of data preprocessing is used for dichotomous outcome.
library(rqt)
data <- data.matrix(read.table(system.file("extdata/test.bin1.dat",
package="rqt"), header=TRUE))
pheno <- data[,1]
geno <- data[, 2:dim(data)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj <- rqt(phenotype=pheno, genotype=geno.obj)
# Not yet supported, sorry!
#res <- geneTest(obj, method="pls", out.type = "D", scale = TRUE)
print(res)## Phenotype:
## [1] 3.422452 2.457394 2.708564 4.589394 5.461723 4.438707
## ...
##
## Genotype:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## [1,] 0 0 0 1 0 2 1 0 0 0 2 1 2 0 0 2 1 1 0 0 1 0 0 0 0 2 1 0
## [2,] 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 2 1 2 1 0 1 0 0 0 2 1 1
## [3,] 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0
## [4,] 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 0 1 0 2
## [5,] 0 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 2 0
## [6,] 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 2 0 0 1 0 1 1 0 0 0 2 0 0
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [1,] 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1
## [2,] 2 0 0 0 1 2 1 2 0 1 1 0 1 0 0 1 0 0 0 2 0 1 0 1 0
## [3,] 0 0 0 0 0 1 1 1 2 2 0 0 0 1 0 1 0 2 1 1 1 0 0 0 1
## [4,] 0 0 0 1 2 2 0 1 1 1 1 0 0 0 1 1 0 0 1 1 2 0 1 0 2
## [5,] 0 0 0 0 1 1 0 0 1 1 0 0 2 2 0 0 0 1 2 1 0 0 0 0 1
## [6,] 0 0 0 0 1 0 1 1 0 0 2 0 1 1 2 0 0 0 1 1 1 0 0 1 1
## 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## [1,] 1 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 2 1 1 1 0 0 2 1
## [2,] 2 1 0 1 0 0 0 1 0 1 2 2 1 1 0 1 0 1 0 0 0 0 1 1 2
## [3,] 2 1 0 0 0 0 0 0 0 1 0 0 1 1 0 1 2 1 0 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0 2 1 0 1 0 0 2 1 0 0 1 1 1 0 0 0 0 0 1
## [5,] 1 1 1 0 0 0 2 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 0 1
## [6,] 0 1 1 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 1 1 1 1 1 1
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## [1,] 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
## [2,] 1 1 2 0 0 0 1 2 0 2 1 0 0 2 1 0 0 0 0 1 1 2
## [3,] 1 1 0 0 0 0 2 0 0 2 0 0 0 0 1 1 0 0 1 2 0 0
## [4,] 0 1 2 0 1 0 1 2 0 0 0 0 0 0 1 0 0 0 1 0 0 0
## [5,] 0 2 0 1 0 0 2 0 0 0 0 1 0 0 1 0 0 1 1 2 1 1
## [6,] 2 1 0 0 1 0 0 0 0 2 0 1 0 0 0 0 0 1 1 0 1 1
## ...
##
## Covariates:
## data frame with 0 columns and 0 rows
##
##
## Results:
##
## $Qstatistic
## Q1 Q2 Q3
## 1 0.1590242 0.0001860569 0.1590242
##
## $pValue
## pVal1 pVal2 pVal3
## 1 0.6900565 1 0.9035676
##
## $beta
## [1] 0.4623658
##
## $var.pooled
## [1] 1.344337
##
## $mean.vif
## [1] 1
##
## $model
## Partial least squares regression , fitted with the kernel algorithm.
## Call:
## plsr(formula = pheno ~ ., ncomp = numcomp, data = dd, validation = "none")Quite often, researchers want to supply not only genetic data but also specific covariates, representic some physiological parameters or environment (for example, to evaluate hyphoteses of gene-environment interactions). In such cases, the package can accept additional covariates, in form of by matrix, as provided below:
library(rqt)
data <- data.matrix(read.table(system.file("extdata/test.bin1.dat",
package="rqt"), header = TRUE))
pheno <- data[,1]
geno <- data[, 2:dim(data)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
covars <- read.table(system.file("extdata/test.cova1.dat",package="rqt"),
header=TRUE)
obj <- rqt(phenotype=pheno, genotype=geno.obj, covariates = covars)
res <- geneTest(obj, method="pca", out.type = "D")
print(res)## Phenotype:
## [1] 1 1 1 1 1 1
## ...
##
## Genotype:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## [1,] 0 0 0 1 0 2 1 0 0 0 2 1 2 0 0 2 1 1 0 0 1 0 0 0 0 2 1 0
## [2,] 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 2 1 2 1 0 1 0 0 0 2 1 1
## [3,] 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0
## [4,] 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 0 1 0 2
## [5,] 0 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 2 0
## [6,] 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 2 0 0 1 0 1 1 0 0 0 2 0 0
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [1,] 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1
## [2,] 2 0 0 0 1 2 1 2 0 1 1 0 1 0 0 1 0 0 0 2 0 1 0 1 0
## [3,] 0 0 0 0 0 1 1 1 2 2 0 0 0 1 0 1 0 2 1 1 1 0 0 0 1
## [4,] 0 0 0 1 2 2 0 1 1 1 1 0 0 0 1 1 0 0 1 1 2 0 1 0 2
## [5,] 0 0 0 0 1 1 0 0 1 1 0 0 2 2 0 0 0 1 2 1 0 0 0 0 1
## [6,] 0 0 0 0 1 0 1 1 0 0 2 0 1 1 2 0 0 0 1 1 1 0 0 1 1
## 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## [1,] 1 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 2 1 1 1 0 0 2 1
## [2,] 2 1 0 1 0 0 0 1 0 1 2 2 1 1 0 1 0 1 0 0 0 0 1 1 2
## [3,] 2 1 0 0 0 0 0 0 0 1 0 0 1 1 0 1 2 1 0 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0 2 1 0 1 0 0 2 1 0 0 1 1 1 0 0 0 0 0 1
## [5,] 1 1 1 0 0 0 2 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 0 1
## [6,] 0 1 1 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 1 1 1 1 1 1
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## [1,] 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
## [2,] 1 1 2 0 0 0 1 2 0 2 1 0 0 2 1 0 0 0 0 1 1 2
## [3,] 1 1 0 0 0 0 2 0 0 2 0 0 0 0 1 1 0 0 1 2 0 0
## [4,] 0 1 2 0 1 0 1 2 0 0 0 0 0 0 1 0 0 0 1 0 0 0
## [5,] 0 2 0 1 0 0 2 0 0 0 0 1 0 0 1 0 0 1 1 2 1 1
## [6,] 2 1 0 0 1 0 0 0 0 2 0 1 0 0 0 0 0 1 1 0 1 1
## ...
##
## Covariates:
## COV1
## 1 -0.612463927
## 2 -0.464158885
## 3 0.006153597
## 4 -0.732109468
## 5 -0.223530136
## 6 -0.744903822
##
##
## Results:
##
## $Qstatistic
## Q1 Q2 Q3
## 1 0.3732768 13.07139 0.3732768
##
## $pValue
## pVal1 pVal2 pVal3
## 1 0.5412235 0.6440566 0.7896537
##
## $beta
## [1] -0.005553824
##
## $var.pooled
## [1] 8.263294e-05
##
## $mean.vif
## [1] 1
##
## $model
## [1] "PCA"For continous phenotype:
library(rqt)
data <- data.matrix(read.table(system.file("extdata/test.cont1.dat",
package="rqt"), header = TRUE))
pheno <- data[,1]
geno <- data[, 2:dim(data)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
covars <- read.table(system.file("extdata/test.cova1.dat",package="rqt"),
header=TRUE)
obj <- rqt(phenotype=pheno, genotype=geno.obj, covariates = covars)
res <- geneTest(obj, method="pca", out.type = "C")
print(res)## Phenotype:
## [1] 3.422452 2.457394 2.708564 4.589394 5.461723 4.438707
## ...
##
## Genotype:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
## [1,] 0 0 0 1 0 2 1 0 0 0 2 1 2 0 0 2 1 1 0 0 1 0 0 0 0 2 1 0
## [2,] 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 2 1 2 1 0 1 0 0 0 2 1 1
## [3,] 0 0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0
## [4,] 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 0 1 0 2
## [5,] 0 0 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 2 0
## [6,] 0 0 1 1 1 0 0 1 0 1 1 0 1 0 0 2 0 0 1 0 1 1 0 0 0 2 0 0
## 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [1,] 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 1
## [2,] 2 0 0 0 1 2 1 2 0 1 1 0 1 0 0 1 0 0 0 2 0 1 0 1 0
## [3,] 0 0 0 0 0 1 1 1 2 2 0 0 0 1 0 1 0 2 1 1 1 0 0 0 1
## [4,] 0 0 0 1 2 2 0 1 1 1 1 0 0 0 1 1 0 0 1 1 2 0 1 0 2
## [5,] 0 0 0 0 1 1 0 0 1 1 0 0 2 2 0 0 0 1 2 1 0 0 0 0 1
## [6,] 0 0 0 0 1 0 1 1 0 0 2 0 1 1 2 0 0 0 1 1 1 0 0 1 1
## 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
## [1,] 1 1 1 1 0 0 0 0 1 1 0 1 1 1 0 1 1 2 1 1 1 0 0 2 1
## [2,] 2 1 0 1 0 0 0 1 0 1 2 2 1 1 0 1 0 1 0 0 0 0 1 1 2
## [3,] 2 1 0 0 0 0 0 0 0 1 0 0 1 1 0 1 2 1 0 1 1 0 0 0 0
## [4,] 1 1 0 0 0 0 2 1 0 1 0 0 2 1 0 0 1 1 1 0 0 0 0 0 1
## [5,] 1 1 1 0 0 0 2 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 0 1
## [6,] 0 1 1 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 1 1 1 1 1 1
## 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## [1,] 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
## [2,] 1 1 2 0 0 0 1 2 0 2 1 0 0 2 1 0 0 0 0 1 1 2
## [3,] 1 1 0 0 0 0 2 0 0 2 0 0 0 0 1 1 0 0 1 2 0 0
## [4,] 0 1 2 0 1 0 1 2 0 0 0 0 0 0 1 0 0 0 1 0 0 0
## [5,] 0 2 0 1 0 0 2 0 0 0 0 1 0 0 1 0 0 1 1 2 1 1
## [6,] 2 1 0 0 1 0 0 0 0 2 0 1 0 0 0 0 0 1 1 0 1 1
## ...
##
## Covariates:
## COV1
## 1 -0.612463927
## 2 -0.464158885
## 3 0.006153597
## 4 -0.732109468
## 5 -0.223530136
## 6 -0.744903822
##
##
## Results:
##
## $Qstatistic
## Q1 Q2 Q3
## 1 0.7761103 44.94641 11.63278
##
## $pValue
## pVal1 pVal2 pVal3
## 1 0.3783334 0.000115889 0.001549197
##
## $beta
## [1] -0.01629069
##
## $var.pooled
## [1] 0.0003419444
##
## $mean.vif
## [1] 1
##
## $model
## [1] "PCA"library(rqt)
data1 <- data.matrix(read.table(system.file("extdata/phengen2.dat",
package="rqt"), skip=1))
pheno <- data1[,1]
geno <- data1[, 2:dim(data1)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj1 <- rqt(phenotype=pheno, genotype=geno.obj)
data2 <- data.matrix(read.table(system.file("extdata/phengen3.dat",
package="rqt"), skip=1))
pheno <- data2[,1]
geno <- data2[, 2:dim(data2)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj2 <- rqt(phenotype=pheno, genotype=geno.obj)
data3 <- data.matrix(read.table(system.file("extdata/phengen.dat",
package="rqt"), skip=1))
pheno <- data3[,1]
geno <- data3[, 2:dim(data3)[2]]
colnames(geno) <- paste(seq(1, dim(geno)[2]))
geno.obj <- SummarizedExperiment(geno)
obj3 <- rqt(phenotype=pheno, genotype=geno.obj)
res.meta <- geneTestMeta(list(obj1, obj2, obj3))
print(res.meta)## $final.pvalue
## [1] 0.005502565
##
## $pvalueList
## [1] 0.001837563 0.225098063 0.179679688sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] rqt_1.16.0 SummarizedExperiment_1.20.0
## [3] Biobase_2.50.0 GenomicRanges_1.42.0
## [5] GenomeInfoDb_1.26.0 IRanges_2.24.0
## [7] S4Vectors_0.28.0 BiocGenerics_0.36.0
## [9] MatrixGenerics_1.2.0 matrixStats_0.57.0
## [11] BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.0.3 foreach_1.5.1 carData_3.0-4
## [4] tmvnsim_1.0-2 sn_1.6-2 Rdpack_2.0
## [7] BiocManager_1.30.10 TFisher_0.2.0 GenomeInfoDbData_1.2.4
## [10] cellranger_1.1.0 yaml_2.2.1 numDeriv_2016.8-1.1
## [13] pillar_1.4.6 lattice_0.20-41 RUnit_0.4.32
## [16] digest_0.6.27 XVector_0.30.0 rbibutils_1.3
## [19] sandwich_3.0-0 htmltools_0.5.0 Matrix_1.2-18
## [22] pkgconfig_2.0.3 haven_2.3.1 bookdown_0.21
## [25] zlibbioc_1.36.0 mvtnorm_1.1-1 openxlsx_4.2.3
## [28] rio_0.5.16 tibble_3.0.4 car_3.0-10
## [31] ellipsis_0.3.1 TH.data_1.0-10 mnormt_2.0.2
## [34] survival_3.2-7 magrittr_1.5 crayon_1.3.4
## [37] readxl_1.3.1 evaluate_0.14 MASS_7.3-53
## [40] forcats_0.5.0 xml2_1.3.2 foreign_0.8-80
## [43] ropls_1.22.0 tools_4.0.3 data.table_1.13.2
## [46] hms_0.5.3 gbRd_0.4-11 lifecycle_0.2.0
## [49] multcomp_1.4-14 mutoss_0.1-12 stringr_1.4.0
## [52] glmnet_4.0-2 plotrix_3.7-8 zip_2.1.1
## [55] DelayedArray_0.16.0 pls_2.7-3 compiler_4.0.3
## [58] rlang_0.4.8 grid_4.0.3 RCurl_1.98-1.2
## [61] iterators_1.0.13 CompQuadForm_1.4.3 bitops_1.0-6
## [64] rmarkdown_2.5 multtest_2.46.0 codetools_0.2-16
## [67] abind_1.4-5 curl_4.3 zoo_1.8-8
## [70] knitr_1.30 mathjaxr_1.0-1 shape_1.4.5
## [73] metap_1.4 stringi_1.5.3 Rcpp_1.0.5
## [76] vctrs_0.3.4 xfun_0.18