mbQTL 1.0.0
knitr::opts_chunk$set(
dpi = 300,
warning = FALSE,
collapse = TRUE,
error = FALSE,
comment = "#>",
echo = TRUE
)
library(mbQTL)
library(tidyr)
data("SnpFile")
data("microbeAbund")
data("CovFile")
data("metagenomeSeqObj")
doctype <- knitr::opts_knit$get("rmarkdown.pandoc.to")
sessionInfo()
#> R version 4.3.0 RC (2023-04-13 r84269)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB 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
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] tidyr_1.3.0 mbQTL_1.0.0 BiocStyle_2.28.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.3 shape_1.4.6 xfun_0.39
#> [4] bslib_0.4.2 ggplot2_3.4.2 caTools_1.18.2
#> [7] Biobase_2.60.0 lattice_0.21-8 vctrs_0.6.2
#> [10] tools_4.3.0 bitops_1.0-7 generics_0.1.3
#> [13] parallel_4.3.0 tibble_3.2.1 fansi_1.0.4
#> [16] pkgconfig_2.0.3 pheatmap_1.0.12 Matrix_1.5-4
#> [19] KernSmooth_2.23-20 RColorBrewer_1.1-3 readxl_1.4.2
#> [22] lifecycle_1.0.3 compiler_4.3.0 stringr_1.5.0
#> [25] gplots_3.1.3 munsell_0.5.0 codetools_0.2-19
#> [28] htmltools_0.5.5 sass_0.4.5 yaml_2.3.7
#> [31] glmnet_4.1-7 pillar_1.9.0 jquerylib_0.1.4
#> [34] MatrixEQTL_2.3 metagenomeSeq_1.42.0 cachem_1.0.7
#> [37] limma_3.56.0 iterators_1.0.14 foreach_1.5.2
#> [40] gtools_3.9.4 tidyselect_1.2.0 locfit_1.5-9.7
#> [43] digest_0.6.31 stringi_1.7.12 dplyr_1.1.2
#> [46] purrr_1.0.1 bookdown_0.33 splines_4.3.0
#> [49] fastmap_1.1.1 grid_4.3.0 colorspace_2.1-0
#> [52] cli_3.6.1 magrittr_2.0.3 survival_3.5-5
#> [55] utf8_1.2.3 broom_1.0.4 scales_1.2.1
#> [58] backports_1.4.1 Wrench_1.18.0 rmarkdown_2.21
#> [61] matrixStats_0.63.0 cellranger_1.1.0 evaluate_0.20
#> [64] knitr_1.42 rlang_1.1.0 Rcpp_1.0.10
#> [67] glue_1.6.2 BiocManager_1.30.20 BiocGenerics_0.46.0
#> [70] jsonlite_1.8.4 R6_2.5.1
mbQTL
is a statistical package for identifying significant taxa (microbial),
snp (genotype) relationships across large sample sets. The package measures
this relationship using three different approaches. A) linear regression analysis
which uses core MatrixeQTL
model. B) rho estimation for correlation of taxa and SNP,
which is a novel method of estimating rho values after R value of taxa-taxa is
estimated then adjusted for association of taxa-SNP. Using this method various
taxa-SNP associations can be estimated simultaneously, and linkage disequilibrium (LD)
snps and regions can be identified in terms of presence and their associations with
snps (see section for more details). C) Utilizing Logistic regression analysis, we
estimate the relationship of one snp and one taxa across the sample cohort but
simultaneously this could be estimated for all samples and proper statistics is
produced from these regressions. Appropriate plots are generated for all three
methods utilized in this vignette and are part of the package to simplify
visualization.
The datasets used for this vignette are all simulated and are not factual. The “Taxa_table_MVP_cor.txt” and “SNP_table_MVP_cor.txt” are used in all three sections mentioned above and the “covariates.txt” file is used for section A)linear regression. The input file format for “SNP_table_MVP_cor.txt” should be colnames with SNP names/locations/rs_numbers and rownames needs to be sample number. For the “Taxa_table_MVP_cor.txt” the colnames is the microbe name and the row name is the sample number and finally for the “covariates.txt” file the colnames need to be samplenames and the rownames need to be the covariates. Please ensure all 3 files are matching in terms of sample number and are appropriately formatted to avoid errors in your analysis.
mbQTL accepts data in dataframe, MRexperiment
objects or table formats. The files
are accessible when the mbQTL
library is loaded.
# microbial count file
# microbeAbund
# SNP file
# SnpFile
# Covariance file
# CovFile
metagenomeSeqObj
#> MRexperiment (storageMode: environment)
#> assayData: 1313 features, 93 samples
#> element names: counts
#> protocolData: none
#> phenoData
#> sampleNames: 2001 2002 ... 2091 (93 total)
#> varLabels: ID Sample_Type ... Spike_in (15 total)
#> varMetadata: labelDescription
#> featureData
#> featureNames: 1595f09bb28c3e40d9779b2ed6e0c1eb
#> cf3950eea7fd94bf5f51936de1c6bd5d ...
#> fabc9190fc0e035265780a5c5a589b62 (1313 total)
#> fvarLabels: OTUname Kingdom ... Species (8 total)
#> fvarMetadata: labelDescription
#> experimentData: use 'experimentData(object)'
#> Annotation:
If the user has an MRexperiment microbiome file (metagenomeSeq
formatted), one can
use the following function to normalize their metagenomeSeq
object or aggregate based
on various taxonomic levels before normalization.
The output will be formatted in compatible data frame format for mbQTL
. If one has
a phyloseq
S4 object, they can use the phyloseq_to_metagenomeSeq()
option of the
phyloseq
package and start from there.
compatible_metagenome_seq <- metagenomeSeqToMbqtl(metagenomeSeqObj,
norm = TRUE, log = TRUE,
aggregate_taxa = "Genus"
)
#> Default value being used.
# Check for the class of compatible_metagenome_seq
class(compatible_metagenome_seq)
#> [1] "data.frame"
Linear regression analysis to identify significant association between taxa and
SNPs across large sample sized. P values FDRs, Pvalue histogram plots and QQ-plots
are provided for various taxa and SNP regression analysis performed. microbeAbund
is the taxa abundance file, SnpFile us file across samples (0,1,2 (reference,
heterozygous and alternate genotype). MatrixeQTL
makes its estimations assuming
three different models of which we employ two, a) assumption of simple linear
regression for each gene and SNV pair The first approach is a very similar approach
to that of the expression quantitative trait loci association (eQTL). b) regression
model with covariates (Shabalin, 2012).
# perform linear regression analysis between taxa and snps across
# samples and produce a data frame of results.(covariate file is
# optional but highly recommended to avoid results which might be
# prone to batch effects and obtaining optimal biological relevance
# for finding snp - taxa relationships.) The input files are
# microbial abundance file (rownames taxa, colnames sample names) and
# SnpFile (rownames snps and colnames sample names) with optional "CovFile"
# (rownames covariates and colnames sample names) ( note all samples should
# be concordant in the three datasets).
LinearAnalysisTaxaSNP <- linearTaxaSnp(microbeAbund,
SnpFile,
Covariate = CovFile
)
#> Processing covariates
#> Task finished in 0.002 seconds
#> Processing gene expression data (imputation, residualization)
#> Task finished in 0.003 seconds
#> Creating output file(s)
#> Task finished in 0.008 seconds
#> Performing eQTL analysis
#> 100.00% done, 169 eQTLs
#> Task finished in 0.012 seconds
#>
# Create a histogram of P values across all snp - taxa
# linear regression estimations
histPvalueLm(LinearAnalysisTaxaSNP)
# Create a QQPlot of expected versus estimated P value for all all
# snp - taxa linear regression estimations
qqPlotLm(microbeAbund, SnpFile, Covariate = CovFile)
#> Processing covariates
#> Task finished in 0.001 seconds
#> Processing gene expression data (imputation, residualization)
#> Task finished in 0.002 seconds
#> Creating output file(s)
#> Task finished in 0.009 seconds
#> Performing eQTL analysis
#> 100.00% done, 0 eQTLs
#> No significant associations were found.
#> 4
#> Task finished in 0.004 seconds
#>
Estimate taxa-taxa correlation (microbeAbund (taxa abundance file)) and estimate R
(Strength of the relationship between taxa x
and y
). In other words we estimate
Goodness of fits, as a measure of a feature’s association with host genetics
(Significant SNPs) association with taxa (matching 16S results). In details we
measure the strength of relationship between taxa x
and taxa y
(R correlation value),
next we measure goodness of fits R2 and ask to what extent the variance in one taxa
explains the variance in a SNP, and finally we find clusters of taxa associated with
various SNPs across the dataset by evaluating and controlling for both these measures
together and finally estimating rho value.
# Estimate taxa SNP - taxa correlation R and estimate R2
# (To what extent variance in one taxa explains the
# variance in snp)
# First estimate R value for all SNP-Taxa relationships
correlationMicrobes <- coringTaxa(microbeAbund)
# Estimate the correlation value (R2) between a specific snp and all taxa.
one_rs_id <- allToAllProduct(SnpFile, microbeAbund, "chr1.171282963_T")
# Estimate the correlation value (R2) between all snps and all taxa.
for_all_rsids <- allToAllProduct(SnpFile, microbeAbund)
# Estimate rho value for all snps and all taxa present in the dataset.
taxa_SNP_Cor <- taxaSnpCor(for_all_rsids, correlationMicrobes)
# Estimate rho value for all snps and all taxa present in the
# dataset and pick the highest and lowest p values in the
# dataset to identify improtant SNP-Taxa relationships.
taxa_SNP_Cor_lim <- taxaSnpCor(for_all_rsids,
correlationMicrobes,
probs = c(0.0001, 0.9999)
)
## Draw heatmap of rho estimates "taxa_SNP_Cor_lim" is the
# output of taxaSnpCor().
# One can use other pheatmap() settings for extra annotation
mbQtlCorHeatmap(taxa_SNP_Cor_lim,
fontsize_col = 5,
fontsize_row = 7
)