Integrative holo-omic data analysis predicts interactions across the 2 host-microbiome axis

Page 1 of 22 Supplementary material 1 Integrative holo-omic data analysis predicts interactions across the 2 host-microbiome axis 3 Merkesvik J a*, Langa JE bc, Pietroni C c, Alberdi A c, Poulsen LL d, Bojesen AM d, Meuronen T 4 e, Turunen S e, Kärkkäinen O ef, Westereng B a, Pope PB agh, Hvidsten TR a. 5 a Faculty of Chemistry, Biotechnology, and Food Science; Norwegian University of Life Sciences; Ås, Norway. 6 b Faculty of Science and Technology; University of the Basque Country; Leioa, Spain. 7 c Centre for Evolutionary Hologenomics; University of Copenhagen; Copenhagen, Denmark. 8 d Department of Veterinary and Animal Science, University of Copenhagen; Frederiksberg C, Denmark. 9 e Afekta Technologies Ltd; Kuopio, Finland. 10 f School of Pharmacy; University of Eastern Finland; Kuopio, Finland. 11 g Faculty of Biosciences; Norwegian University of Life Sciences; Ås, Norway. 12 h Centre for Microbiome Research; Queensland University of Technology; Woolloongabba, Australia. 13 * Corresponding author, email: jenny.merkes[email protected]. 14 Contents 15 Suppl. 1: Porcine solid feed composition ......................................................... 2 16 Suppl. 2: Method descriptions ....................................................................... 3 17 Suppl. 3: Additional visualisations ................................................................. 7 18 S3.1: Trial metadata on animal performance and feed administration 19 S3.2: Phylogenetic tree of MAG catalogue 20 S3.3: Differentially abundant MAGs across diet groups 21 S3.4: Differentially abundant MAGs across developmental stages 22 S3.5: Differentially abundant metabolites across developmental stages 23 S3.6: Differentially expressed host genes across diet groups 24 S3.7: Differentially expressed host genes across developmental stages 25 References ............................................................................................... 21 26 Abbreviations 27 AcGGM acetylated galactoglucomannan 28 FDR false discovery rate 29 HPLC high-performance liquid chromatography 30 LC-MS liquid chromatography mass spectrometry 31 LFC log2 fold change 32 MAG metagenome-assembled genome 33 MCFA multiset correlation and factor analysis 34 MS/MS tandem mass spectrometry 35 PCR polymerase chain reaction 36 rD-ratio robust dispersion ratio 37 rRSD robust relative standard devioation 38 UHPLC ultrahigh-performance liquid chromatography 39 Page 2 of 22 Suppl. 1: Porcine solid feed composition 40 Composition of the AgroSoft solid feed given to all animals during the trial. 41 nutrient unit per kg per energy nutrient unit per kg per energy FEsvin ny FEsv 1.219 1.0 potassium g 7.570 6.210 FE so ny FEso 1.199 0.983 chloride g 6.370 5.220 crude protein % 18.280 14.990 vitamin A, additive 1000 i.e 16.0 13.120 nitrogen % 2.930 2.40 vitamin D3, additive 1000 i.e 2.0 1.640 lignin % 2.870 2.350 vitamin E (synthetic), additive i.e 100.0 82.0 crude fat % 6.650 5.450 vitamin E/DL alpha tocopherol, additive mg 91.0 74.620 crude ash % 5.310 4.350 vitamin E eq. (from Proviox) mg 100.0 82.0 dry matter % 88.670 88.670 vitamin K3, additive mg 2.0 1.640 starch g 367.020 300.960 vitamin B1/thiamine, additive mg 10.0 8.20 sugar g 63.40 51.990 vitamin B2/riboflavin, additive mg 5.0 4.10 st.dig. crude protein g 166.630 136.640 vitamin B6/pyridoxin, additive mg 5.0 4.10 st.dig. lysine g 13.170 10.80 vitamin B12, additive mg 0.050 0.040 st.dig. methionine g 5.160 4.230 D-pantothenic acid, additive mg 35.0 28.70 st.dig. methionine+cystine g 7.770 6.370 niacin, additive mg 25.0 20.50 st.dig. threonine g 9.030 7.410 biotin vitamin H, additive mg 0.30 0.250 st.dig. tryptophan g 2.770 2.270 choline chloride, additive mg 350.0 287.0 st.dig. isoleucine g 6.530 5.350 folinic acid, additive mg 3.0 2.460 st.dig. valine g 8.670 7.110 L-carnitine mg 23.0 18.860 lysine g 13.970 11.450 betaine hydrochloride, additive mg 100.19 82.160 methionine g 5.380 4.410 Fe, iron II sulphate mg 120.0 98.40 threonine g 9.820 8.050 Cu, copper II sulphate mg 130.0 106.60 tryptophan g 2.990 2.450 Mn, manganese oxide mg 70.0 57.40 calcium g 6.460 5.30 Zn, zink sulphate mg 85.0 69.70 calcium formate mg 0.0 0.0 I, calcium iodate mg 3.06 2.510 phosphorus g 6.720 5.510 Se, sodium selenite mg 0.25 0.210 dig. phospgorus, 0% phytase g 3.820 3.130 6-phytase (EC 3.1.3.26) OTU 265.0 217.30 dig. phospgorus, 60% phytase g 4.080 3.350 endo-1,4-beta-xylanase EPU 1500.0 1230.020 dig. phospgorus, 100% phytase g 4.190 3.430 beta glucanase (EC 3.2.1.5) 1000U 0.10 0.080 dig. phospgorus, 150% phytase g 2.280 3.510 BHT (antioxidant) mg 10.0 8.20 dig. phospgorus, 200% phytase g 4.350 3.560 antioxidant mg 21.52 17.650 sodium g 3.060 2.510 propyle gallate mg 5.0 4.10 magnesium g 1.450 1.190 ProHacidⓇ mg 7000.0 5740.080 magnesium sulphate g 0.0 0.0 42 Page 3 of 22 Suppl. 2: Method descriptions 43 Acetylated galactoglucomannan fibres 44 AcGGM was produced following the protocol described in Michalak et al. [1]. In 45 summary, dry wood from Norway spruce (Picea abies) was milled into chips and 46 soaked in a solution with sodium citrate and potassium phosphate buffers. The wood 47 chips were steam exploded at 200°C and 14.5 bar in a 20 L pressure vessel, using a 48 25 kW electric boiler (Parat, Norway) for steam production. Water was then added to 49 create a slurry, which was filtered to separate the water-insoluble material from the 50 wood chips. The solid mass was dried at 100°C for 48 hours and cut into smaller 51 pieces using a cutting mill (Retsch, Germany) with a 0.5 mm sieve. 52 Animal trial and sampling 53 The feeding trial was conducted at the University of Copenhagen during spring 2022 54 under the experimental animal license number 2020-15-0201-00520. Twelve ten-day 55 old male piglets (2.8-4.9 kg) of the same litter from Bøgholt farm (Blangstrupsvej 17, 56 DK-5610 Assens) were housed individually and divided into three weight-balanced 57 experimental groups with specific feeding regimens over a period of six weeks. All 58 received the same basal diet: four weeks of DanMilk Supreme 1.0 milk replacement 59 formula (AB Neo, Denmark), and two weeks of solid feed (AgroSoft, Denmark) (Suppl. 60 1). All piglets were weaned at age 24 days. Two of three groups received AcGGM 61 fibres as a supplement at 4% inclusion, but starting at different time points: one for 62 the final four weeks of the trial, spanning two weeks before and two weeks after 63 weaning; the second only for the final two trial weeks, starting after weaning. Body 64 weight, feed intake, and clinical scores are visualised in Fig. S3.1. After the trial’s 65 sixth week, the piglets were euthanized by Zoletil anaesthesia (0.1 mL/kg) followed 66 by a lethal intracardial pentobarbital dose (0.25 mL/kg). 67 Sample collection and preparation 68 Endpoint samples were collected from the caecum of each animal. Digesta (100 mg) 69 was sampled in all 12 animals for metagenomic, metatranscriptomic, and untargeted 70 metabolomics analyses, whereas tissue samples (1x1 cm) were collected from nine 71 individuals for transcriptomic and untargeted metabolomics analyses. Samples 72 intended for metagenomics and (meta)transcriptomics were stabilised in DNA/RNA 73 Shield (1 mL) (Zymo, USA) to prevent nucleic acid degradation and enable long-term 74 storage at -20°C until further processing for data generation. For untargeted 75 metabolomics, digesta and tissue samples were snap frozen in liquid nitrogen and 76 stored at -80°C until further processing. 77 Samples for metagenomics, metatranscriptomics, and host transcriptomics were 78 prepared following the laboratory workflow of the Earth Hologenome Initiative 79 (https://www.earthhologenome.org/laboratory). Briefly, thawed samples were 80 mechanically lysed using a Lysing Matrix E (MP Biomedicals, USA) with two 6-minute 81 bead-beating cycles at 30 Hz on a TissueLyser II (Qiagen, Germany). Silica magnetic 82 Page 4 of 22 beads and solid-phase reversible immobilisation were used to isolate RNA and DNA in 83 separate fractions and to remove inhibitors. DNA and RNA concentrations were 84 quantified with a Quibit 3 Fluorometer using Quibit HS or BR Assay Kits (Thermo 85 Fisher Scientific, USA). For metagenomics, DNA was fragmented to ~400 bp by 86 ultrasonication using a Covaris LE220 platform. Library preparation was conducted 87 following a blunt-end adapter ligation single-tube protocol [2] with optimised 88 conditions as described by Mak et al. [3]. qPCR screening was used to determine the 89 number of PCR cycles needed for library indexing. Samples were then amplified using 90 unique dual index primers before the libraries were analysed through capillary 91 electrophoresis with Fragment Analyzer (Agilent, USA). Metagenomic libraries were 92 generated using the Illumina Stranded Total RNA Prep kit (Illumina, USA) with 93 enzymatic rRNA depletion via Ribo-Zero Plus (Illumina, USA). Transcriptomic libraries 94 employed bead-based rRNA depletion using the TruSeq Stranded Total RNA Library 95 Prep kit for Human/Mouse/Rat (Illumina, USA). Sequencing libraries were pooled and 96 sequenced with an Illumina NovaSeq6000 system using S4 flow cells, generating 97 approx. 10 GB per (meta)transcriptomic library, and 5 GB per metagenomic library. 98 For untargeted metabolomics, digesta and tissue samples were transferred to cooled 99 Bead Ruptor homogeniser tubes (Omni International, USA) with four metal beads and 100 cold methanol (80%, 400 µL/100 mg tissue, 1000 µL/100 g digesta). Digesta samples 101 were homogenised using a Bead Ruptor 24 Elite (Omni International, USA) (2±2°C, 102 7.0 m/s, 2x30 seconds, 20 seconds between cycles). Tissue samples followed the 103 same homogenisation protocol repeated once. Samples were incubated on ice for 15 104 minutes before being vortexed (10 seconds) and centrifuged (4°C, 17,000 G, 10 105 minutes). Supernatants were decanted through a 0.2 µm Acrodisc filter with 106 polytetrafluoroethylene membrane (Pall, USA) into HPLC vials (Agilent, USA) for LC107 MS analysis. Quality control samples – one for digesta and one for tissue samples – 108 were prepared by pooling 150 µL of supernatants from each sample of each respective 109 sample type into another HPLC vial. The two pooled samples were vortexed and 110 filtered like the other samples. LC-MS analysis utilised a 1290 Infinity II UHPLC 111 (Agilent, USA) coupled with a high-resolution quadrupole time of flight Agilent 6546 112 mass spectrometer with a Dual Jet Stream ion source (Agilent, USA) and followed the 113 protocol described in Hanheniva et al. [4] and Klåvus et al. [5]. In brief, a Zorbax 114 Eclipse XDB-C18 column (2.1×100 mm, 1.8 µm) (Agilent, USA) was used for 115 reversed-phase separation, while hydrophilic interaction liquid chromatography 116 separation was conducted with an UHPLC ethylene bridged hybrid amide column 117 (Waters, USA). After each chromatographic run, samples were ionised through Jet 118 Stream electrospray in positive and negative mode, yielding four data files per sample. 119 Collision energies for MS/MS analysis were selected as 10, 20, and 40 V for 120 compatibility with spectral databases. 121 Omic data generation 122 Page 5 of 22 Metagenomics, metatranscriptomics, and host transcriptomics sequencing data were 123 processed using Snakemake [6] workflows developed within the 3D’omics project, 124 which are available on GitHub (https://github.com/3d-omics). 125 For untargeted metabolite data generation, peak detection and alignment was 126 performed in MS-DIAL (ver. 4.92) [7]. Mass-to-charge ratios between 50 and 1,500 127 at any retention time were considered for peak collection. Minimum peak height 128 amplitude was set to 2,000, and detection used the linear weighted moving average 129 algorithm. To align peaks across samples, tolerances were set to 0.2 min retention 130 time and 0.015 Da m/z. Solvent background was removed using solvent blank 131 samples with the condition that the maximum signal abundance across samples was 132 at least five times higher than the mean signal in the blanks. A total of 62,526 133 metabolite features remained after peak picking, and were subjected to further 134 preprocessing and clean-up with the R-package notame (ver. 0.2.0) [5] for each 135 sample type separately. Metabolomic features were considered high quality based on 136 presence (>70% of quality control samples and >50% of samples in at least one 137 study group); rRSD (<20%); and rD-ratio (<10%). In addition, features with too high 138 rRSD or rD-ratio were still considered good quality if their classic RSD, rRSD, and 139 basic D-ratio were all low (<10%). Low-quality features were flagged and discarded, 140 leaving 16,328 metabolite features in the dataset. 141 Compound identification was conducted with MS-DIAL by comparing chromatographic 142 and spectrometric characteristics with both in-house and publicly available databases 143 (MassBank, ReSpect, RIKEN, GNPS, Fiehn libraries, CASMI2016, MetaboBASE, 144 BMDMS-NP, and PFAS). Additional identifiers were obtained using MS-FINDER (ver. 145 3.60) [8] to compare acquired MS/MS spectra with in silico-generated records for 146 known compounds. All metabolite annotations followed the Metabolomics Standards 147 Initiative by the Chemical Analysis Working Group [9]. 148 The metagenome count matrix was filtered for genome coverage (>50% in at least 149 two samples) and presence in more than one sample, reducing the MAG catalogue 150 from 439 to 373 populations, spanning 815,027 microbial genes. The host genome 151 saw 27,376 mapped genes within these samples, of which 22,170 were expressed in 152 at least two samples and thus included in the analysis. From the 16,284 high-quality 153 metabolic features, 6,388 with acquired MS/MS data were kept in the metabolomic 154 dataset before initiating omic data analyses. 155 Omic data analyses 156 All statistical analyses were conducted in R (ver. 4.4.2) [10] on an aarch64 Apple 157 Darwin20 platform running macOS Sequoia 15.7, unless stated otherwise. Relevant 158 data and metadata files are available in the GitHub repository 159 jennymerkesvik/3domics_wp3-2, which also contains an RMarkdown report 160 demonstrating the full analysis approach. 161 Page 6 of 22 Read count-based omic datasets (metagenomics, metatranscriptomics, and host 162 transcriptomics) were compared across sample groups with DESeq2 (ver. 1.46.0) 163 [11], using default settings and significance thresholds of absolute LFC > 1; FDR p164 value < 0.05; and base mean > 50. Molecular features from metabolomics were 165 compared across pairs of sample groups using two-tailed t-tests assuming equal 166 variance, false discovery rate-corrected p-values, and log2 fold changes. Variance167 stabilising transformations were used to normalise the read count data before they 168 were used for holo-omic integration. Metabolomic features were log2-transformed for 169 the same purpose. R packages ggplot2 (ver. 3.5.1) [12], ggtree (ver. 3.14.0) [13], 170 and ggtreeExtra (ver. 1.16.0) [14] were used to visualise analysis results. 171 Using Python (ver. 3.12.9) [15], we used MCFA (ver. 1.0.2) [16] to model the holo172 omic dataset. The utilised script is found on GitHub, in which standardised omic data 173 layers across common samples were used as input for the MCFA function fit(), while 174 querying a one-dimensional shared space (d=1). Shared and private components and 175 weights, along with files informing on variance explained by each model space, was 176 imported into R to inspect and visualise model performance. Features with weight for 177 the shared model component in the 93rd percentile were selected as the most relevant 178 features for identifying host-microbiome interactions, while yielding a manageable 179 number of features for predicting cross-omic relationships. A co-occurrence network 180 of these omic data were created with FlashWeave (ver. 0.19.2) [17], which was 181 implemented in Julia (ver. 1.11.6) [18]. Specifically, the function learn_network() 182 was provided a matrix consisting of the selected omic data and one-hot encoded 183 metadata on diet group and developmental stage for dietary fibre introduction; a 184 mask to inform on which features were experimental covariates; and additional 185 settings sensitive = true, heterogenous = false, n_obs_min = 7, normalize 186 = false. The learned network was imported into Cytoscape (ver. 3.10.3) [19] and R 187 for inspection and visualisation. 188 Page 7 of 22 Suppl. 3: Additional visualisations 189 Metadata on animal performance and feed administration (Fig. S3.1). MAG catalogue 190 composition and differential transcriptomic activity (Fig. S3.2). Differentially 191 abundant features across metagenomic (Figs. S3.3-2), metabolomic (Fig. S3.5), 192 and host transcriptomic (Figs. S3.6-5) data layers. 193 194 Figure S3.1. Animal performance and administered feed throughout the feeding trial. 195 Page 8 of 22 196 Figure S3.2. Phylogenetic tree of the 370 bacterial populations in the MAG catalogue. The most specific 197 taxonomic classification is used as branch labels. Repeated taxon identifiers are numerated by order of 198 appearance. A) Sum of metatranscriptomic reads mapped to each MAG in animals fed AcGGM, split into 199 fibre introduction with or without AcGGM (outer ring) and before and after weaning (inner ring). 200 Significant differences (LFC>1, FDR p-value<0.05) in activity levels are marked with asterisks below 201 each bar, with colour indicating which of diet group had the most transcriptionally active population of 202 the MAG. B) Genome sizes (bar height) and completeness (colour) of MAGs. C) Taxonomic classification 203 of MAGs on the order level. 204 Page 9 of 22 The remainer of Suppl. 3 contains overviews of differentially abundant features from 205 metagenomic, metabolomic, and host transcriptomic layers. Due to the high number 206 of differentially expressed microbial genes, metatranscriptomics is not included. Lists 207 of these features are instead found on GitHub (jennymerkesvik/3domics_wp3-2). 208 209 Figure S3.3. Differentially abundant metagenome-assembled genomes in microbiomes of animals given 210 different diets. Significance is indicated with horizontal bars, reporting log2 fold change and false 211 discovery rate-adjusted p-values (* < 0.05, ** < 0.01, *** < 0.001). Continued on next page. 212 Page 16 of 22 228 Figure S3.5. Differentially abundant metabolites in digesta or tissue samples of animals introduced to 229 acetylated galactoglucomannan fibres at different developmental stages. Significance is indicated with 230 horizontal bars, reporting log2 fold change and false discovery rate-adjusted p-values (* < 0.05, ** < 231 0.01, *** < 0.001). 232 233 234 Figure S3.6. Differentially expressed genes in caecal gut wall samples of animals given different diets. 235 Significance is indicated with horizontal bars, reporting log2 fold change and false discovery rate236 adjusted p-values (* < 0.05, ** < 0.01, *** < 0.001). 237 Page 17 of 22 238 Figure S3.7. Differentially expressed genes in caecal gut wall samples of animals introduced to 239 acetylated galactoglucomannan fibres starting at different developmental stages. Significance is 240 indicated with horizontal bars, reporting log2 fold change and false discovery rate-adjusted p-values (* 241 < 0.05, ** < 0.01, *** < 0.001). Continued on next page. 242 Page 18 of 22 243 Figure S3.7 continued. 244 Page 19 of 22 245 Figure S3.7 continued. 246 Page 20 of 22 247 Figure S3.7 continued (end). 248 Page 21 of 22 References 249 1. Michalak L, Knutsen SH, Aarum I, Westereng B. Effects of pH on steam explosion 250 extraction of acetylated galactoglucomannan from Norway spruce. Biotechnol Biofuels. 251 2018;11(1):311. doi: 10.1186/s13068-018-1300-z. 252 2. Carøe C, Gopalakrishnan S, Vinner L, Mak SST, Sinding MHS, Samaniego JA et al. 253 Single-tube library preparation for degraded DNA. Methods Ecol Evol. 2018;9(2):410– 254 9. doi: 10.1111/2041-210x.12871. 255 3. Mak SST, Gopalakrishnan S, Carøe C, Geng C, Liu S, Sinding M-HS et al. 256 Comparative performance of the BGISEQ-500 vs Illumina HiSeq2500 sequencing 257 platforms for palaeogenomic sequencing. GigaScience. 2017;6(8):gix049. doi: 258 10.1093/gigascience/gix049. 259 4. Hanhineva K, Lankinen MA, Pedret A, Schwab U, Kolehmainen M, Paananen J et al. 260 Nontargeted Metabolite Profiling Discriminates Diet-Specific Biomarkers for 261 Consumption of Whole Grains, Fatty Fish, and Bilberries in a Randomized Controlled 262 Trial. J Nutr. 2015;145(1):7–17. doi: 10.3945/jn.114.196840. 263 5. Klåvus A, Kokla M, Noerman S, Koistinen VM, Tuomainen M, Zarei I et al. “Notame”: 264 Workflow for Non-Targeted LC–MS Metabolic Profiling. Metabolites. 2020;10(4):135. 265 doi: 10.3390/metabo10040135. 266 6. Köster J, Rahmann S. Snakemake—a scalable bioinformatics workflow engine 267 [software]. Bioinformatics. 2012;28(19):2520–2. doi: 268 10.1093/bioinformatics/bts480. 269 7. Tsugawa H, Cajka T, Kind T, Ma Y, Higgins B, Ikeda K et al. MS-DIAL: data270 independent MS/MS deconvolution for comprehensive metabolome analysis 271 [software]. Nat Methods. 2015;12(6):523–6. doi: 10.1038/nmeth.3393. 272 8. Tsugawa H, Kind T, Nakabayashi R, Yukihira D, Tanaka W, Cajka T et al. Hydrogen 273 Rearrangement Rules: Computational MS/MS Fragmentation and Structure 274 Elucidation Using MS-FINDER Software [software]. Anal Chem. 2016;88(16):7946– 275 58. doi: 10.1021/acs.analchem.6b00770. 276 9. Sumner LW, Amberg A, Barrett D, Beale MH, Beger R, Daykin CA et al. Proposed 277 minimum reporting standards for chemical analysis: Chemical Analysis Working 278 Group (CAWG) Metabolomics Standards Initiative (MSI). Metabolomics. 279 2007;3(3):211–21. doi: 10.1007/s11306-007-0082-2. 280 10. R Core Team. R: A language and environment for statistical computing [software]. 281 2022. https://www.R-project.org/. 282 11. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion 283 for RNA-seq data with DESeq2 [software]. Genome Biol. 2014;15(12):550. doi: 284 10.1186/s13059-014-0550-8. 285 Page 22 of 22 12. Wickham H. ggplot2 [software]. Cham: Springer International Publishing. 286 2016:2nd ed. doi: 10.1007/978-3-319-24277-4. 287 13. Yu G, Smith DK, Zhu H, Guan Y, Lam TT. ggtree: an R package for visualization 288 and annotation of phylogenetic trees with their covariates and other associated data 289 [software]. Methods Ecol Evol. 2017;8(1):28–36. doi: 10.1111/2041-210X.12628. 290 14. Xu S, Dai Z, Guo P, Fu X, Liu S, Zhou L et al. ggtreeExtra: Compact Visualization 291 of Richly Annotated Phylogenetic Data [software]. Mol Biol Evol. 2021;38(9):4039– 292 42. doi: 10.1093/molbev/msab166. 293 15. Guido van R, Jelke de B. Interactively testing remote servers using the Python 294 programming language [software]. CWI Quarterly. 1991;4(4):283–304. 295 16. Brown BC, Wang C, Kasela S, Aguet F, Nachun DC, Taylor KD et al. Multiset 296 correlation and factor analysis enables exploration of multi-omics data. Cell Genomics. 297 2023:100359. doi: 10.1016/j.xgen.2023.100359. 298 17. Tackmann J, Matias Rodrigues JF, Von Mering C. Rapid Inference of Direct 299 Interactions in Large-Scale Ecological Networks from Heterogeneous Microbial 300 Sequencing Data [software]. Cell Systems. 2019;9(3):286-296.e8. doi: 301 10.1016/j.cels.2019.08.002. 302 18. Bezanson J, Edelman A, Karpinski S, Shah VB. Julia: A Fresh Approach to 303 Numerical Computing [software]. SIAM Rev. 2017;59(1):65–98. doi: 304 10.1137/141000671. 305 19. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D et al. Cytoscape: 306 A Software Environment for Integrated Models of Biomolecular Interaction Networks 307 [software]. Genome Res. 2003;13(11):2498–504. doi: 10.1101/gr.1239303. 308