{"id":442,"date":"2022-04-19T11:06:25","date_gmt":"2022-04-19T03:06:25","guid":{"rendered":"https:\/\/linguopeng.top\/?p=442"},"modified":"2022-04-19T11:13:54","modified_gmt":"2022-04-19T03:13:54","slug":"%e8%a1%a8%e5%9e%8b-%e5%be%ae%e7%94%9f%e7%89%a9-%e4%bb%a3%e8%b0%a2%e7%89%a9-%e7%bb%84%e5%ad%a6%e5%85%b3%e8%81%94%e5%88%86%e6%9e%90","status":"publish","type":"post","link":"https:\/\/linguopeng.top\/?p=442","title":{"rendered":"\u8868\u578b-\u5fae\u751f\u7269-\u4ee3\u8c22\u7269 \u7ec4\u5b66\u5173\u8054\u5206\u6790"},"content":{"rendered":"\n
### \r\n#Pedersen, H.K., Forslund, S.K., Gudmundsdottir, V. et al. A computational #framework to integrate high-throughput \u2018-omics\u2019 datasets for the identification #of potential mechanistic links. Nat Protoc 13, 2781\u20132800 (2018). #https:\/\/doi.org\/10.1038\/s41596-018-0064-z\n        <\/g>\n        <\/svg><\/span>\r\n### Procedure for three-pronged host-microbiome analysis pipeline\r\n###\r\n### This R script essentially reproduces results from the analysis of MetaHIT\r\n### data in Pedersen, H. K. et al. 'Human gut microbes impact host serum\r\n### metabolome and insulin sensitivity'. Nature 535, 376\u2013381 (2016)., and is\r\n### intended to allow adaptations for analogous work on other datasets.\r\n###\r\n### The sections in this script corresponds to procedural sections in the parent\r\n### protocol: 'A computational framework for conducting a three-pronged\r\n### association study to integrate high-throughput \u2018-omics\u2019 datasets for the\r\n### identification of potential mechanistic links'.\r\n### It was adapted from initial scripts by Helle Krogh Pedersen (HKP), Valborg\r\n### Gudmundsdottir (VG) and Henrik Bj\u00f8rn Nielsen (BN) by HKP, Sofia K. Forslund\r\n### (SKF), VG and Anders Petersen (AP). \r\n### Last version 26\/07\/2018.\r\n###\r\n\r\n###\r\n### Step 1 - Set working directory\r\n###\r\n### The code will look for input files and deposit output files in\r\n### subdirectories relative to a main directory. On your computer, create a\r\n### directory for the analysis, to remain generic (and without assumption on\r\n### user operating system) we here call it top\/. Then create subdirectories to\r\n### obtain the following directory hierarchy:\r\n###\r\n###   top\/ \r\n###     r-code\/   \/\/ containing this script and any others\r\n###     data\/     \/\/ containing input files (i.e. all files (unzipped) from the 'example_input.zip'-file).\r\n###     results\/  \/\/ location of output files\r\n###\r\n### Copy all R-code from the git repository to your top\/r-code\/ subdirectory (i.e.\r\n### all .R files). Then open the main R-script called protocol_main.R in e.g.\r\n### RStudio. \r\n### Finally, modify the following command to set the working directory to\r\n### your top\/ directory: setwd(\u201c~\/top\u201d)\r\n###\r\n\r\n#setwd (\"~\/top\") ### Change the path to your main working directory. \r\n\r\n###\r\n### Step 2 - Ensure availability of software packages and satisfaction of\r\n### dependencies\r\n###\r\n### The code in this section simply makes all script-wise dependencies and\r\n### settings available in the present workspace. The analysis was tested using R\r\n### 3.3.3. If packages are available under another version, it should run, but\r\n### specifics of the implementation of each package may change results slightly.\r\n### The provided r-code completes in ~1 hour on a MacBook Pro (2.9GHz quad- core\r\n### 7th-generation Intel Core i7 processor, 16GB 2133MHz LPDDR3 memory) and can\r\n### be paused at any point.\r\n###\r\n### Ensure the following packages (including their indirect dependencies other\r\n### packages needed for their compilation and operation, some of which should be\r\n### installed via Bioconductor) are installed and loaded correctly correctly (in\r\n### every case available via the built-in R and Bioconductor package managers).\r\n###\r\n\r\nlibrary (xlsx) ### Saving to spreadsheet\r\nlibrary (data.table) ### Fast read of large files into R\r\nlibrary (WGCNA) ### -\tClustering software. Previously reported work done using v1.34\r\nlibrary (flashClust) ### Clustering software\r\nlibrary (ppcor) ### Partial Spearman correlations, for confounder analysis. Previously reported work done using v1.0\r\nlibrary (gplots) ### Plotting\r\nlibrary (cowplot) ### Plotting; to arrange several plots on the same page\r\nlibrary (ggplot2) ### Plotting\r\nlibrary (plyr) ### Data transformations\r\n\r\n### If you have problems installing some packages (e.g. WGCNA), it is likely \r\n### because you are missing one or more of the Bioconductor packages. \r\n### First install the dependencies with the code below, then reattempt the \r\n### installation of WGCNA\r\n# source(\"https:\/\/bioconductor.org\/biocLite.R\")\r\n# biocLite(\"impute\")\r\n# biocLite(\"GO.db\")\r\n# install.packages(\"WGCNA)\r\n\r\n### The script was testing using the following version of packages:\r\n# R version 3.3.3 (2017-03-06)\r\n# Platform: x86_64-apple-darwin13.4.0 (64-bit)\r\n# Running under: macOS Sierra 10.12.6\r\n# \r\n# locale:\r\n# [1] da_DK.UTF-8\/da_DK.UTF-8\/da_DK.UTF-8\/C\/da_DK.UTF-8\/da_DK.UTF-8\r\n# \r\n# attached base packages:\r\n# [1] stats     graphics  grDevices utils     datasets  methods   base     \r\n# \r\n# other attached packages:\r\n# [1] plyr_1.8.4            cowplot_0.9.1         ggplot2_2.2.1         gplots_3.0.1          ppcor_1.1            \r\n# [6] MASS_7.3-47           flashClust_1.01-2     WGCNA_1.61            fastcluster_1.1.24    dynamicTreeCut_1.63-1\r\n# [11] data.table_1.10.4-3  xlsx_0.5.7            xlsxjars_0.6.1        rJava_0.9-9          \r\n# \r\n# loaded via a namespace (and not attached):\r\n# [1] B iobase_2.34.0        bit64_0.9-7           splines_3.3.3         foreach_1.4.3         gtools_3.5.0         \r\n# [6]  Formula_1.2-2         stats4_3.3.3          latticeExtra_0.6-28   blob_1.1.0            fit.models_0.5-14    \r\n# [11] yaml_2.1.14           robustbase_0.92-8     impute_1.48.0         RSQLite_2.0           backports_1.1.1      \r\n# [16] lattice_0.20-35       digest_0.6.12         RColorBrewer_1.1-2    checkmate_1.8.5       colorspace_1.3-2     \r\n# [21] htmltools_0.3.6       preprocessCore_1.36.0 Matrix_1.2-12         pcaPP_1.9-72          pkgconfig_2.0.1      \r\n# [26] GO.db_3.4.0           mvtnorm_1.0-6         scales_0.5.0          gdata_2.18.0          htmlTable_1.11.2     \r\n# [31] tibble_1.3.4          IRanges_2.8.2         nnet_7.3-12           BiocGenerics_0.20.0   lazyeval_0.2.1       \r\n# [36] survival_2.41-3       magrittr_1.5          memoise_1.1.0         doParallel_1.0.11     foreign_0.8-69       \r\n# [41] tools_3.3.3           matrixStats_0.52.2    stringr_1.2.0         S4Vectors_0.12.2      munsell_0.4.3        \r\n# [46] cluster_2.0.6         AnnotationDbi_1.36.2  caTools_1.17.1        rlang_0.1.4           grid_3.3.3           \r\n# [51] iterators_1.0.8       rstudioapi_0.7        htmlwidgets_0.9       robust_0.4-18         labeling_0.3         \r\n# [56] bitops_1.0-6          base64enc_0.1-3       gtable_0.2.0          codetools_0.2-15      DBI_0.7              \r\n# [61] reshape2_1.4.2        rrcov_1.4-3           gridExtra_2.3         knitr_1.17            bit_1.1-12           \r\n# [66] Hmisc_4.1-1           KernSmooth_2.23-15    stringi_1.1.6         parallel_3.3.3        Rcpp_0.12.13         \r\n# [71] rpart_4.1-11          acepack_1.4.1         DEoptimR_1.0-8       \r\n\r\n###\r\n### Step 3 - Import input files\r\n### \r\n### The following input data is assumed to exist within the top\/data subdirectory:\r\n###\r\n### Data-files: \r\n### -\tphenotypes.tab: \r\n###   File with clinical phenotypes (columns) per\r\n###   individual (rows). Used to test for associations with, or for confounder\r\n###   analysis. Individuals are labeled \u2018idv\u2019 followed by a number, e.g. \u2018idv001\u2019.\r\n###\r\n### -\tmetabolomic.tab \/ lipidomic.tab: \r\n###   Input data matrix for abundance of 325\r\n###   polar metabolites or 876 molecular lipids per individual. Note that no\r\n###   additional normalization is done in this script, so data is assumed to be\r\n###   comparable in these regards. Such different data types are eventually merged\r\n###   into a single set of metabolite cluster abundances. Individual\r\n###   metabolite\/lipids are named M or L (for specifying a polar metabolite or\r\n###   molecular lipid, respectively), followed by a number and lastly the\r\n###   annotation or \u2018unknown\u2019 in case of unannotated metabolites\/lipids, e.g.\r\n###   \u2018M_20_Valine\u2019.\r\n###\r\n### -\tMGS_abundance.tab: \r\n###   File with the abundance (e.g. median gene abundance) of\r\n###   MGSs (columns) per individual (rows). These are assumed to have been\r\n###   rarefied to comparable depth or otherwise normalized. For historical\r\n###   reasons, the MGSs are labeled \u2018T2DCAG\u2019 followed by a number, e.g.\r\n###   \u2018T2DCAG00001\u2019.\r\n###\r\n### -\tKO_abundance.tab: \r\n###   File with the abundance of each KO (columns) per\r\n###   individual (rows). The data is assumed to be rarefied to comparable depth or\r\n###   otherwise normalized. For this, the software tool rtk58 can be used.\r\n###\r\n### - gene_abundance_sub.tab: \r\n###   File with the abundance of each catalog gene\r\n###   (subset-version) in each individual assumed to be rarefied to comparable\r\n###   depth or otherwise normalized.\r\n### \r\n### Annotation-files:\r\n### -\tcluster_mapping_file.tab: \r\n###   Input file with annotation for metabolite\r\n###   clusters, as available from curation of data in the specific dataset. The\r\n###   WGCNA clustering algorithm names the generated clusters with color codes.\r\n###   This mapping file simply facilitate renaming to more meaningful cluster\r\n###   descriptions. Here, the serum polar metabolite and serum molecular lipid\r\n###   clusters are labelled M01\u2013M35 and L01\u2013L39, respectively, and collectively\r\n###   termed metabolite clusters.\r\n###\r\n### -\tMGS_taxonomy.tab: \r\n###   File with taxonomic annotation of the MGSs (rows) used\r\n###   in the analysis. Each row contains the following information:\r\n###   species_taxonomy, species_pct, genus_taxonomy, genus_pct, family_taxonomy,\r\n###   family_pct, order_taxonomy, order_pct, phylum_taxonomy, phylum_pct, where\r\n###   x_pct is the percentage of the MGS genes that can be annotated (by sequence\r\n###   similarity) to the taxonomy of the MGS. If no taxonomy can be assigned to\r\n###   the MGS, the value will be NA.\r\n###\r\n### -\tKEGG_modules.tab: \r\n###   File containing definition of KEGG gene functional\r\n###   modules (rows) specifying which KO gene groups constitute each KEGG module.\r\n###   The first two columns contain KEGG module entry (number) and name, the third\r\n###   column list all KOs separated by semicolon. Any other functional annotation\r\n###   used analogously could be swapped in instead of the name. This file can be\r\n###   obtained by downloading KEGG modules from\r\n###   www.genome.jp\/kegg-bin\/get_htext?ko00002.keg (www.genome.jp\/kegg\/ --> KEGG\r\n###   MODULE --> KEGG modules); download htext and then running the provided\r\n###   script \u2018parse_kegg.pl\u2019 (after changing input filename).\r\n###\r\n### -\tKO_to_MGS.tab: \r\n###   File listing for each KO (rows) the MGSs (space-separated)\r\n###   it is a member of. This is based on the KO annotation of the gene catalog\r\n###   (gene_to_KO.tab) and the information, what gene from the gene catalog is\r\n###   within each MGS as given in the list MGS_to_gene.tab.\r\n###\r\n### -\tgene_to_KO.tab: \r\n###   File containing the KO annotation (if any) per gene (rows)\r\n###   in the gene catalog (constituting 7,328,469 genes). Genes are labeled\r\n###   \u2018RefCat620\u2019 followed by a number (1\u20267328469), e.g. \u2018RefCat620.1\u2019. Used to\r\n###   create KO_to_MGS.tab. Note, for the purpose of this protocol and to\r\n###   considerably reduce size of input data, only the subset of catalogue genes\r\n###   with KO annotation (n = 2,205,769) are provided in files with gene abundance\r\n###   or annotation (gene_to_KO.tab, MGS_to_gene.tab and gene_abundance_sub.tab)\r\n###   as only those genes are used in the driver-species analysis.\r\n###\r\n### -\tMGS_to_gene.tab: \r\n###   List of genes binned into a given MGS (rows). Used to\r\n###   create KO_to_MGS.tab. Rather than gene names the file contains the index\r\n###   value of the gene. i.e. position of the gene in the gene catalogue\r\n###   (subset-version, thus 1\u20262205769).\r\n###\r\n### All demonstration files are tab-delimited text files, but other formats\r\n### would equally work after modifying the respective file-import commands in\r\n### the R-script. \r\n###\r\n### For all demonstration data, pseudonymised sample names were re-randomised to\r\n### generate anonymised data.\r\n###\r\n\r\noptions (stringsAsFactors = FALSE)\r\n\r\n### - phenotypes.tab\r\n\r\nphenotypes = read.table (file = \"data\/phenotypes.tab\", row.names = 1, header = T, sep = \"\\t\")\r\nctrl = rownames (subset (phenotypes, Diabetes == \"nonDiabetic\")) ### set of control samples for analyses restricted to non-diabetic individuals.\r\n\r\n### - metabolomic.tab\r\n\r\nmetabolomic = read.table (file = \"data\/metabolomic.tab\", row.names = 1, header = T, sep = \"\\t\")\r\n\r\n### - lipidomic.tab\r\n\r\nlipidomic = read.table (file = \"data\/lipidomic.tab\", row.names = 1, header = T, sep = \"\\t\")\r\n\r\n### - cluster_mapping_file.tab\r\n\r\ncluster_mapping_file = read.table (file = \"data\/cluster_mapping_file.tab\", header = T, sep = \"\\t\", row.names = 1)\r\ncluster_mapping_file$label =  sapply (rownames (cluster_mapping_file),  function (x) paste (cluster_mapping_file [x, \"New_Name\"], cluster_mapping_file [x, \"Description\"], sep = \": \"))\r\n\r\n### - MGS_abundance.tab\r\n\r\nmgs_abundance = read.table (file = \"data\/MGS_abundance.tab\", sep = \"\\t\", row.names = 1, header = T)\r\n\r\n### - MGS_taxonomy.tab\r\n\r\nmgs_taxonomy  = read.table (file = \"data\/MGS_taxonomy.tab\", sep = \"\\t\", row.names = 1, header = T)\r\n\r\n### - KEGG_modules.tab\r\n\r\ntmp = read.table (\"data\/KEGG_modules.tab\", sep = \"\\t\")\r\nkoann = strsplit (tmp[,3], split = \";\")\r\nnames (koann) = tmp[,1]\r\n\r\nmodule_mapping = tmp[,2] ### description of Kegg modules\r\nnames (module_mapping) = tmp[,1] ; rm (tmp)\r\n\r\n### remove KEGG references in square brackets for more clean names for plotting\r\nmodule_mapping_clean = sapply(module_mapping, function(x) strsplit(x, \" \\\\[\")[[1]][1])\r\n\r\n### - KO_abundance.tab\r\n\r\nko_abundance = read.table (file = \"data\/KO_abundance.tab\", sep = \"\\t\", row.names = 1, header = T)\r\n\r\n### - KO_to_MGS.tab\r\n\r\ntmp = read.table (\"data\/KO_to_MGS.tab\", sep = \"\\t\", strip.white = T)\r\nKO2MGS = strsplit (tmp[,2], split = \" \")\r\nnames (KO2MGS) = tmp[,1] ; rm (tmp)\r\n\r\n### - gene_to_KO.tab\r\n\r\ngene2KO = read.table (\"data\/gene_to_KO.tab\", sep = \"\\t\", row.names = 1, header = T)\r\nKO2gene = tapply (1:nrow (gene2KO), gene2KO[,1], c) ; rm (gene2KO)\r\n\r\n### - MGS_to_gene.tab\r\n\r\ntmp = read.table (\"data\/MGS_to_gene.tab\", sep = \"\\t\", strip.white = T)\r\nMGS2gene = strsplit (tmp[,2], split = \" \")\r\nnames (MGS2gene) = tmp[,1] ; rm (tmp)\r\n\r\n### - gene_abundance_sub.tab\r\n\r\ngene_abundance_sub = data.frame( fread (\"data\/gene_abundance_sub.tab\", sep = \"\\t\", header = T), row.names = 1)\r\n\r\n### \r\n### Step 4 - Preprocessing\/cleanup of loaded data for sparsity and domain\r\n### limitation\r\n###\r\n### The following commands restrict the input data to account for method\r\n### limitations.\r\n###\r\n\r\n### MGS sparsity filter step - exclude MGSs that occur in <3 of the control individuals\r\ntest = apply (mgs_abundance [intersect (ctrl, rownames (mgs_abundance)),], 2, function(x) length (x [x != 0])) \r\nincl = names (test [test >= 3])\r\nlength (incl)\r\nMGSs = incl\r\nmgs_abundance = mgs_abundance [,incl]\r\nmgs_taxonomy = mgs_taxonomy [incl,]\r\nrm (test, incl)\r\n### MGS sparsity filter step done\r\n\r\n### Filtering out KEGG modules that are eukaryotic only\r\neuk = unique (c (\"M00352\", \"M00355\", \"M00354\", \"M00285\", \"M00295\", \"M00341\", \r\n                 \"M00177\", \"M00160\", \"M00359\", \"M00391\", \"M00182\", \"M00340\", \r\n                 \"M00359\", \"M00182\", \"M00160\", \"M00177\", \"M00391\", \"M00180\",  ### \"eukaryotes\" in name\r\n                 \"M00351\", \"M00352\", \"M00355\", \"M00354\", \"M00353\", \"M00427\")) ### spliceosome or nuclear functions\r\nincl = setdiff (names (koann), euk)\r\nkoann = koann [incl] ; rm (incl)\r\n### Eukaryote filtering done\r\n\r\n### KO sparsity filter step - exclude KOs that occur in <3 of the control individuals\r\ntest = apply (ko_abundance [intersect (ctrl, rownames (ko_abundance)),], 2, function (x) length (x [x != 0]))  \r\nincl = names (test [test >= 3])\r\nlength (incl)\r\nko_abundance = ko_abundance [, incl]\r\nrm (test, incl)\r\n### KO sparsity filter step done\r\n\r\n### Determine control sample set and ensure there are no missing values for\r\n### phenotypes tested, as algorithms used cannot handle this well.\r\n\r\nctrl.no.na = ctrl [! is.na (phenotypes [ctrl, \"Homa.IR\"])]\r\n\r\n###\r\n### Step 5 - Identify clusters of polar metabolites\r\n###\r\n### A central part of the protocol is dimensionality reduction of the high-\r\n### dimensional input data. Here the space of metabolite measurements is reduced\r\n### by identifying clusters of co-abundant metabolite peaks, analogously to the\r\n### identification of MGSs as co-abundant gene groups. This is done using\r\n### algorithms developed for gene expression analysis, specifically WGCNA.\r\n###\r\n### Importantly, effective WGCNA performance depends on parameter settings for\r\n### cluster reconstruction which will vary to some extent between datasets.\r\n### Establishing these for a particular dataset requires some inspection of\r\n### resulting clustering behavior under parameter exploration. The full\r\n### procedure for this lies beyond the scope of the present Protocols, but the\r\n### reader is referred to the WGCNA main documentation.\r\n###\r\n### Thus, for optimal performance, the steps below should be performed for\r\n### different values and the resulting curves inspected, hereafter parameters\r\n### should be specified as appropriate. The examples below corresponds to\r\n### optimal choices for the dataset and analysis used in Pedersen et al., 2016.\r\n###\r\n### In this step, WGCNA clustering is performed separately on lipidomic and\r\n### metabolomic measurements to detect cluster of densely connected\r\n### metabolites\/lipids. The metabolite\/lipid profiles constituting a given\r\n### cluster are summarized by the 1st principal component of the\r\n### metabolite\/lipid abundance matrix (\u2018Module Eigen-metabolite\/lipid\u2019); i.e.\r\n### basically a weighted average abundance profile. Prior to this, optimal\r\n### parameters for WGCNA should be established for the dataset being analyzed.\r\n###\r\n\r\n### Settings for WGCNA generally\r\ncor_method          = \"spearman\" ### for association with clinical parameters\r\ncorFun_tmp          = \"bicor\"\r\ncluster_method      = \"average\"\r\ncorOptions_list     = list (use = 'pairwise.complete.obs') \r\ncorOptions_str      = \"use = 'pairwise.complete.obs'\"\r\nBH_pval_asso_cutoff = 0.05\r\nNetworkType         = \"signed\" ### Signed-network (as the PC1 and median profile does not make sense as a summary measure of a cluster with anticorrelated metabolites.)\r\n\r\n### Specify data and parameters \r\nID_common_all_samples_meta = intersect (rownames (phenotypes), rownames (metabolomic)) ### only include individuals with information for both domains.\r\ndat = as.data.frame (metabolomic [ID_common_all_samples_meta,])\r\ndat_tmp = log2 (dat) ### work in logarithmic space\r\n\r\n### Settings for WGCNA on polar metabolite measurements\r\n### Once these are established the steps below can be run\r\nRsquareCut_val      = 0.89 ### usually ranges 0.80-0.95 but requires inspecting curves\r\nmergingThresh       = 0.20 ### Maximum dissimilarity of module eigengenes (i.e. 1-correlation) for merging modules.\r\nminModuleSize       = 3 ### minimum number of metabolites constituting a cluster\r\nSoftPower           = 13 ### beta-value, main parameter to optimize\r\n\r\n### Calculate weighted adjacency matrix\r\nA = adjacency (dat_tmp, power = SoftPower, type = NetworkType, corFnc = corFun_tmp, corOptions = corOptions_str)\r\ncolnames (A) = rownames (A) = colnames (dat_tmp)\r\n### Define dissimilarity based on topological overlap\r\ndissTOM = TOMdist (A, TOMType = NetworkType)\r\ncolnames (dissTOM) = rownames (dissTOM) = colnames (dat_tmp)\r\n### Hierarchical clustering\r\nmetaTree = flashClust (as.dist (dissTOM), method = cluster_method)\r\n### Define modules by cutting branches\r\nmoduleLabels1 = cutreeDynamic (dendro = metaTree, distM = dissTOM, method = \"hybrid\", deepSplit = 4, pamRespectsDendro = T, minClusterSize = minModuleSize)\r\nmoduleLabels1 = labels2colors (moduleLabels1)\r\n### Automatically merge highly correlated modules\r\nmerge = mergeCloseModules (dat_tmp, moduleLabels1, corFnc = corFun_tmp, corOptions = corOptions_list, cutHeight = mergingThresh)\r\n### Determine resulting merged module colors\r\nmoduleLabels2 = merge$colors\r\n### Establish eigengenes of the newly merged modules, used for cluster overall abundances\r\nMEs = merge$newMEs\r\n### Choose final module assignments\r\nmoduleColorsMeta = moduleLabels2\r\nnames (moduleColorsMeta) = colnames (dat_tmp)\r\nMEsMeta = orderMEs (MEs)\r\nrownames (MEsMeta) = rownames (dat_tmp)\r\n\r\n### Determine relevant descriptive statistics of established clusters\r\n### kIN: within-module connectivity, determined by summing connectivity with all\r\n###      other metabolites in the given cluster.\r\n### kME: bicor-correlation between the metabolite profile and module eigenvector; \r\n### both measures of intramodular hub-metabolite status.\r\nkIN <-      vector (length = ncol (dat_tmp)); names (kIN) = colnames (dat_tmp)\r\nkME <-      vector (length = ncol (dat_tmp)); names (kME) = colnames (dat_tmp)\r\nmodules <-  vector (length = ncol (dat_tmp)); names (modules) = colnames (dat_tmp)\r\n\r\nfor (module in names (table (moduleColorsMeta))) {   \r\n\r\n\tall.metabolites = names (dat_tmp)\r\n\tinModule = (moduleColorsMeta == module)\r\n\tmodule.metabolites = names (moduleColorsMeta [inModule])\r\n  modules [module.metabolites] = module \r\n\tkIN [module.metabolites] = sapply (module.metabolites, function (x) sum (A [x, module.metabolites]) - 1)\r\n  datKME = signedKME (dat_tmp, MEsMeta, corFnc = corFun_tmp, corOptions = corOptions_str)\r\n\trownames (datKME) = colnames (dat_tmp)\r\n\tkME [module.metabolites] = datKME [module.metabolites, paste (\"kME\", module, sep = \"\")]   \r\n\r\n}\r\noutput = data.frame (\"module\" = modules, \"kME\" = kME, \"kIN\" = kIN, \"cluster_name\" = sapply (modules, function (m) cluster_mapping_file [paste0 (\"M_ME\", m), \"New_Name\"]))\r\n\r\n###\r\n### Step 6 - Link individual polar metabolites to phenotype of interest\r\n###\r\n### The dimensionality reduction approach (detailed in Stage 3 (Step 8-12))\r\n### hinges on identifying features (e.g. metabolites and lipids clusters, etc.)\r\n### that are associated with a host phenotype of interest. In the example work\r\n### we describe here, this was insulin resistance as assessed by the HOMA-IR\r\n### measurement.\r\n###\r\n### In this step, for reference, the individual polar metabolites are linked\r\n### with HOMA-IR, both directly and under adjustment for a potential confounder\r\n### variable. In this case, such de-confounding was done for body mass index\r\n### (BMI). In case of a binary phenotype variable, one can, for example,\r\n### substitute the Spearman correlation test with a Mann-Whitney U (MWU) test.\r\n###\r\n\r\ntmpMat = array (NA, c (ncol (metabolomic), 2, 2))\r\ndimnames (tmpMat) [[1]] = colnames (metabolomic)\r\ndimnames (tmpMat) [[2]] = c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\") \r\ndimnames (tmpMat) [[3]] = c (\"estimate\", \"p.value\")\r\n\r\n### Associating individual polar metabolites with HOMA-IR without de-confounding\r\n### for BMI\r\ntmpMat [, \"HOMA.ir\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (metabolomic [ctrl.no.na,], MARGIN = 2, FUN = function (x) \r\n\t  unlist (cor.test (phenotypes [ctrl.no.na, \"Homa.IR\"], x,\r\n\t\tmethod = cor_method, use = \"pairwise.complete.obs\")[c (\"estimate\", \"p.value\")])))       \r\n\r\n### Associating individual polar metabolites with HOMA-IR while de-confounding\r\n### for BMI\r\ntmpMat [, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (metabolomic [ctrl.no.na,], MARGIN = 2, FUN = function (x) \r\n\t  unlist (pcor.test (x = x, y = phenotypes [ctrl.no.na, \"Homa.IR\"], z = phenotypes [ctrl.no.na, \"BMI.kg.m2\"],\r\n\t\tmethod = cor_method) [c (\"estimate\", \"p.value\")])))       \r\n\r\n### Sort by cluster_name and then decreasing values of kIN\r\noutput2 = cbind (tmpMat[, \"HOMA.ir\", c (\"estimate\", \"p.value\")], \"p.adjust\" = p.adjust (tmpMat [, \"HOMA.ir\", \"p.value\"], method = \"BH\"), \r\n                 tmpMat[, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")], \"p.adjust\" = p.adjust (tmpMat [, \"HOMA.ir_BMI.adj.partial\", \"p.value\"], method = \"BH\"))\r\ncolnames (output2) = paste (rep (c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\"), each = 3), colnames (output2), sep = \"_\")\r\noutput3 = cbind (output, output2)\r\noutput3 = output3 [with (output3, order (cluster_name, -kIN)), ]\r\n\r\n### Write to file\r\nwrite.table (output3, file = \"results\/individual_metabolites.txt\", sep = \"\\t\", col.names = NA, quote = F, row.names = T)\r\nrm (output, output2, output3, tmpMat, dat, dat_tmp)\r\n\r\n###\r\n### Step 7a - Identify clusters of molecular lipid\r\n###\r\n### Analogous to Step 5. First identify optimal parameters for WGCNA, then\r\n### establish clusters.\r\n###\r\n\r\n### Specify data and parameters\r\ndat = as.data.frame (lipidomic [rownames (phenotypes),])\r\ndat_tmp = log2 (dat)\r\n\r\n### Settings for WGCNA on molecular lipids measurements\r\n### Once these are established the steps below can be run\r\nRsquareCut_val = 0.90\r\nmergingThresh = 0.25\r\nminModuleSize = 5\r\nSoftPower = 14\r\n\r\n### Calculate weighted adjacency matrix\r\nA = adjacency (dat_tmp, power = SoftPower, type = NetworkType, corFnc = corFun_tmp, corOptions = corOptions_str)\r\n### Define dissimilarity based on topological overlap\r\ndissTOM = TOMdist (A, TOMType = NetworkType)\r\n### Hierarchical clustering\r\nlipidTree = flashClust (as.dist (dissTOM), method = cluster_method)\r\n### Define modules by cutting branches\r\nmoduleLabels1 = cutreeDynamic (dendro = lipidTree, distM = dissTOM, method = \"hybrid\", deepSplit = 4, pamRespectsDendro = T, minClusterSize = minModuleSize)\r\nmoduleLabels1 = labels2colors (moduleLabels1)\r\n### Automatically merge highly correlated modules\r\nmerge = mergeCloseModules (dat_tmp, moduleLabels1, corFnc = corFun_tmp, corOptions = corOptions_list, cutHeight = mergingThresh)\r\n### Determine resulting merged module colors\r\nmoduleLabels2 = merge$colors\r\n### Establish eigengenes of the newly merged modules, used for cluster overall\r\n### abundances\r\nMEs = merge$newMEs\r\n### Choose module assignments\r\nmoduleColorsLipid = moduleLabels2\r\nnames (moduleColorsLipid) = colnames (dat_tmp)\r\nMEsLipid = orderMEs (MEs)\r\nrownames (MEsLipid) = rownames (dat_tmp)\r\n\r\n### Determine relevant descriptive statistics of established clusters\r\nkIN <-      vector (length = ncol (dat_tmp)); names (kIN) = colnames (dat_tmp)\r\nkME <-      vector (length = ncol (dat_tmp)); names (kME) = colnames (dat_tmp)\r\nmodules <-  vector (length = ncol (dat_tmp)); names (modules) = colnames (dat_tmp)\r\n \r\nfor (module in names (table (moduleColorsLipid))) {   \r\n\r\n\tall.lipids = names (dat_tmp)\r\n\tinModule = (moduleColorsLipid == module)\r\n\tmodule.lipids = names (moduleColorsLipid [inModule])\r\n  modules [module.lipids] = module \r\n\tkIN [module.lipids] = sapply (module.lipids, function (x) sum (A [x, module.lipids]) -1)\r\n\tdatKME = signedKME (dat_tmp, MEsLipid, corFnc = corFun_tmp, corOptions = corOptions_str)\r\n\trownames (datKME) = colnames (dat_tmp)\r\n\tkME [module.lipids] = datKME [module.lipids, paste (\"kME\", module, sep = \"\")]   \r\n\r\n}\r\noutput = data.frame (\"module\" = modules, \"kME\" = kME, \"kIN\" = kIN, \"cluster_name\" = sapply (modules, function (l) cluster_mapping_file [paste0 (\"L_ME\", l), \"New_Name\"]))\r\n\r\n###\r\n### Step 7b - Link individual molecular lipids to phenotype of interest\r\n###\r\n### Analogous to Step 6.\r\n### The resulting metabolite and lipid clusters are thereafter merged into a\r\n### combined dataset, collectively termed 'metabolite clusters', for downstream\r\n### analyses.\r\n###\r\n\r\ntmpMat = array (NA, c (ncol (lipidomic), 2,  2))\r\ndimnames (tmpMat) [[1]] = colnames (lipidomic)\r\ndimnames (tmpMat) [[2]] = c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\") \r\ndimnames (tmpMat) [[3]] = c (\"estimate\", \"p.value\")\r\n\r\n### Associating individual molecular lipids with HOMA-IR without de-confounding\r\n### for BMI\r\ntmpMat [, \"HOMA.ir\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (lipidomic [ctrl.no.na,], MARGIN = 2, FUN = function (x) \r\n\t  unlist (cor.test (phenotypes [ctrl.no.na, \"Homa.IR\"], x,\r\n\t\tmethod = cor_method, use = \"pairwise.complete.obs\") [c (\"estimate\", \"p.value\")])))       \r\n\r\n### Associating individual molecular lipids with HOMA-IR while de-confounding\r\n### for BMI\r\ntmpMat[, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (lipidomic [ctrl.no.na,], MARGIN = 2, FUN = function (x) \r\n\t  unlist (pcor.test (x = x, y = phenotypes [ctrl.no.na, \"Homa.IR\"], z = phenotypes [ctrl.no.na, \"BMI.kg.m2\"],\r\n\t\tmethod = cor_method) [c (\"estimate\", \"p.value\")])))       \r\n\r\n### Sort by cluster_name and then decreasing values of kIN\r\noutput2 = cbind (tmpMat[, \"HOMA.ir\", c (\"estimate\", \"p.value\")], \"p.adjust\" = p.adjust (tmpMat [, \"HOMA.ir\", \"p.value\"], method = \"BH\"), \r\n                 tmpMat[, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")], \"p.adjust\" = p.adjust (tmpMat [, \"HOMA.ir_BMI.adj.partial\", \"p.value\"], method = \"BH\"))\r\ncolnames (output2) = paste (rep (c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\"), each = 3), colnames (output2), sep = \"_\")\r\noutput3 = cbind (output, output2)\r\noutput3 = output3 [with (output3, order (cluster_name, -kIN)), ]\r\n\r\n### Write to file\r\nwrite.table (output3, file = \"results\/individual_lipids.txt\", sep = \"\\t\", col.names = NA, quote = F, row.names = T)\r\nrm (output, output2, output3, tmpMat, dat, dat_tmp)\r\n\r\n### Create a joint data frame of metabolite\/lipid cluster eigengene\r\n### equivalents\/effective abundances ('MEsMetLip')\r\nMEsMetLip = cbind (MEsMeta [rownames (MEsMeta),], MEsLipid [rownames (MEsMeta),])\r\n### Rename module names (columnames) from 'colors' to numbers\r\ncolnames (MEsMetLip) = cluster_mapping_file [c (paste0 (\"M_\", colnames (MEsMeta)), paste0 (\"L_\", colnames (MEsLipid))), \"New_Name\"]\r\n### Exclude the two bin-clusters (i.e. \"M_remaining\" and \"L_remaining\")\r\nMEsMetLip = subset (MEsMetLip, select = c (-M_remaining, -L_remaining))\r\n### Save MEsMetLip to file, order by module number\r\nwrite.table (MEsMetLip [ , order (colnames (MEsMetLip))], file = \"results\/MEs_metabolite_clusters.txt\", sep = \"\\t\", row.names = T, col.names = NA, quote = F)\r\n\r\n###\r\n### Step 8 - Link metabolite clusters to phenotype of interest\r\n###\r\n### Analogous to Step 6 (and 7b). \r\n### This is a core analysis step generating associations between the\r\n### integrated\/clustered \u2013omics data and a clinically interesting phenotype. In\r\n### Pedersen et al., 2016, this was insulin resistance (HOMA-IR measurement),\r\n### but any phenotype is possible, as is checking against other \u2013omics spaces or\r\n### overall \u2013omics measurements such as gut diversity or enterotype. This\r\n### analysis can further be conducted controlling for confounders such as BMI in\r\n### the analysis we previously reported, by performing tests with partial\r\n### correlations or extended to binary phenotype variables by substituting tests\r\n### of Spearman correlation with e.g. MWU tests.\r\n###\r\n\r\ncor_HOMA.IR <- list () ### Data structure for storing results of correlation tests under different setups\r\n\r\ntmpMat = array (NA, c (ncol (MEsMetLip), 2, 2))\r\ndimnames (tmpMat) [[1]] = names (MEsMetLip)\r\ndimnames (tmpMat) [[2]] = c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\") \r\ndimnames (tmpMat) [[3]] = c (\"estimate\", \"p.value\")\r\n\r\n### Associating metabolite clusters with HOMA-IR without de-confounding for BMI    \r\ntmpMat [, \"HOMA.ir\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (MEsMetLip [ctrl.no.na,], MARGIN = 2, FUN = function (x) \r\n\t  unlist (cor.test (phenotypes [ctrl.no.na, \"Homa.IR\"], x,\r\n\t\tmethod = cor_method, use = \"pairwise.complete.obs\") [c (\"estimate\", \"p.value\")])))       \r\n\r\n### Associating metabolite clusters with HOMA-IR while de-confounding for BMI    \r\ntmpMat [, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (MEsMetLip [ctrl.no.na,], MARGIN = 2, FUN = function (x) \r\n\t  unlist (pcor.test (x = x, y = phenotypes [ctrl.no.na, \"Homa.IR\"], z = phenotypes [ctrl.no.na, \"BMI.kg.m2\"],\r\n\t\tmethod = cor_method) [c (\"estimate\", \"p.value\")])))       \r\n\r\ncor_HOMA.IR [[\"metlip\"]] <- tmpMat\r\nrm (tmpMat)\r\n\r\n###\r\n### Step 9 - Link MGS metagenomic entities to phenotype of interest\r\n###\r\n### Analogous to Step 6 but for metagenomic taxonomic data.\r\n###\r\n\r\ntmpMat = array (NA, c (ncol (mgs_abundance), 2, 2))\r\ndimnames (tmpMat) [[1]] = colnames (mgs_abundance)\r\ndimnames (tmpMat) [[2]] = c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\")\r\ndimnames (tmpMat) [[3]] = c (\"estimate\", \"p.value\")\r\n\r\n### MGS sparsity filter step - exclude MGSs that occur in <3 of the control\r\n### individuals that are actually used in the analysis. Repeat the step since we\r\n### have excluded control individuals with missing information for HOMA-IR\r\nctrl.no.na.4MGS = intersect (ctrl.no.na, rownames (mgs_abundance)) # length=275\r\ntest = apply (mgs_abundance [intersect (ctrl.no.na.4MGS, rownames (mgs_abundance)),], 2, function (x) length (x [x != 0]))\r\nmgs.subset = names (test[test>=3])\r\nrm (test)\r\n### Sparsity filter done\r\n\r\n### Associating MGSs with HOMA-IR without de-confounding for BMI    \r\ntmpMat [mgs.subset, \"HOMA.ir\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (mgs_abundance [ctrl.no.na.4MGS, mgs.subset], MARGIN = 2, FUN = function (x) \r\n\t  unlist (cor.test (phenotypes [ctrl.no.na.4MGS, \"Homa.IR\"], x,\r\n\t\tmethod = cor_method, use = \"pairwise.complete.obs\") [c (\"estimate\", \"p.value\")])))       \r\n\r\n### Associating MGSs with HOMA-IR while de-confounding for BMI    \r\ntmpMat [mgs.subset, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")] =\r\n\tt (apply (mgs_abundance [ctrl.no.na.4MGS, mgs.subset], MARGIN = 2, FUN = function (x) \r\n\t  unlist (pcor.test (x = x, y = phenotypes [ctrl.no.na.4MGS, \"Homa.IR\"], z = phenotypes [ctrl.no.na.4MGS, \"BMI.kg.m2\"],\r\n\t\tmethod = cor_method) [c (\"estimate\", \"p.value\")])))       \r\n\r\ncor_HOMA.IR [[\"MGSs\"]] <- tmpMat\r\nrm (tmpMat)\r\n\r\n###\r\n### Step 10 - Link KEGG functions to phenotype of interest\r\n###\r\n### Analogous in goal to Step 6 but for metagenomics functional data. \r\n### Here we use KEGG modules, but any other groupings of genes into functional\r\n### modules could similarly be used (see examples in Table 1 in the accompanying\r\n### Protocol). However, it is more complex in that each KEGG module is\r\n### constituted by multiple KOs. Thus, to generate results on the level of\r\n### modules, a test is made if correlations between the phenotype and the\r\n### abundances of KOs in the module is significantly higher or lower (MWU test)\r\n### for the module member KOs than for all other KOs, thus also considering\r\n### module completeness beyond the single gene level. In case of a binary\r\n### phenotype variable, the KOs can instead (of Spearman correlation\r\n### coefficients) be ranked based on Wald statistics for testing differentially\r\n### abundant KOs with a negative binomial test with the DESeq2 R package using\r\n### non-rarefied gene counts.\r\n###\r\n\r\ntmpMat = array (NA, c (length (koann), 3, 2))\r\ndimnames (tmpMat) [[1]] = names (koann)\r\ndimnames (tmpMat) [[2]] = c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\", \"HOMA.ir_RichnessGenes7M.adj.partial\")\r\ndimnames (tmpMat) [[3]] = c (\"estimate\", \"p.value\")\r\n\r\n### KOs sparsity filter step - exclude KOs that occur in <3 of the control\r\n### individuals that are actually used in the analysis. Repeat the step since we\r\n### have excluded control individuals with missing information for HOMA-IR\r\nctrl.no.na.4KOs = intersect (ctrl.no.na, rownames (ko_abundance)) # length=275\r\ntest = apply (ko_abundance [intersect (ctrl.no.na.4KOs, rownames (ko_abundance)),], 2, function (x) length (x [x != 0]))\r\nko.subset = names (test [test >= 3])\r\nrm (test)\r\n### Sparsity filtering done\r\n\r\n### Associating KOs with HOMA-IR without de-confounding for BMI    \r\n### First, obtain the correlation coefficients for Spearman correlation between\r\n### KOs and HOMA-IR\r\nKO_cor = apply (ko_abundance [ctrl.no.na.4KOs, ko.subset], MARGIN = 2, FUN = function (x) \r\n\tcor (phenotypes [ctrl.no.na.4KOs, \"Homa.IR\"], x,\r\n\tmethod = cor_method, use = \"pairwise.complete.obs\"))      \r\nKO_cor_HOMA = KO_cor\r\n### Then, test for difference in correlation coefficients between KOs in the\r\n### KEGG module and all other KOs.\r\nfor (k in names (koann)) {\r\n\tincat =    na.omit (KO_cor [   names (KO_cor) %in% koann [[k]] ]) ### select all correlations between HOMA-IR and KOs in the KEGG module\r\n\tnotincat = na.omit (KO_cor [! (names (KO_cor) %in% koann [[k]])]) ### select all correlations between HOMA-IR and KOs NOT in the KEGG module\r\n\tif (length (incat) > 0 & length (notincat) > 0) {\r\n\t\tx = wilcox.test (incat, notincat)\r\n\t\ttmpMat [k, \"HOMA.ir\", \"p.value\"] = x$p.value\r\n\t\ttmpMat [k, \"HOMA.ir\", \"estimate\"] = (median (incat, na.rm = T) - median (notincat, na.rm = T))\r\n\t}\r\n}\r\nrm (KO_cor)\r\n\r\n### Associating KOs with HOMA-IR while de-confounding for BMI    \r\n### First, obtain the correlation coefficients for partial Spearman correlation\r\n### between KOs and HOMA-IR\r\nKO_cor = apply (ko_abundance [ctrl.no.na.4KOs, ko.subset], MARGIN = 2, FUN = function (x) \r\n\tpcor.test (x = x, y = phenotypes [ctrl.no.na.4KOs, \"Homa.IR\"], z = phenotypes [ctrl.no.na.4KOs, \"BMI.kg.m2\"],\r\n\tmethod = cor_method) [1, \"estimate\"])      \r\nKO_cor_HOMAadjBMI = KO_cor\r\n### Then, test for difference in correlation coefficients between KOs in the \r\n### KEGG module and all other KOs.\r\nfor (k in names (koann)) {\r\n\tincat =    na.omit (KO_cor [   names (KO_cor) %in% koann [[k]] ]) ### select all correlations between HOMA-IR.BMI.adj and KOs in the KEGG module\r\n\tnotincat = na.omit (KO_cor [! (names (KO_cor) %in% koann [[k]])]) ### select all correlations between HOMA-IR.BMI.adj and KOs NOT in the KEGG module\r\n\tif (length (incat) > 0 & length (notincat) > 0) {\r\n\t\tx = wilcox.test (incat, notincat)\r\n\t\ttmpMat [k, \"HOMA.ir_BMI.adj.partial\", \"p.value\"]  = x$p.value\r\n\t\ttmpMat [k, \"HOMA.ir_BMI.adj.partial\", \"estimate\"] = (median (incat, na.rm = T) - median (notincat, na.rm = T))\r\n\t}\r\n}\r\nrm (KO_cor)\r\n\r\n### Associating KOs with HOMA-IR while de-confounding for gut gene richness    \r\n### First, obtain the correlation coefficients for partial Spearman correlation\r\n### between KOs and HOMA.IR\r\nKO_cor = apply (ko_abundance [ctrl.no.na.4KOs, ko.subset], MARGIN = 2, FUN = function (x) \r\n\tpcor.test (x = x, y = phenotypes [ctrl.no.na.4KOs, \"Homa.IR\"], z = phenotypes [ctrl.no.na.4KOs, \"richnessGenes7M\"],\r\n\tmethod = cor_method) [1, \"estimate\"])      \r\nKO_cor_HOMAadjRichnessGenes7M = KO_cor\r\n### Then, test for difference in correlation coefficients between KOs in the \r\n### KEGG module and all other KOs.\r\nfor (k in names (koann)) {\r\n\tincat =    na.omit (KO_cor [   names (KO_cor) %in% koann [[k]] ]) ### select all correlations between HOMA-IR.RichnessGenes.adj and KOs in the KEGG module\r\n\tnotincat = na.omit (KO_cor [! (names (KO_cor) %in% koann [[k]])]) ### select all correlations between HOMA-IR.RichnessGenes.adj and KOs NOT in the KEGG module\r\n\tif (length (incat) > 0 & length (notincat) > 0) {\r\n\t\tx = wilcox.test (incat, notincat)\r\n\t\ttmpMat [k, \"HOMA.ir_RichnessGenes7M.adj.partial\", \"p.value\"]  = x$p.value\r\n\t\ttmpMat [k, \"HOMA.ir_RichnessGenes7M.adj.partial\", \"estimate\"] = (median (incat, na.rm = T) - median (notincat, na.rm = T))\r\n\t}\r\n}\r\nrm (KO_cor)\r\n\r\ncor_HOMA.IR [[\"keggmodules\"]] <- tmpMat\r\nrm (tmpMat, ko.subset, ctrl.no.na.4KOs)\r\n\r\n###\r\n### Step 11 - Save phenotype associations\r\n###\r\n### Save the (BMI corrected) HOMA-IR association of metabolite clusters, MGSs\r\n### and KEGG modules calculated in Step 8-10.\r\n###\r\n\r\ntmp.excel.file = \"results\/HOMA.IR_associations.xlsx\"\r\nwrite.xlsx (paste (\"Associations with HOMA-IR and HOMA-IR adjusted for BMI\"), \r\n            sheetName = \"info\", file = tmp.excel.file, row.names = F, col.names = F)\r\n\r\nfor (tmp_name in names (cor_HOMA.IR)) {\r\n  tmpMat = cor_HOMA.IR [[paste (tmp_name)]]\r\n  out = cbind (tmpMat [, \"HOMA.ir\", c(\"estimate\", \"p.value\")], \"p.adjust\" = p.adjust (tmpMat [, \"HOMA.ir\",\"p.value\"], method = \"BH\"), \r\n               tmpMat [, \"HOMA.ir_BMI.adj.partial\", c (\"estimate\", \"p.value\")], \"p.adjust\" = p.adjust (tmpMat [, \"HOMA.ir_BMI.adj.partial\", \"p.value\"], method = \"BH\"))\r\n  colnames (out) = paste (rep (c (\"HOMA.IR\", \"HOMA.IR_BMI.adj.partial\"), each = 3), colnames (out), sep = \"_\")\r\n  if (tmp_name == \"metlip\") { rownames (out) = cluster_mapping_file [match( rownames (out), cluster_mapping_file$New_Name), \"label\"] }\r\n  write.xlsx (x = out, file = tmp.excel.file, append = T, sheetName = paste (tmp_name, sep = \"_\"), row.names = T, col.names = T)            \r\n}\r\n\r\n###\r\n### Step 12 - Select features with significant differences\r\n###\r\n### Here, combine and integrate those functional, taxonomic and metabolomics\r\n### features which reliably correspond to the host phenotype of interest, then\r\n### later determine their inter-correlations.\r\n###\r\n\r\nfinal.fdr.cutoffs = 0.1 ### FDR thresholds, change to your likings.\r\n\r\nvartotest_union <- list ()\r\nvartotest_intersect <- list ()\r\nfor (tmp_name in names (cor_HOMA.IR)) {\r\n\ttmpMat = cor_HOMA.IR [[paste (tmp_name)]]\r\n\tvartotest_union [[paste (\"fdr\", final.fdr.cutoffs, sep = \"_\")]] [[tmp_name]] =\r\n\tunique (c (names (which (p.adjust (tmpMat [ , \"HOMA.ir\", \"p.value\"], method = \"BH\") < final.fdr.cutoffs)), \r\n\t\t   names (which (p.adjust (tmpMat [ , \"HOMA.ir_BMI.adj.partial\", \"p.value\"], method = \"BH\") < final.fdr.cutoffs))))\r\n\tvartotest_intersect [[paste(\"fdr\", final.fdr.cutoffs, sep = \"_\")]] [[tmp_name]] =\r\n\tintersect (names (which (p.adjust (tmpMat[ , \"HOMA.ir\", \"p.value\"], method = \"BH\") < final.fdr.cutoffs)),\r\n\t\t   names (which (p.adjust (tmpMat[ , \"HOMA.ir_BMI.adj.partial\", \"p.value\"], method = \"BH\") < final.fdr.cutoffs)))\r\n}\r\n\r\n###\r\n### Step 13 - Correlate metabolite clusters to functional metagenomic potentials\r\n### (KOs)\r\n###\r\n### The set of metabolite clusters associated with the phenotype of interest\r\n### should next be tested for association with the set of functional metagenome\r\n### features likewise so associated (identified in Step 12). Here once again the\r\n### complex nature of KEGG modules (consisting of multiple KOs) must be taken\r\n### into account.\r\n###\r\n### In this step, the correlations between each metabolite\/lipid cluster and\r\n### each KO are determined.\r\n###\r\n\r\nctrl.4KOs = intersect (ctrl, rownames (ko_abundance))\r\n\r\nKO_MetLip_cor = matrix (NA, nrow = ncol (ko_abundance), ncol = ncol (MEsMetLip))\r\nrownames (KO_MetLip_cor) = colnames (ko_abundance)\r\ncolnames (KO_MetLip_cor) = colnames (MEsMetLip)\r\n\r\nfor (m in colnames (MEsMetLip)) {\r\n\t\tKO_MetLip_cor [ , m] = apply (ko_abundance [ctrl.4KOs, ], MARGIN = 2, FUN = function (x) \r\n\t\t  cor (x, MEsMetLip [ctrl.4KOs, m],\r\n\t\t       method = cor_method, use = \"pairwise.complete.obs\"))\r\n}\r\n\r\n###\r\n### Step 14 - Associate metabolite clusters to functional metagenomic potentials\r\n### (KEGG modules)\r\n###\r\n### Using the KO level data generated in Step 13, to calculate module-level\r\n### associations between a KEGG module and a metabolite cluster.\r\n###\r\n\r\n### Select variables to include in the heatmap\r\n### note that for the metabolite clusters (metlip) we take the intersect, to \r\n### reduce the number of clusters in the heatmap, but for the KEGG modules \r\n### (keggmodules) we use the union.\r\nmetlip2test = vartotest_intersect [[paste(\"fdr\", final.fdr.cutoffs, sep = \"_\")]][[\"metlip\"]] \r\nkeggmodules2test = vartotest_union [[paste(\"fdr\", final.fdr.cutoffs, sep = \"_\")]][[\"keggmodules\"]] \r\n\r\nkegg_metlip_est = matrix (data = NA, nrow = length (keggmodules2test), ncol = length (metlip2test))\r\nrownames (kegg_metlip_est) = keggmodules2test\r\ncolnames (kegg_metlip_est) = metlip2test\r\nkegg_metlip_p = kegg_metlip_est\r\n\r\nfor (m in metlip2test) {\r\n\tfor (k in keggmodules2test) {\r\n\t\tincat =    na.omit (KO_MetLip_cor [   rownames (KO_MetLip_cor) %in% koann [[k]],  m]) ### select all correlations between metlip group and KOs in the KEGG module\r\n\t\tnotincat = na.omit (KO_MetLip_cor [! (rownames (KO_MetLip_cor) %in% koann [[k]]), m]) ### select all correlations between metlip group and KOs NOT in the KEGG module\r\n\t\tif (length (incat) > 0 & length (notincat) > 0) {\r\n\t\t\tx = wilcox.test (incat, notincat)\r\n\t\t\tkegg_metlip_p [k, m] = x$p.value\r\n\t\t\tkegg_metlip_est [k, m] = median (incat, na.rm = T) - median (notincat, na.rm = T)\r\n\t\t}\r\n\t}\r\n}\r\nrm (KO_MetLip_cor)\r\n\r\n###\r\n### Step 15 - Plot metabolome-microbiome functional analysis results\r\n###\r\n### Create visual representation of the generated results.\r\n###\r\n\r\n### load modified heatmap function from the gplots R-package, which allows for\r\n### both multi-column row-sidebar and showing text within cells (here\r\n### significance). Note this function (heatmap.3) is different from the\r\n### 'heatmap3' R package and the heatmap.3 function in GMD package The latest\r\n### version can be obtained from https:\/\/gist.github.com\/amcdavid\/5439787 We\r\n### have tested the script with the version loaded below.\r\nsource(\"r-code\/heatmap_3.R\")\r\n\r\n### BH adjust p-values for associations between KEGG modules and metabolite\r\n### modules\r\nkegg_metlip_p_adj = p.adjust (kegg_metlip_p, method = \"BH\")\r\ndim (kegg_metlip_p_adj) = dim (kegg_metlip_p)\r\nrownames (kegg_metlip_p_adj) = rownames (kegg_metlip_p)\r\ncolnames (kegg_metlip_p_adj) = colnames (kegg_metlip_p)\r\n\r\n### Select only significant associations to show in heatmap (to make it visually\r\n### comprehensible) I.e. exclude KEGG modules and metabolite modules with not at\r\n### least one significant association. In this example, all rows\/columns are\r\n### kept.\r\ntmp = kegg_metlip_p_adj < final.fdr.cutoffs\r\nissig = rowSums (tmp, na.rm = T)\r\nkegg = na.omit (names (issig [issig > 0]))\r\nlength (kegg)\r\nissig = colSums (tmp, na.rm=T)\r\nmetlip = na.omit (names (issig [issig > 0]))\r\nlength (metlip)\r\n\r\n### create matrix for heatmap\r\nplotmat = kegg_metlip_est [kegg, metlip]\r\n\r\n### create matrix with significance stars\r\nplotmat_p = kegg_metlip_p_adj [kegg, metlip] ### p-values to make stars for heatmap\r\nstars = matrix (\"\", ncol = ncol (plotmat_p), nrow = nrow (plotmat_p))\r\nrownames (stars) = rownames (plotmat_p)\r\ncolnames (stars) = colnames (plotmat_p)\r\nfor (z in 1:ncol (stars)) {\r\n  for (j in 1:nrow (stars)) {\r\n    if (plotmat_p [j, z] < 0.1) {\r\n      stars [j, z] = \"+\"\r\n    }\r\n    if (plotmat_p [j, z] < 0.01) {\r\n      stars [j, z] = \"*\"\r\n    }\r\n    if (plotmat_p [j, z] < 0.001) {\r\n      stars [j, z] = \"**\"\r\n    }\r\n  }\r\n}\r\nrm (plotmat_p)\r\n\r\n### make sidebar with phenotype associations\r\npheno = c (\"HOMA.ir\", \"HOMA.ir_BMI.adj.partial\")\r\nphenobar = matrix (NA, ncol = length (pheno), nrow = length (kegg))\r\ncolnames (phenobar) = pheno\r\nrownames (phenobar) = kegg\r\n\r\nfor (j in pheno) {\r\n  tmp = p.adjust (cor_HOMA.IR$keggmodules [ , j, \"p.value\"], method=\"BH\") ### adjust column wise on phenotypes\r\n  for (k in kegg) {\r\n    if (tmp [k] >= 0.1) {\r\n      phenobar [k, j] = \"grey\"\r\n    } else {\r\n      if (cor_HOMA.IR$keggmodules [k , j, \"estimate\"] > 0) {\r\n        phenobar [k, j] = \"darkblue\"\r\n      } else if (cor_HOMA.IR$keggmodules [k , j, \"estimate\"] < 0){\r\n        phenobar[k, j] = \"darkred\"\r\n      } else {\r\n        phenobar[k, j] = \"white\"\r\n      }\r\n    }\r\n  }\r\n}\r\n\r\n### Rename column names of phenobar\r\ncolnames (phenobar) = c(\"HOMA-IR\", \"HOMA-IR.BMI.adj\")\r\n\r\n### Note, normally one would cluster the rows and columns in the heatmap as\r\n### shown below but for the paper (Pedersen et al, 2016) we needed to group the\r\n### KEGG modules by biological similarity to make a higher-level annotation for\r\n### the figure and thus made a manual arrangement.\r\n\r\n### cluster rows of matrix \r\nd = dist (plotmat)\r\ndend = as.dendrogram (hclust (d, method = \"average\"))\r\nrm (d)\r\n\r\n### cluster columns of matrix \r\nd2 = dist (t (plotmat))\r\ndend2 = as.dendrogram (hclust (d2, method = \"average\"))\r\nrm (d2)\r\n\r\n### make heatmap with all the selected KEGG modules\r\npdf (\"results\/heatmap_KEGG_vs_metabolite_clusters.pdf\", height = 10, width = 10)\r\nheatmap.3 (plotmat,\r\n          Rowv = dend, \r\n          Colv = dend2, \r\n          dendrogram = \"column\",\r\n          col = rev (bluered (100)),\r\n          symbreaks = T,\r\n          key = T,\r\n          symkey = T,\r\n          keysize = 1,\r\n          KeyValueName = \"SCC.bg.adj\",\r\n          lhei = c(0.6, 4),\r\n          cellnote = stars, \r\n          notecol = \"black\",\r\n          notecex = 1.1,\r\n          trace = \"none\",\r\n          labRow = module_mapping_clean [rownames (plotmat)],\r\n          labCol = cluster_mapping_file [match( colnames (plotmat), cluster_mapping_file$New_Name), \"label\"],\r\n          RowSideColors = t (phenobar [rownames (plotmat), ]), \r\n          side.height.fraction = 0.35,    \r\n          NumColSideColors = dim (phenobar) [2], \r\n          margins = c (23, 33),\r\n          cexRow = 1,\r\n          cexCol = 1\r\n)\r\ndev.off () \r\n\r\n###\r\n### Step 16 (optional) - Export metabolome-microbiome associations for network\r\n### analysis\r\n###\r\n### Further exploration of such high dimensional association data can be\r\n### performed using software for network analysis, e.g. Cytoscape\r\n### (www.cytoscape.org) or igraph (www.igraph.org). Here we provide code for\r\n### exporting an edge-file with pairwise association scores and FDR-values\r\n### between metabolite clusters and KEGG modules as well as a corresponding node\r\n### attribute file; both in .txt-format. Several tutorials for importing,\r\n### visualizing and analyzing networks in Cytoscape can be found here:\r\n### https:\/\/github.com\/cytoscape\/cytoscape-tutorials\/wiki.\r\n###\r\n\r\n### Note, for simplicity, the node and edge files are made for a network\r\n### representing the heatmap plotted in Step 15 (i.e. KEGG modules in\r\n### \"keggmodules2test\" and metabolite clusters in \"metlip2test\")\r\n\r\n### Make edge file with the following four columns: \r\n### - KEGG_module, \r\n### - Metabolite_cluster, \r\n### - association estimate (SCC.bg.adj), \r\n### - FDR adjusted p-values\r\ntmp.edge.file.p_adj = melt (kegg_metlip_p_adj [keggmodules2test, metlip2test])\r\ntmp.edge.file.est   = melt (kegg_metlip_est [keggmodules2test, metlip2test])\r\ncolnames (tmp.edge.file.p_adj) = c (\"KEGG_module\", \"Metabolite_cluster\", \"p.adjust\")\r\ncolnames (tmp.edge.file.est)   = c (\"KEGG_module\", \"Metabolite_cluster\", \"SCC.bg.adj\")\r\nedge.file = cbind (tmp.edge.file.est, \"p.adjust\" = tmp.edge.file.p_adj$p.adjust)\r\nrm (tmp.edge.file.p_adj, tmp.edge.file.est)\r\n\r\n### Save the edge-file as a tab-seperated txt.file\r\nwrite.table (edge.file, file = \"results\/edge_file.txt\", sep = \"\\t\", row.names = F, col.names = T, quote = F)\r\n\r\n### Make node annotation file with the following three columns: \r\n### - Node_name (corresponds to names in the edge-file), \r\n### - Description (where existing), \r\n### - Examples (example metabolites, only for metabolite clusters)\r\n### Get metabolite annotation\r\nnode.file = cluster_mapping_file [cluster_mapping_file$New_Name %in% metlip2test, c (\"New_Name\", \"Description\", \"Examples\")]\r\ncolnames (node.file) = c (\"Node_name\", \"Description\", \"Examples\")\r\n### Add KEGG module annotation\r\nnode.file = rbind (node.file, data.frame (\"Node_name\" = keggmodules2test, \"Description\" = module_mapping_clean [keggmodules2test], \"Examples\" = rep (\"\", time = length (keggmodules2test))))\r\n\r\n### Save the node-file as a tab-seperated txt.file\r\nwrite.table (node.file, file = \"results\/node_file.txt\", sep = \"\\t\", row.names = F, col.names = T, quote = F)\r\n\r\n###\r\n### Step 17 - Leave-one-MGS-out analysis\r\n###\r\n### This part of the approach allows testing which bacterial taxa (sensu\r\n### metagenomics species (MGS) as defined from the metagenomics datasets\r\n### themselves) are driving the functional effects seen. I.e. enabling\r\n### assessment of the extent to which different taxa explain a functional\r\n### potential association to a phenotype of interest, such as insulin resistance\r\n### in the case of Pedersen et al., 2016. It is done by testing for each\r\n### functional feature (here: KEGG modules) to what extent leaving out each MGS\r\n### and the genes it contains causes a change in the association between those\r\n### modules and the target phenotype. Thus, the first step is setting which\r\n### modules and taxa will be tested, then performing this test for each\r\n### combination, after which results are plotted and reported.\r\n###\r\n\r\n### Specify the MGSs to be tested in the leave-on-MGS-out analysis. \r\n### Here it is set to all 788 MGSs.\r\nMGS.list = MGSs\r\n### Specify the KEGG modules to test in the leave-one-MGS-out analysis, and get\r\n### the KOs for those modules.\r\n### Here it is set to all modules defined as significant associated with HOMA-IR\r\n### and\/or HOMA-IR.bmi.adjusted (i.e. 'keggmodules2test').\r\nKOsets = koann [keggmodules2test] \r\n### Manually add the bcaa_biosynthesis 'mega-module' (relevant only for Pedersen\r\n### et al., 2016)\r\nKOsets [[\"bcaa_biosyn\"]] = unique (c (koann [[\"M00019\"]], koann [[\"M00570\"]], koann [[\"M00535\"]], koann [[\"M00432\"]]))\r\n\r\n### The following code can take long time.\r\n### for testing purposes, especially on older systems, we recommend subsetting \r\n### to a few KEGG modules and maybe also MGSs.\r\n# KOsets = list (\"M00060\"= KOsets[[\"M00060\"]],  \"bcaa_biosyn\"=KOsets[[\"bcaa_biosyn\"]])\r\n# MGS.list = c (\"T2DCAG00004\", \"T2DCAG00005\", \"T2DCAG00385\", \" T2DCAG00011\")\r\n\r\n### Functions: \r\n\r\n### Calculate abundance sum of ALL genes annotated to a KO\r\nKOprofile_FUN <- function (KO) {\r\n\tKOgenes_i <- unique (unlist (KO2gene [KO])) ### get all genes that are annotated to the KO\r\n\treturn (colSums (gene_abundance_sub [KOgenes_i,]))\r\n}\r\n\r\n### Calculate abundance sum of genes annotated to a KO while excluding the genes\r\n### that are part of 'LeftOutMGS'\r\nLOOKOprofile <- function (KO, LeftOutMGS) {\r\n\tif ( any (LeftOutMGS == unlist (KO2MGS [KO])) ){\tKOgenes_i <- setdiff (unique (unlist (KO2gene [KO])), MGS2gene [[LeftOutMGS]])} \r\n\telse { KOgenes_i <- unique (unlist (KO2gene [KO])) }\r\n\treturn (colSums (gene_abundance_sub [KOgenes_i,]))\r\n}\r\n\r\n### The main function calculating the effect of leaving-one-MGS-out on the \r\n### association between a KEGG module and a phenotype of interest\r\nLeaveOneMGSOut_FUN <- function (KOsets = KOsets, MGS.list = MGS.list, sample.list = sample.list, y, z = NULL) {\r\n\t### number of KOs kinds per MGS\r\n\tKO_types_per_MGS <- lapply (KOsets, function (KOset) {\r\n\t\tmgses <- sapply (KOset, function (KO) { intersect (unique (unlist (KO2MGS [KO])), MGS.list)})\r\n\t\ttable (unlist (mgses))})\r\n  \r\n\tif  (is.null (z) == T) { \r\n    \r\n\t\tprint (\"No 'z' is provided, calculating spearman correlation\")\r\n    \r\n\t\t### Code for spearman correlation, i.e. NOT adjusting for a confounding factor\r\n\t\ty = y [sample.list,]\r\n    \r\n\t\t### correlating the mean KO signal (calculated with 'KOprofile_FUN') to \r\n\t\t### HOMA-IR, and then taking the median SCC for the module\r\n\t\tprint (\"calculating spearman correlation, using all genes\")\r\n\t\tSCCallMGS<-sapply (KOsets, function (KOset) { \r\n      \t\t\tKOprofiles<-t (sapply (KOset, KOprofile_FUN))\r\n\t\t\treturn (median (apply (KOprofiles [rowSums (KOprofiles) > 0, sample.list], 1, function (x) { \r\n\t\t\t\tcor.test (x = x, y = y, method = \"spearman\", use = \"pairwise.complete.obs\", exact = FALSE)$estimate }))) \r\n\t\t})\r\n    \r\n\t\t### correlating the mean KO signal while excluding contribution from 'MGS' \r\n\t\t### (calculated with 'LOOKOprofile') to HOMA-IR, and then taking the median \r\n\t\t### SCC for the module.\r\n\t\t### Repeating for all MGS in 'MGS.list'\r\n\t\tprint (\"calculating spearman correlation, leaving-one-MGS-out\")\r\n\t\tSCComitingMGS <- lapply (KOsets, function (KOset) {\r\n\t\t\tsapply (intersect (unique (unlist (KO2MGS[KOset])),MGS.list), function (MGS){ \r\n\t\t\t\tKOprofiles <- (sapply (KOset, function (KO) { LOOKOprofile (KO, MGS) }))\r\n\t\t\t\treturn (median (cor (x = KOprofiles [sample.list, colSums (KOprofiles) > 0], y = y, method = \"spearman\", use = \"pairwise.complete.obs\"))) }) \r\n\t\t})\r\n    \r\n\t\t### Summarizing outout\r\n\t\tDeltaSCCperMGS <- lapply (names (SCComitingMGS), function (N) {\r\n\t\t\tSCC = SCCallMGS [[N]]\r\n\t\t\tSCC.bgadj =  (median (na.omit (KO_cor_HOMA [names (KO_cor_HOMA) %in% KOsets [[N]]]), na.rm = T) \r\n\t\t\t\t          - median (na.omit (KO_cor_HOMA[! (names (KO_cor_HOMA) %in% KOsets [[N]])]), na.rm = T)) \r\n\t\t\t# summary(KO_cor_HOMA); adjust % SCC effect for background distribution (i.e. the fact that median(SCC for HOMA vs modules) is negative and not = 0)\r\n\t\t\tdata.frame (SCC = SCC, SCC.bgadj = SCC.bgadj, SCC_omiting_MGS = SCComitingMGS [[N]], DeltaMGS_SCC = SCC - SCComitingMGS [[N]], \r\n\t\t\t\tpctSCCeffect = 100 * (SCC - SCComitingMGS [[N]]) \/ SCC, pctSCCeffect.bgadj = 100 * (SCC - SCComitingMGS [[N]]) \/ SCC.bgadj,\r\n\t\t\t\tDistinct_KOs_in_MGS = as.vector (KO_types_per_MGS [[N]][names (SCComitingMGS [[N]])]), row.names = names (SCComitingMGS [[N]])\r\n\t\t\t) [rev (order (pctSCCeffect = 100 * (SCC - SCComitingMGS[[N]]) \/ SCC)),] \r\n\t\t})\r\n\r\n\t\tnames (DeltaSCCperMGS) <- names (KOsets)\r\n\t\treturn (DeltaSCCperMGS)\r\n\r\n\t} \r\n\r\n\telse {\r\n    \r\n\t\tprint (\"'z' is provided, calculating partial spearman correlation, adjusting for 'z'\")\r\n\r\n\t\t### Code for partial spearman correlation, i.e. adjusting for a confounding factor\r\n\t\ty = y [sample.list,]\r\n\t\tz = z [sample.list,]\r\n\r\n\t\t### correlating the mean KO signal (calculated with 'KOprofile_FUN') to \r\n\t\t### HOMA-IR, and then taking the median SCC for the module\r\n\t\tprint (\"calculating partial spearman correlation, using all genes\")\r\n\t\tpartialSCCallMGS <- sapply (KOsets, function (KOset) {\t\r\n\t\t\tKOprofiles <- t (sapply (KOset, KOprofile_FUN))\r\n\t\t\treturn (median (apply (KOprofiles [rowSums (KOprofiles) > 0,], 1, function (x) { pcor.test (x = x [sample.list], y = y, z = z, method = \"spearman\")$estimate })))\r\n\t\t})\r\n\r\n\t\t### correlating the mean KO signal while excluding contribution from 'MGS' \r\n\t\t### (calculated with 'LOOKOprofile') to HOMA-IR, and then taking the median \r\n\t\t### SCC for the module.\r\n\t\t### Repeating for all MGS in 'MGS.list'\r\n\t\tprint (\"calculating partial spearman correlation, leaving-one-MGS-out\")\r\n\t\tpartialSCComitingMGS <- lapply (KOsets, function (KOset) {\r\n\t\t\tsapply (intersect (unique (unlist (KO2MGS [KOset])), MGS.list), function (MGS) { # looping over all MGGs in 'MGS.list'\r\n\t\t\t\tKOprofiles <- t (sapply (KOset, function (KO) {LOOKOprofile (KO, MGS) }))\r\n\t\t\t\treturn (median (apply (KOprofiles [rowSums (KOprofiles) > 0,], 1, function (x) { pcor.test (x = x [sample.list], y = y, z = z, method=\"spearman\")$estimate })))\r\n\t\t\t})\r\n\t\t})\r\n    \r\n\t\t### Summarizing outout\r\n\t\tpartialDeltaSCCperMGS<-lapply (names (partialSCComitingMGS), function (N){\r\n\t\t\tSCC = partialSCCallMGS[[N]]\r\n\t\t\tSCC.bgadj =  (median (na.omit (KO_cor_HOMAadjBMI[names (KO_cor_HOMAadjBMI) %in% KOsets[[N]]]), na.rm=T) \r\n\t\t\t\t          - median (na.omit (KO_cor_HOMAadjBMI[! (names (KO_cor_HOMAadjBMI) %in% KOsets[[N]])]), na.rm=T)) \r\n\t\t\t### summary(KO_cor_HOMA); adjust % SCC effect for background distribution (i.e. the fact that median(SCC for HOMA vs modules) is negative and not = 0) \r\n\t\t\tdata.frame (SCC =SCC, SCC.bgadj = SCC.bgadj, SCC_omiting_MGS = partialSCComitingMGS [[N]], DeltaMGS_SCC = partialSCCallMGS [[N]] - partialSCComitingMGS [[N]], \r\n\t\t\t\tpctSCCeffect = 100 * (SCC - partialSCComitingMGS [[N]]) \/ SCC, pctSCCeffect.bgadj = 100 * (SCC - partialSCComitingMGS [[N]]) \/ SCC.bgadj,\r\n\t\t\t\tDistinct_KOs_in_MGS = as.vector (KO_types_per_MGS [[N]][names (partialSCComitingMGS [[N]])]), \r\n\t\t\t\trow.names = names (partialSCComitingMGS [[N]])\r\n\t\t\t) [rev (order (pctSCCeffect = 100 * (SCC - partialSCComitingMGS [[N]]) \/ SCC)),] \r\n\t\t})\r\n\r\n\t\tnames (partialDeltaSCCperMGS) <- names (KOsets)\r\n\t\treturn (partialDeltaSCCperMGS)\r\n\r\n\t}  \r\n\r\n}\r\n### End functions\r\n\r\n### Next, run these functions to compute the delta SCC values for the test in\r\n### question - change here if partial correlation to account for a confounder is\r\n### desired.\r\n\r\n### Performing the Leave-one-MGS-out for HOMA-IR\r\nDeltaSCCperMGS <- LeaveOneMGSOut_FUN (KOsets = KOsets, MGS.list = MGS.list,\r\n\tsample.list = ctrl.no.na.4MGS, ### the subset of individuals with complete information for phenotype and metagenomic data.\r\n\ty = phenotypes [,\"Homa.IR\", drop = F]) \r\nsave (DeltaSCCperMGS, file = \"results\/delta_SCC_per_MGS.RData\")\r\n\r\n### Performing the Leave-one-MGS-out for HOMA-IR.bmi.adjusted\r\n# partialDeltaSCCperMGS = LeaveOneMGSOut_FUN (KOsets = KOsets, MGS.list = MGS.list,\r\n#\tsample.list = ctrl.no.na.4MGS, y = phenotypes [,\"Homa.IR\", drop = F], z = phenotypes [,\"BMI.kg.m2\", drop = F])\r\n# save (partialDeltaSCCperMGS, file = \"results\/partial_delta_SCC_per_MGS.RData\")\r\n\r\n###\r\n### Step 18 - Extracting top driver species from leave-one-MGS-out analysis for\r\n### each microbiome functional module\r\n###\r\n### Here, having computed the contribution per taxon to each module, we extract\r\n### top five drivers for interpretation.\r\n###\r\n\r\ntopX = 5\r\noutput = as.data.frame (matrix (NA, nrow = length (DeltaSCCperMGS), ncol= (2 + 4 * topX)))\r\ncolnames (output) = c (\"Module name\",\"Number of genes in KEGG_module\",\r\n\t\t\t\"top1\", \"top1_Number of module genes in MGS\", \"top1_DeltaMGS_SCC\",  \"top1_pctSCCeffect_bg.adj\",\r\n\t\t\t\"top2\", \"top2_Number of module genes in MGS\", \"top2_DeltaMGS_SCC\",  \"top2_pctSCCeffect_bg.adj\",\r\n\t\t\t\"top3\", \"top3_Number of module genes in MGS\", \"top3_DeltaMGS_SCC\",  \"top3_pctSCCeffect_bg.adj\",\r\n\t\t\t\"top4\", \"top4_Number of module genes in MGS\", \"top4_DeltaMGS_SCC\",  \"top4_pctSCCeffect_bg.adj\",\r\n\t\t\t\"top5\", \"top5_Number of module genes in MGS\", \"top5_DeltaMGS_SCC\",  \"top5_pctSCCeffect_bg.adj\")\r\nrownames (output) = names (DeltaSCCperMGS)\r\n\r\noutput [,\"Number of genes in KEGG_module\"] = sapply (rownames (output), function (i) length (KOsets [[i]]))\r\noutput [,\"Module name\"] = module_mapping [rownames (output)]\r\n\r\nfor (i in names (DeltaSCCperMGS)) {\r\n\r\n\tif (length (rownames (DeltaSCCperMGS[[i]])) <= topX ) {\r\n\r\n\t\ttmp.MGSs.top = rownames (DeltaSCCperMGS [[i]]) [1:topX]\r\n\t\ttmp.MGSs.top [DeltaSCCperMGS [[i]][tmp.MGSs.top, \"pctSCCeffect\"] <= 0] <- NA ### set top.MGSs to NA if removing them are NOT INCREASING the effect\r\n\t\ttmp.MGSs = c (tmp.MGSs.top)\r\n\r\n\t} \r\n\r\n\telse {\r\n\r\n\t\ttmp.MGSs.top = rownames (DeltaSCCperMGS [[i]]) [1:topX]\r\n\t\ttmp.MGSs.top [DeltaSCCperMGS [[i]][tmp.MGSs.top, \"pctSCCeffect\"] <= 0] <- NA ### set top.MGSs to NA if removing them are NOT INCREASING the effect\r\n\t\ttmp.MGSs = c (tmp.MGSs.top)\r\n\r\n\t}\r\n\r\n\t## Add taxonomically annotation (for those MGSs that have one)\r\n\ttmp.MGSs.names = sapply (tmp.MGSs, function (i) ifelse (i %in% rownames (mgs_taxonomy),\r\n\t\tpaste (\" \", i, \" : \", mgs_taxonomy [i, \"genus\"], \" sp \", sep = \"\"), # true\r\n\t\tpaste (\" \", i, \" \", sep = \"\") )) # false\r\n\ttmp.KOsInMGS      = DeltaSCCperMGS [[i]][ tmp.MGSs, \"Distinct_KOs_in_MGS\"]\r\n\ttmp.delta         = DeltaSCCperMGS [[i]][ tmp.MGSs, \"DeltaMGS_SCC\"]\r\n\ttmp.pct_bg.adj    = DeltaSCCperMGS [[i]][ tmp.MGSs, \"pctSCCeffect.bgadj\"]\r\n\toutput [i, 3:ncol (output)] = c (rbind (tmp.MGSs.names, tmp.KOsInMGS, tmp.delta, tmp.pct_bg.adj))\r\n\r\n\trm (tmp.delta, tmp.pct_bg.adj, tmp.MGSs.names, tmp.KOsInMGS, tmp.MGSs, tmp.MGSs.top)\r\n\r\n}\r\n\r\n### write 'output' to tab-delimited text file\r\nwrite.table (x = output, file = \"results\/top_driver_species.txt\", row.names = T, col.names = NA, quote = F, sep = \"\\t\")         \r\n\r\n###\r\n### Step 19 - Plotting leave-one-MGS-out results\r\n###\r\n### Having computed these results, we plot the top driver taxa for the\r\n### phenotype-associated gene functions (Figure 4).\r\n### These include sub-plots for:\r\n### - Distribution\/density plots of correlations for KOs in a KEGG module vs \r\n###   all other KOs\r\n### - Distribution of correlations when leave-one-MGS-out. \r\n### - Distribution of correlations when leave-one-MGS-out. bg.adj. SCC \r\n###   (i.e. what was shown in Figure 3c+d in Pedersen et al., 2016)\r\n### Median SCC for KO within a module (red) and all other remaining KOs (green)\r\n### are indicated in the first two plots.\r\n###\r\n\r\n### Specify pdf-file for output.\r\npdf (file = \"results\/density_plot_SCC_HOMA.IR.pdf\", width = 6 * 1, height = 3 * 3)\r\n\r\nfor (m in names (KOsets)) {\r\n  \r\n  ### Distribution of correlations for KOs in modules vs all other KOs\r\n  incat =     na.omit (KO_cor_HOMA [names (KO_cor_HOMA) %in% KOsets [[m]] ]) ### select all correlations between HOMA-IR and KOs in the KEGG module\r\n  incat2 =    data.frame (\"KO_cor_HOMA\" = incat, \"cat\" = \"KOs in module\")\r\n  notincat =  na.omit (KO_cor_HOMA [! (names (KO_cor_HOMA) %in% KOsets [[m]])]) ### select all correlations between HOMA-IR and KOs NOT in the KEGG module   \r\n  notincat2 = data.frame (\"KO_cor_HOMA\" = notincat, \"cat\" = \"KOs not in module\")\r\n  tmp.df = rbind (incat2, notincat2)\r\n  cdf <- ddply (tmp.df, \"cat\", summarise, rating.median = median (KO_cor_HOMA)) ### find median for each group\r\n  g1 <- ggplot (tmp.df, aes (x = KO_cor_HOMA, fill = cat)) + \r\n    geom_density (alpha = 0.4) +\r\n    geom_vline (data = cdf, aes (xintercept = rating.median, colour = cat), linetype = \"dashed\", size = 0.8) +\r\n    ggtitle (paste0 (\"KEGG module: \", m, \r\n                     \"\\n\", strsplit (module_mapping [m], split = \"\\\\[|,\")[[1]][1], \r\n                     \"\\n\", length (KOsets [[m]]), \" KOs in module vs all remaining \", (length (KO_cor_HOMA) - length (KOsets[[m]])), \" KOs\", \"\\n\")) +\r\n    xlab (\"SCC for KOs and HOMA-IR\") +\r\n    theme_classic () + theme (panel.grid.major = element_blank (), panel.grid.minor = element_blank ())\r\n  \r\n  ### Plot Distribution of correlations when leave-one-MGS-out.   \r\n  if (nrow (DeltaSCCperMGS [[m]]) > 1) {\r\n    g2 <- ggplot (DeltaSCCperMGS [[m]], aes (x = SCC_omiting_MGS)) + \r\n      geom_density (alpha = 1, fill = \"grey\", aes (y = ..scaled..)) + \r\n      geom_segment (aes (y = -0.1, yend = -0.02, x = SCC_omiting_MGS, xend = SCC_omiting_MGS)) +\r\n      geom_vline (data = cdf, aes (xintercept = rating.median, colour = cat), linetype = \"dashed\", size = 0.8) +\r\n      ggtitle (paste0 (\"KEGG module: \", m, \r\n                       \"\\n\", strsplit (module_mapping [m], split = \"\\\\[|,\")[[1]][1], \r\n                       \"\\nSCC for leave-1-MGS-out\",\r\n                       \"\\nNumber of MGSs = \", nrow (DeltaSCCperMGS [[m]]))) +\r\n      xlab (\"SCC for KOs and HOMA-IR\") +\r\n      theme_classic () + theme (panel.grid.major = element_blank (), panel.grid.minor = element_blank ())\r\n    \r\n  } else {\r\n    \r\n    ### Special circumstance when there is only one MGS (then one cannot make a \r\n    ### density plot, instead a black line is plotted at the value when that \r\n    ### one MGS is left out).\r\n    \r\n    g2 <- ggplot () +\r\n      scale_x_continuous (limits = range (c (DeltaSCCperMGS [[m]][, \"SCC_omiting_MGS\"], cdf$rating.median))) +\r\n      scale_y_continuous (name = \"\", limits = c (0, 1)) +\r\n      geom_vline (data = cdf, aes (xintercept = rating.median, colour = cat), linetype = \"dashed\", size = 0.8) +\r\n      geom_vline (data = DeltaSCCperMGS [[m]], aes (xintercept = SCC_omiting_MGS), linetype = \"longdash\", color = \"black\",size = 1) + \r\n      ggtitle (paste0 (\"KEGG module: \", m, \r\n                       \"\\n\", strsplit (module_mapping [m], split = \"\\\\[|,\") [[1]][1], \r\n                       \"\\nSCC for leave-1-MGS-out\",\r\n                       \"\\nNumber of MGSs = \", nrow (DeltaSCCperMGS [[m]]))) +\r\n      xlab (\"SCC for KOs and HOMA-IR\") +\r\n      theme_classic () + theme (panel.grid.major = element_blank (), panel.grid.minor = element_blank ()) \r\n    \r\n  }\r\n  \r\n  ### Plot distribution of correlations when leave-one-MGS-out. Bg.adjust SCC \r\n  ### (i.e. what we show in Figure 3c+d)\r\n  tmp.DeltaSCCperMGS = DeltaSCCperMGS [[m]]\r\n  ### Calculate delta-SCC in relation to the background adjusted SCC. \r\n  ### i.e the 'plottedvalue.s.m' shown as the second last equation in the \r\n  ### methods section in Pedersen et al., 2016.\r\n  tmp.DeltaSCCperMGS$DeltaMGS_SCC.bgadj = tmp.DeltaSCCperMGS$SCC.bgadj - tmp.DeltaSCCperMGS$DeltaMGS_SCC\r\n  ### specify the range of the x-axis to include 0, i.e. [min:0] for negative \r\n  ### correlations and [0:max] for positive correlations\r\n  x.range = c (min (0, tmp.DeltaSCCperMGS$DeltaMGS_SCC.bgadj), max (0, tmp.DeltaSCCperMGS$DeltaMGS_SCC.bgadj))\r\n  \r\n  if (nrow (tmp.DeltaSCCperMGS) > 1 ) {\r\n    \r\n    g3 <- ggplot (tmp.DeltaSCCperMGS, aes (x = DeltaMGS_SCC.bgadj)) + \r\n      geom_density (alpha = 1, fill = \"grey\", aes (y = ..scaled..)) + \r\n      geom_segment (aes (y = -0.1, yend = -0.02, x = DeltaMGS_SCC.bgadj, xend = DeltaMGS_SCC.bgadj)) +\r\n      ggtitle (paste0 (\"KEGG module: \", m, \r\n                       \"\\n\", strsplit (module_mapping [m], split = \"\\\\[|,\")[[1]][1], \r\n                       \"\\nbg.adj.SCC for leave-1-MGS-out\",\r\n                       \"\\nNumber of MGSs = \", nrow (tmp.DeltaSCCperMGS))) +\r\n      xlab (\"bg.adj.SCC for KOs and HOMA-IR\") + xlim (x.range) +\r\n      theme_classic () + theme (panel.grid.major = element_blank (), panel.grid.minor = element_blank ())\r\n    \r\n  } else { \r\n    \r\n    ### Special circumstance when there is only one MGS (then one cannot make a \r\n    ### density plot)\r\n    \r\n    g3 <- ggplot () +\r\n      ggtitle (paste0 (m, \r\n                       \"\\n\", strsplit (module_mapping [m], split=\"\\\\[|,\")[[1]][1], \r\n                       \"\\nbg.adj.SCC for leave-1-MGS-out\",\r\n                       \"\\nNumber of MGSs = \", nrow (tmp.DeltaSCCperMGS), \" - consequently not showing a density\/rug plot\")) \r\n    \r\n  }\r\n  \r\n  g = plot_grid (g1, g2, g3, ncol = 1, nrow = 3, align = \"v\", labels = c (\"a\", \"b\", \"c\"), axis = \"rl\")\r\n  # g = plot_grid (g1, g3, ncol = 1, nrow = 2, align = \"v\", labels = c (\"a\", \"b\", \"c\"), axis = \"rl\") ### use this if deleting g2\r\n  plot (g)\r\n  \r\n}\r\n\r\ndev.off ()  \r\n\r\n### END\r\n\r\n#print (sessionInfo ())\r\n\r\n# R version 3.3.3 (2017-03-06)\r\n# Platform: x86_64-apple-darwin13.4.0 (64-bit)\r\n# Running under: macOS Sierra 10.12.6\r\n# \r\n# locale:\r\n# [1] da_DK.UTF-8\/da_DK.UTF-8\/da_DK.UTF-8\/C\/da_DK.UTF-8\/da_DK.UTF-8\r\n# \r\n# attached base packages:\r\n# [1] stats     graphics  grDevices utils     datasets  methods   base     \r\n# \r\n# other attached packages:\r\n# [1] plyr_1.8.4            cowplot_0.9.1         ggplot2_2.2.1         gplots_3.0.1          ppcor_1.1            \r\n# [6] MASS_7.3-47           flashClust_1.01-2     WGCNA_1.61            fastcluster_1.1.24    dynamicTreeCut_1.63-1\r\n# [11] data.table_1.10.4-3  xlsx_0.5.7            xlsxjars_0.6.1        rJava_0.9-9          \r\n# \r\n# loaded via a namespace (and not attached):\r\n# [1]  Biobase_2.34.0        bit64_0.9-7           splines_3.3.3         foreach_1.4.3         gtools_3.5.0         \r\n# [6]  Formula_1.2-2         stats4_3.3.3          latticeExtra_0.6-28   blob_1.1.0            fit.models_0.5-14    \r\n# [11] yaml_2.1.14           robustbase_0.92-8     impute_1.48.0         RSQLite_2.0           backports_1.1.1      \r\n# [16] lattice_0.20-35       digest_0.6.12         RColorBrewer_1.1-2    checkmate_1.8.5       colorspace_1.3-2     \r\n# [21] htmltools_0.3.6       preprocessCore_1.36.0 Matrix_1.2-12         pcaPP_1.9-72          pkgconfig_2.0.1      \r\n# [26] GO.db_3.4.0           mvtnorm_1.0-6         scales_0.5.0          gdata_2.18.0          htmlTable_1.11.2     \r\n# [31] tibble_1.3.4          IRanges_2.8.2         nnet_7.3-12           BiocGenerics_0.20.0   lazyeval_0.2.1       \r\n# [36] survival_2.41-3       magrittr_1.5          memoise_1.1.0         doParallel_1.0.11     foreign_0.8-69       \r\n# [41] tools_3.3.3           matrixStats_0.52.2    stringr_1.2.0         S4Vectors_0.12.2      munsell_0.4.3        \r\n# [46] cluster_2.0.6         AnnotationDbi_1.36.2  caTools_1.17.1        rlang_0.1.4           grid_3.3.3           \r\n# [51] iterators_1.0.8       rstudioapi_0.7        htmlwidgets_0.9       robust_0.4-18         labeling_0.3         \r\n# [56] bitops_1.0-6          base64enc_0.1-3       gtable_0.2.0          codetools_0.2-15      DBI_0.7              \r\n# [61] reshape2_1.4.2        rrcov_1.4-3           gridExtra_2.3         knitr_1.17            bit_1.1-12           \r\n# [66] Hmisc_4.1-1           KernSmooth_2.23-15    stringi_1.1.6         parallel_3.3.3        Rcpp_0.12.13         \r\n# [71] rpart_4.1-11          acepack_1.4.1         DEoptimR_1.0-8       \r\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n","protected":false},"excerpt":{"rendered":"","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-442","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/442","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=442"}],"version-history":[{"count":3,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/442\/revisions"}],"predecessor-version":[{"id":449,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/442\/revisions\/449"}],"wp:attachment":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=442"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=442"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=442"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}