MetaMap: An atlas of metatranscriptomic reads in human disease-related RNA-seq data

GigaScience, 7, 2018, 1–8 doi: 10.1093/gigascience/giy070 Advance Access Publication Date: 12 June 2018 Data Note DATA NOTE MetaMap: an atlas of metatranscriptomic reads in human disease-related RNA-seq data L.M. Simon 1,*,S.Karg 1, A.J. Westermann2,3,M.Engel 1,4, A.H.A. Elbehery5, B. Hense1,M.Heinig 1,L.Deng 5and F.J. Theis 1,6,* 1Helmholtz Zentrum M ¨ unchen, German Research Center for Environmental Health, Institute of Computational Biology, Neuherberg, Germany, 2Institute of Molecular Infection Biology, University of W¨ urzburg, W ¨ urzburg, Germany, 3Helmholtz Institute for RNA-Based Infection Research, W¨ urzburg, Germany, 4Helmholtz Zentrum M ¨ unchen, German Research Center for Environmental Health, Scientific Computing Research Unit, Neuherberg, Germany, 5Helmholtz Zentrum M ¨ unchen, German Research Center for Environmental Health, Institute of Virology, Neuherberg, Germany and 6Department of Mathematics, Technische Universit¨ at M ¨ unchen, Munich, Germany ∗Correspondence address. L.M. Simon, Helmholtz Zentrum M ¨ unchen German Research Center for Environmental Health, Institute of Computational Biology, Ingolst¨ adter Landstraße, 185764, Neuherberg; E-mail: [email protected] http://orcid.org/0000-0001-6148-8861; F.J. Theis; E-mail: [email protected] http://orcid.org/0000-0002-2419-1943 Abstract Background: With the advent of the age of big data in bioinformatics, large volumes of data and high-performance computing power enable researchers to perform re-analyses of publicly available datasets at an unprecedented scale. Ever more studies imply the microbiome in both normal human physiology and a wide range of diseases. RNA sequencing technology (RNA-seq) is commonly used to infer global eukaryotic gene expression patterns under defined conditions, including human disease-related contexts; however, its generic nature also enables the detection of microbial and viral transcripts. Findings: We developed a bioinformatic pipeline to screen existing human RNA-seq datasets for the presence of microbial and viral reads by re-inspecting the non-human-mapping read fraction. We validated this approach by recapitulating outcomes from six independent, controlled infection experiments of cell line models and compared them with an alternative metatranscriptomic mapping strategy. We then applied the pipeline to close to 150 terabytes of publicly available raw RNA-seq data from more than 17,000 samples from more than 400 studies relevant to human disease using state-of-the-art high-performance computing systems. The resulting data from this large-scale re-analysis are made available in the presented MetaMap resource. Conclusions: Our results demonstrate that common human RNA-seq data, including those archived in public repositories, might contain valuable information to correlate microbial and viral detection patterns with diverse diseases. The presented MetaMap database thus provides a rich resource for hypothesis generation toward the role of the microbiome in human disease. Additionally, codes to process new datasets and perform statistical analyses are made available. Keywords: high-performance computing; big data; RNA-seq; sequence read archive; metatranscriptomics; microbiome; virome; human disease; infection Received: 22 February 2018; Revised: 1 June 2018; Accepted: 4 June 2018 C The Author(s) 2018. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 1 Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 2MetaMap Data Description Context Recent studies have demonstrated the paramount importance of the microbiome for human health and disease [1]. For example, imbalance of the human gut microbiome was linked to noncommunicable diseases such as obesity [2,3], diabetes [4], cardiovascular disease [5], chronic obstructive pulmonary disease [6], and colorectal carcinoma [7,8], to name just a few. The advent of high-throughput sequencing technologies has revolutionized the life sciences. RNA sequencing (RNA-seq) technology produces one of the most frequent next-generation sequencing data types and has been applied to the study of a large number of biological samples relevant to human disease. The majority of the underlying raw data are freely accessible from data repositories such as the Gene Expression Omnibus (>1,700 human RNA-seq datasets as of January 2018) and the Sequence Read Archive (SRA) [9]. However, these data are typically exclusively used for single species (i.e., human) transcriptomics such as differential gene expression and alternative splicing analysis [9,10]. Reads that do not map onto the human genome are considered noise or contamination and therefore are generally ignored [11,12](collectively about 9% of total reads, Fig. 1). Five years ago, it was postulated that interspecies interactions might be studied by simultaneous detection and quantification of RNA transcripts from a given host and a microbe via “dual” RNA-seq [13]. Meanwhile, this approach has been successfully applied to the interaction of mammalian cells with diverse bacterial [14] and viral pathogens [15-19]. Inspired by dual RNA-seq, in this study we hypothesize that reads in archived RNA-seq datasets derived from human primary cells or tissue samples that fail to map against the human reference genome may contain valuable information about thepresenceofcertainmicrobesintherespectivebodyniches and/or under defined disease conditions. To enable metatranscriptomic study of these data, we combined existing read alignment and metagenomic classification software into a two-step “omni” RNA-seq pipeline to comprehensively quantify archaeal, bacterial, and viral reads in human RNA-seq data (Fig.1). In the first step of this so-called MetaMap pipeline, all reads are aligned against the human genome using the ultra-fast RNAseq aligner Spliced Transcripts Alignment to a Reference software (STAR) [20]. Subsequently, only the fraction of unmapped reads is subjected to metatranscriptomic classification using CLARK-S [21] (see Methods for details). The combination between scalability and accuracy was the main motivation behind choosing these two software packages over competing methods [22,23]. It is important to note that CLARK-S uses a set of uniquely discriminative short sequences at the species level to classify reads. Therefore, reads containing nondiscriminative sequences that fail to be uniquely assigned to a single species, e.g., reads originating from the bacterial ribosomal 16S rRNA gene will be considered “unclassified” (altogether 8.6% in Fig.1). The output of CLARK-S is an operational taxonomic units (OTUs) count matrix, where rows correspond to viral, bacterial, and archeal species and columns correspond to (human) samples. Each entry corresponds to the number of non-human reads classified to the respective species. For convenience, in the following, we refer to the set of microbial and viral species profiled using our approach as “metafeatures.” By screening the study abstracts of the SRA for search terms prioritizing human clinical datasets derived from polyAindependent sequencing protocols (see Methods), we identified more than 400 studies relevant to human disease comprising more than 17,000 cDNA libraries (close to 150 terabytes of raw sequencing data). Raw sequencing reads from these studies were downloaded and analyzed using the high-performance computing system of the Leibniz Supercomputing Centre (LRZ) of the Bavarian Academy of Sciences and Humanities, which facilitated ultra-fast processing with median speeds of 25 and 21 million reads per hour per core per run for the STAR and CLARKS steps, respectively. Of the more than 500 billion RNA-seq reads processed, around 91% could be mapped to the human genome. A fraction of 8.6% of all reads remained nondiscriminative at the species level and defined as “unclassified.” In addition, 0.03%, 0.20%, and 0.39% of all reads were assigned to archaeal, bacterial, or viral metafeatures, respectively. Despite these relatively low percentages, the absolute numbers of reads classified were in the hundred millions to billions, enabling statistical analyses. Methods High-performance computing environment Project computations including download, alignment of reads onto the human genome, and metafeature quantification were made on the high-performance Linux Cluster at the LRZ [24]. RNA-seq data retrieval Raw next-generation sequencing data were downloaded from the SRA. The R package SRAdb was downloaded on 23 May 2017 and used to query the SRA database. To identify SRA projects that contain transcriptomic analyses of human RNA-seq data, the SRA attributes “taxon id,” “library source,” “library strategy,” and “platform” were searched for the terms “9606,” “TRANSCRIPT,” “RNA-seq,” and “ILLUMINA,” respectively. To remove potential bias derived from different sequencing technologies, we also restricted the query to SRA runs annotated with “ILLUMINA” in SRA attribute “platform.” To exclude studies with insufficient sample size for statistical analysis, the query was restricted to SRA projects containing more than five runs. To avoid concentrating the analysis on a small number of large projects, the query was restricted to SRA projects with fewer than 500 runs. To identify studies focusing on phenotypes relevant to human disease, we restricted the query to runs containing at least one or more of the terms “disease,” “patient,” “primary,” and “clinical” in the SRA attribute “study abstract.” To exclude in vitro (cell-culture) experiments but focus on primary (clinical) samples, SRA runs containing the terms “mutant” or “cell-line” were removed from our selection. Furthermore, SRA runs containing the terms “single cell” and “GTEx” were removed. Finally, samples with fewer than 1 million total reads or read lengths <50 bp were excluded. The described query resulted in 484 short read projects (SRPs) containing 21,659 RNA-seq runs. Due to technical problems (i.e., missing URLs, restricted access), we were unable to download a fraction of 4,078 samples. Human alignment Alignment of reads against the human reference genome (hg38) and simultaneous human gene expression quantification was conducted with STAR (version 2.5.2). To increase mapping speed of a large number of samples, we used the –genomeLoad LoadAndKeep function to load the STAR index once and keep it in memory for subsequent alignments. The parameter –quantmode GeneCounts was used to generate the human gene expression count tables. Unmapped reads were saved with the – outReadsUnmapped Fastx parameter. To further increase mapping speed, multiple threads were used as implemented with the paDownloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 Simon et al. 3 Figure 1: Schematic of the MetaMap pipeline. More than 400 projects from studies relevant to human disease were identified in the SRA database. More than 500 billion RNA-seq reads were downloaded and first filtered by mapping them onto the human genome. The remaining reads underwent metafeature classification. It is noted that 90.7% of all reads mapped to the human genome; 0.03%, 0.20%, and 0.39% of all reads were assigned to archaeal, bacterial, or viral metafeatures, respectively; and 8.6% of all reads remain nondiscriminative at the species level (“unclassified”). rameter –runThreadN 28. Runs with fewer than 30% reads mapping to the human genome were excluded from downstream analysis. All human alignments were conducted on the LRZ “CoolMUC2” Linux-Cluster. This cluster contains 384 nodes with 64 GB random access memory (RAM) memory and 28 cores each. Metafeature quantification Metafeature quantification was conducted with CLARK-S (version 1.2.3). CLARK-S is a software method for fast and accurate sequence classification of metagenomic next-generation sequencing data, including RNA-seq data. One major issue during the classification of metagenomic data is the rising number of targets to align against. CLARK-S solves this issue by building a large index file consisting of discriminative k-mers. The metagenomic reference database was generated following the description of the CLARK website using the following two commands: set targets.sh bacteria virus –species and buildSpacedDB.sh. This database contained 16,551 genome sequences corresponding to 6,979 unique species (Additional File 2). To allow uniform processing, paired-end sequencing experiments were analyzed independently. Each single unmapped read file was used as input for CLARK-S with the following parameters: classify metagenome.sh –spaced –O list of FASTQ files. To increase classification speed, the CLARK-S express mode was selected and multiple threads were used with parameters –m 2 and –n 32, respectively. The output files of this step contain all input read identifiers with the corresponding metafeature classification. In the subsequent step, total counts are summarized for each feature with the estimate abundance.sh command. To enable comparison across single-end and paired-end experiments, metafeature counts from paired-end experiments were averaged and subsequently rounded to conserve count distribution. To account for varying sequencing depths, metafeature abundance was estimated as the number of reads per million total reads sequenced. Metafeature quantification was conducted on the LRZ “Teramem” Linux-Cluster. This cluster contains one node with 6,144 GB RAM memory and 96 cores. BLAST-based metafeature classification To validate results generated by the MetaMap pipeline, the Basic Local Alignment Search Tool (BLAST) [25] was used as follows. A BLAST database was created from the same genome sequences used in the CLARK-S approach. Then, reads were aligned to this database using BLASTN with a threshold E-value of 1e-10. Produced counts from paired-end experiments were averaged. For each file, BLAST was done by running approximately 10 kb chunks (record separator “>”) in parallel using parallel [41] (28 jobs), each with eight threads using one node on the LRZ “CoolMUC3” Linux Cluster. This cluster contains 148 nodes with 96 GB RAM memory and 64 cores each. Output was parsed to exclusively keep reads that could be assigned at the species level. Differential metafeature abundance Differential metafeature abundance analysis was performed using the R package DESeq2 [26]. DESeq2 models differential gene expression by fitting a negative binomial distribution to the raw counts underlying RNA-seq data. This framework can account for confounding variables such as sequencing depth. Therefore, the data need not be normalized prior to statistical inter-sample comparisons. For each of the four published bona fide dual RNAseq studies, we classified samples into the following two groups based on the provided annotations: samples expected to contain the known pathogen, such as human papillomavirus-positive tumors in the Zhang et al. study [28], and pathogen-free controls, such as mock-treated cells in the Westermann et al. [27] study. Using this binary outcome, we performed differential expression analysis across all detected metafeatures. To account for sequencing depth, library size factors were estimated from the total number of sequenced reads. The dispersion for the negative binomial distribution was estimated using a local linear regression as implemented in the DESeq() function via the fitType parameter “local.” Data validation and quality control We validated our approach by recovering the ground truth in bona fide dual RNA-seq experiments performed with human cell lines and samples from patients with well-known infection status. Of the four selected studies, one analyzed an infection model based on a bacterial (Salmonella enterica serovar Typhimurium) and three based on distinct viral pathogens (human papillomavirus, herpes simplex virus, rhinovirus). As expected, MetaMap detected the known pathogen at higher levels in the respective study compared to the other studies and pathogens (Table1). However, comparisons across studies and metafeatures may be biased by technical confounders (discussed in detail in the Re-use potential section). Therefore, we focused our analysis on the comparison of a single metafeature across subjects within a study. Using the annotation provided in the respective study, we performed differential metafeature abundance analysis to identify those metafeatures that show the largest relative difference in abundance levels between the infected and control Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 4MetaMap Table 1: Overview of four dual RNA-seq studies used to validate the MetaMap pipeline. Study Infection agent Total reads Salmonella enterica Alphapapillomavirus 9 Human alphaherpesvirus 1 Rhinovirus A Westermann et al. [27] Salmonella enterica serovar Typhimurium 1.0e+07 6.3e+03 1.2e-01 1.5e-01 1.2e-01 Zhang et al. [28] Human papillomavirus 4.6e+07 3.0e-02 5.1e+01 2.2e-02 2.2e-02 Rutkowski et al. [29] Herpes simplex virus 3.5e+07 1.1e+00 3.1e-02 3.1e+04 3.0e-02 Bai et al. [30] Rhinovirus 6.6e+06 2.0e-01 1.5e-01 1.5e-01 4.4e+01 Total reads column depicts the average read depth per sample for each study. Average metafeature abundance for alphapapillomavirus 9, Salmonella enterica, human alphaherpesvirus 1, and rhinovirus A are shown in reads per million. The correct infection agent for the respective study is highlighted in bold font 0 25 50 75 0510 −log10 p−value Westermann et al 0 5000 10000 15000 infected (n=36) mock (n=6) Salmonella enterica A 0 10 20 04812 Zhang et al Alphapapillomavirus 9 B 0.0 5.0 10.0 −5.0 −2.5 0.0 2.5 5.0 7.5 Rutkowski et al 0 25000 50000 75000 infected (n=8) mock (n=2) H. Alphaherpesvirus 1 C 0 20 40 048 Fold change (log2) Bai et al 0 50 100 150 200 infected (n=12) vehicle (n=12) Rhinovirus A D 0 100 200 300 400 500 negative (n=18) positive (n=18) Figure 2: Differential metafeature abundance analysis of controlled infection experiments recovers ground truth. “Volcano” plots show fold change and inverted Pvalue on the xand yaxes, respectively. Each dot represents a metafeature. The most significant metafeature is colored in red. Insets display box plots of the abundance levels in reads per million of the top hit metafeature across conditions for each study. For all box plots, the box represents the interquartile range, the horizontal line in the box is the median, and the whiskers represent 1.5 times the interquartile range. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 Simon et al. 5 Figure 3: Analysis of lymphoblast cell line experiments further supports the MetaMap pipeline. (A and B) Mean abundance levels across all samples of the top five metafeatures for projects SRP041338 and SRP091453, respectively. (C) Relative proportion of reads mapping to EBV, phiX, and all other metafeatures across RNA-seq samples. (D) Cumulative distribution plot of the average proportion of bacterial metafeature reads across all projects. Purple and pink vertical lines highlight projects SRP041338 and SRP091453, respectively. samples (see Methods for details). The correct infection agent showed the most significant difference across all metafeatures between infected and control samples for each study (Fig.2). For example, Westermann et al. [27] generated dual RNA-seq data from HeLa cells infected with the enteric bacterial pathogen S. enterica serovar Typhimurium and compared them to mocktreated control samples. Accordingly, we observed S. enterica as the most differentially abundant metafeature between the infected and the control samples (P<1e-75, Fig. 2A). Likewise, we recovered alphapapillomavirus 9,human alphaherpesvirus 1 (also known as herpes simplex virus 1), and rhinovirus A as the most differentially abundant metafeatures in the data from Zhang et al. [28], Rutkowski et al. [29], and Bai et al. [30], respectively. In the Westermann et al. [27] and Rutkowski et al. [29] studies, several additional metafeatures showed a strong differential abundance effect (Fig. 2A and 2C). These metafeatures were closely related to the true infection agent, i.e., Salmonella bongori (P<1e-67) and Panine alphaherpesvirus 3 (P<1e-9) for the Westermann et al. [27] or Rutkowski et al. [29] study, respectively. These findings confirm that our MetaMap pipeline recapitulates results from dedicated dual RNA-seq studies, i.e., studies based on known infectious agents. Therefore, MetaMap may be equally suited to detect previously unknown microbial and viral species in human primary samples. As an additional control, we re-analyzed two projects contained in our data collection that are derived from the B lymphoblast cell line under noninfectious conditions. However, since Epstein-Barr virus (EBV) is used for transfection and transformation of lymphocytes to lymphoblasts, we expected to detect reads from this virus in these projects [31], but no further viral or microbial reads [32]. Indeed, the most abundant metafeatures in each project were dominated by reads classified to gammaherpesvirus 4 (also known as EBV) and Enterobacteria phage phiX174 sensu lato (phiX), commonly used as spike-in in Illumina sequencing runs [33](Fig.3A and 3B). On average, 95% and 97% of all metafeature reads were classified as phiX or EBV for projects SRP041338 and SRP091453, respectively (Fig. 3C). Conversely, the abundance of reads mapping to bacterial species for these two projects corresponds to the bottom percentile as compared to all other projects in the MetaMap database, supporting sterility of this cell line (Fig. 3D). This demonstrates that MetaMap not only is capable of rediscovering known pathogenic species (true positives) in controlled infection experiments (Fig. 2) but it also minimizes the detection of false positives or, at Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 6MetaMap 0 5000 10000 15000 0 5000 10000 15000 Salmonella enterica levels CLARK-S based abundance (RPM) BLAST based abundance (RPM) Spearman Rho 0.98 P < 2.2e-16 1 5 50 500 5000 1 10 100 1000 10000 Average metafeature abundance CLARK-S based (RPM) BLAST based (RPM) Spearman Rho 0.16 P 3.1e−10 AB C Salmonella enterica CLARK−S BLAST 0246 Reads per hour per thread (log10) 8 Figure 4: Alternative BLAST-based classification method validates metafeature abundance estimates by MetaMap. (A) Average metafeature reads per million levels derived using the CLARK-S software, as implemented in the MetaMap pipeline, and a BLAST-based alternative approach on the xand yaxes, respectively. (B) Correlation in S. enterica abundance levels between the two classification approaches. (C) Difference in classification speed between the BLAST and CLARK-S metatranscriptomic classification. The yaxis shows the number of reads processed per hour per thread in log10 space. least, provides measures such as abundance and significance, allowing the user to identify and counterselect those species. As a technical validation, we compared our approach to an alternative metatranscriptomic classification strategy for the Westermann et al. [34] study. All non-human reads were aligned using BLASTN to a BLAST database consisting of the same genomic sequences used by CLARK-S (see Methods for details). The average metafeature abundances across all 42 samples derived from the BLAST-based approach and CLARK-S correlated significantly (Spearman correlation, Rho: 0.16, P: 3.1e-10) (Fig. 4A). BLAST showed higher sensitivity and detected more metafeatures compared to CLARK-S (indicated by the accumulation of dots at value 0 on the xaxis in Fig. 4A).Thisismostlyobserved for low abundance metafeatures that could represent low counts derived from sequencing and/or mapping errors. However, most importantly, the true pathogen metafeature “Salmonella enterica” showed very high correlation across samples between the BLASTand CLARK-based abundance estimates (Fig. 4B). Noteworthy, the MetaMap pipeline processed reads more than three orders of magnitude faster than BLAST, demonstrating a significant speed advantage while generating comparable results (Fig. 4C). Re-use potential Microbial and viral contamination in next-generation sequencing data has been observed. It can be caused by mapping errors due to genome sequence similarity between different species [35,36]. In addition, technical confounders can obstruct the analysis and potentially generate artificial differences if not considered properly. For example, different types of human samples may contain different amounts of non-human material due to varying sterility of the tissues. Furthermore, sequencing depth may introduce a detection floor for metafeatures that are not abundant. Therefore, comparisons across different tissues and sequencing depths may generate artificial differences. Additionally, given that only uniquely discriminative sequences are counted, the absolute abundance levels may not be comparable across metafeatures. Finally, the MetaMap pipeline captures metafeature abundance at the RNA level, which may not necessarily correspond to genomic abundance levels. Metafeatures Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 Simon et al. 7 may not be abundant at the DNA level but highly transcriptionally active and thus abundantly detected at the RNA level, or the inverse. These potential challenges need to be taken into consideration when comparing across metafeatures. To minimize these effects, we encourage focusing on studies that include intraproject comparisons that test one metafeature at a time, as exemplified in the differential metafeature abundance analysis. Our rationale is that technical confounders, in contrast to biologically meaningful changes, should affect all runs within a project to the same extent and therefore not show condition-specific effects. For example, in the Westermann et al. study [34], we detected substantial levels of phiX in both conditions (infected samples and mock-treated controls), but only the “Salmonella” metafeature showed a condition-specific effect. We aim to address the challenges inherent to interproject and intermetafeature comparisons in future work. All the raw data described in the present study were publicly available, yet have been very cumbersome to extract individually. The presented MetaMap database makes these data easily accessible for a very broad community, thereby allowing for global comparisons over hundreds of individual studies and thousands of sampled conditions. While we attempted to minimize the risk of detecting false positives (Fig. 3), it should be noted that not all metafeatures classified by MetaMap will necessarily refer to true biological factors. Noteworthy, our approach reveals a correlation between metafeatures and disease, not causality, and cannot discriminate disease-associated effects from potential treatment effects. However, our pipeline provides the user with a scientific starting ground to validate the presence/absence of defined microbial and viral species under defined conditions and explore the underlying biology and significance in greater detail. As a potential use case of these data, users can test for associations of microbial or viral metafeatures with a plethora of human diseases or between themselves. In addition, users with interest in a specific bacterial or viral species can easily identify studies and, consequently, disease contexts in which reads from this organism were detected. This could give an important first hint to assess whether the respective species might be implicated in a given human disease etiology. Furthermore, this resource provides the opportunity to support findings derived from standard microbiome profiling technologies, such as 16S rRNA gene based or shotgun metagenomics [37]. Finally, metafeature detection in human clinical RNA-seq samples may provide a diagnostic advantage when studying microbes or viruses that are challenging to isolate. The composite metafeature OTU count table, derived from 17 278 cDNA libraries from 436 SRA projects, including annotations is provided for download [38]. Availability of source code and requirements Project name: MetaMap Project home: https://github.com/theislab/MetaMap Operating system(s): Platform-independent Programming language: Unix command line, R Other requirements: STAR and CLARK-S may require large amounts of memory (>100 GB) License: GNU GPL Availability of supporting data The datasets supporting the results presented here are available in the GigaScience Database repository [38]. The protocols are also available at [40]. Additional file Additional File1.csv Abbreviations BLAST: basic local alignment search tool; EBV: Epstein-Barr virus; LRZ: Leibniz Supercomputing Centre; OTU: operational taxonomic unit; phiX: Enterobacteria phage phiX174 sensu lato; RNAseq: RNA sequencing; SRA: Sequencing Read Archive; SRP: short read project; STAR: Spliced Transcripts Alignment to a Reference software. Competing interests The authors declare that they have no competing interests. Funding L.S. acknowledges funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement (753039). The operation of the LRZ Linux Cluster is funded via the Bavarian State Ministry of Education, Science, and the Arts. Author contributions Conceptualization: L.S., M.E., L.D., and B.H.; formal analysis: L.S., M.H., S.K., and A.E.; investigation: L.S., A.J.W., and M.E.; methodology: L.S., S.K, M.H.; writing the original draft: L.S. and A.J.W; writing, reviewing, and editing: L.S., A.J.W., M.E., A.E., L.D., M.H., and F.T.; supervision: L.D., M.H., and F.T. Acknowledgments The authors thank Yu Wang and Ferdinand Jamitzky from the LRZ for their support. References 1. Young VB. The role of the microbiome in human health and disease: an introduction for clinicians. BMJ 2017;356:j831. 2. Turnbaugh PJ, Ley RE, Mahowald MA, et al. An obesityassociated gut microbiome with increased capacity for energy harvest. Nature 2006;444:1027–31. 3. Henao-Mejia J, Elinav E, Jin C, et al. Inflammasome-mediated dysbiosis regulates progression of NAFLD and obesity. Nature 2012;482:179–85. 4. Cani PD, Bibiloni R, Knauf C, et al. Changes in gut microbiota control metabolic endotoxemia-induced inflammation in high-fat diet-induced obesity and diabetes in mice. Diabetes 2008;57:1470–81. 5. Wang Z, Klipfell E, Bennett BJ, et al. Gut flora metabolism of phosphatidylcholine promotes cardiovascular disease. Nature 2011;472:57–63. 6. Engel M, Endesfelder D, Schloter-Hai B, et al. Influence of lung CT changes in chronic obstructive pulmonary disease (COPD) on the human lung microbiome. PLoS One 2017;12:e0180859. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019 8MetaMap 7. Kostic AD, Gevers D, Pedamallu CS, et al. Genomic analysis identifies association of Fusobacterium with colorectal carcinoma. Genome Res 2012;22:292–8. 8. Castellarin M, Warren RL, Freeman JD, et al. Fusobacterium nucleatum infection is prevalent in human colorectal carcinoma. Genome Res 2012;22:299–306. 9. Kodama Y, Shumway M, Leinonen R, International Nucleotide Sequence Database Collaboration. The Sequence Read Archive: explosive growth of sequencing data. Nucleic Acids Res 2012;40:D54–6. 10.Conesa A, Madrigal P, Tarazona S, et al. A survey of best practices for RNA-seq data analysis. Genome Biol 2016;17:13. 11.Gouin A, Legeai F, Nouhaud P, et al. Whole-genome resequencing of non-model organisms: lessons from unmapped reads. Heredity 2015;114:494–501. 12.Peng X, Wang J, Zhang Z, et al. Re-alignment of the unmapped reads with base quality score. BMC Bioinformatics 2015;16(Suppl 5):S8. 13.Westermann AJ, Gorski SA, Vogel J. Dual RNA-seq of pathogen and host. Nat Rev Microbiol 2012;10:618–30. 14.Westermann AJ, Barquist L, Vogel J. Resolving host-pathogen interactions by dual RNA-seq. PLoS Pathog 2017;13:e1006033. 15.Juranic Lisnic V, Babic Cac M, Lisnic B, et al. Dual analysis of the murine cytomegalovirus and host cell transcriptomes reveal new aspects of the virus-host cell interface. PLoS Pathog 2013;9:e1003611. 16.Xu G, Strong MJ, Lacey MR, et al. RNA CoMPASS: a dual approach for pathogen and host transcriptome analysis of RNAseq datasets. PLoS One 2014;9:e89445. 17.Park S-J, Kumar M, Kwon H-I, et al. Dynamic changes in host gene expression associated with H5N8 avian influenza virus infection in mice. Sci Rep 2015;5:16512. 18.Saxena K, Simon LM, Zeng X-L, et al. A paradox of transcriptional and functional innate interferon responses of human intestinal enteroids to enteric virus infection. Proc Natl Acad Sci 2017;114:E570–9. 19.Wesolowska-Andersen A, Everman JL, Davidson R, et al. Dual RNA-seq reveals viral infections in asthmatic children without respiratory illness which are associated with changes in the airway transcriptome. Genome Biol 2017;18:12. 20.Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2012;29:15–21. 21.Ounit R, Lonardi S. Higher classification sensitivity of short metagenomic reads with CLARK-S. Bioinformatics 2016;32:3823–5. 22.Lindgreen S, Adair KL, Gardner PP. An evaluation of the accuracy and speed of metagenome analysis tools. Sci Rep 2016;6:19233. 23.Engstr ¨ om PG, Steijger T, Sipos B, et al. Systematic evaluation of spliced alignment programs for RNA-seq data. Nat Methods 2013;10:1185–91. 24. www.lrz.de/services/compute/linux-cluster, Leibniz Supercomputing Centre 25.Altschul S. Basic local alignment search tool. J Mol Biol 1990;215:403–10. 26.Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 [Internet]. Genome Biol 2014;15:550. 27.Westermann AJ, F¨ orstner KU, Amman F, et al. Dual RNA-seq unveils noncoding RNA functions in host–pathogen interactions. Nature 2016;529:496–501. 28.Zhang Y, Koneva LA, Virani S, et al. Subtypes of HPV-positive head and neck cancers are associated with HPV characteristics, copy number alterations, PIK3CA mutation, and pathway signatures. Clin Cancer Res 2016;22:4735–45. 29.Rutkowski AJ, Erhard F, L’Hernault A, et al. Widespread disruption of host transcription termination in HSV-1 infection. Nat Commun 2015;6:7126. 30.Bai J, Smock SL, Jackson GR, Jr, et al. Phenotypic responses of differentiated asthmatic human airway epithelial cultures to rhinovirus. PLoS One 2015;10:e0118286. 31.Santpere G, Darre F, Blanco S, et al. Genome-wide analysis of wild-type Epstein–Barr virus genomes derived from healthy individuals of the 1000 Genomes Project. Genome Biol Evol 2014;6:846–60. 32.Mangul S, Olde Loohuis LM, Ori A, et al. Total RNA sequencing reveals microbial communities in human blood and disease specific effects , bioRxiv. 2016. doi:10.1101/057570. 33.Mukherjee S, Huntemann M, Ivanova N, et al. Large-scale contamination of microbial isolate genomes by Illumina PhiX control. Stand Genomic Sci 2015;10:18. 34.Westermann AJ, F¨ orstner KU, Amman F, et al. Dual RNA-seq unveils noncoding RNA functions in host-pathogen interactions. Nature 2016;529:496–501. 35.Strong MJ, Xu G, Morici L, et al. Microbial contamination in next generation sequencing: implications for sequence-based analysis of clinical samples. PLoS Pathog 2014;10:e1004437. 36.Bonfert T, Csaba G, Zimmer R, et al. Mining RNA–seq data for infections and contaminations. PLoS One 2013;8:e73071. 37.Cox MJ, WO C, Moffatt MF. Sequencing the human microbiome in health and disease. Hum Mol Genet 2013;22:R88–94. 38.Simon LM, Karg S, Westermann A, et al. Supporting data for “MetaMap: an atlas of metatranscriptomic reads in human disease-related RNA-seq data.” GigaScience Database 2018. http://dx.doi.org/10.5524/100456. 40.Simon LM, Karg S. MetaMap pipeline. protocols.io 2018; doi:dx.doi.org/10.17504/protocols.io.msec6be. 41 Tange O, GNU Parallel - The Command-Line Power tool. The USENIX Magazine 2011;36:42–47. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy070/5036539 by Ges f Biotechnologische user on 11 December 2019