The workflow is illustrated in Figure 1. Raw reads can be uploaded in compressed FASTQ format. Optionally, read runs from the NCBI Sequence Read Archive (SRA) can be added by either selecting from the 700 existing cohorts of the DAnIEL database or by entering their accessions directly. Metadata about the samples can be uploaded in CSV or Excel file format if statistical analysis is needed. Parameter sets for tweaking the workflow can be created, e.g., to filter features by abundance or to trim custom primer sequences. A comprehensive documentation about parameters to tweak the workflow and a tutorial is available on the DAnIEL web server. To facilitate biological insights in significant features, we constructed a manually curated database containing 1,669 fungal interactions with diseases, bacteria species and immune components retrieved from 761 published papers. This database is used by the web server for biological interpretation of significantly different abundant or correlated taxa. Furthermore, we incorporated a list of clinical samples of species involved in fungal infections from the German National Reference Center for Invasive Fungal Infections (NRZMyk).

Features are generated from the raw reads provided. Samples are demultiplexed, if necessary, according to the barcode mapping provided in the metadata table. External samples are downloaded from the NCBI Sequence Read Archive using grabseqs (Taylor et al., 2020). Quality control (QC) is performed afterwards. FastQC and MultiQC are used to monitor sequencing errors (Ewels et al., 2016). Cutadapt is used to trim primer and adapter sequences (Martin, 2011). Samples can be excluded from downstream analysis using various criteria such as minimum number of quality-controlled reads or base quality tests specified by FastQC. Representative biological sequences are created from quality-controlled reads via denoising. Either OTUs or amplicon sequence variants (ASVs) can be called using PIPITS (Gweon et al., 2015) or DADA2 (Callahan et al., 2016), respectively. Taxonomy of denoised sequences is assigned using either Naive Bayes or BLAST consensus approach of QIIME2 (Bolyen et al., 2019). Abundance counts are pooled at any given taxonomic rank and filtered by abundance and prevalence. Lastly, pooled counts are normalised using the methods aware of different library sizes like rarefaction or cumulative sum scaling (CSS), as implemented in the R packages vegan and metagenomeSeq, respectively (Dixon, 2003; Paulson et al., 2013). Centered log-ratio (CLR) normalisation is used by default to account for the compositionality. The generated features are used in downstream analysis to infer biological insights.

DAnIEL was initially run on three cohorts to retrieve a list of fungal species relevant for analyses when studying human samples: Faecal samples from mycobiome datasets of cancer patients (N=71, ITS2, PRJEB33756; Mirhakkak et al., 2021), antibiotics intervention (N=59, ITS2, PRJNA579284; Seelbinder et al., 2020) and human skin swab samples (N=203, ITS1, PRJNA286273; Leung et al., 2016). For each species found, we constructed a NCBI Entrez query to search for PubMed abstracts. Terms “disease”, “cytokine”, “immune system” and “prokaryote” and a limit of 20 papers per species were used to narrow down the focus of our subsequent manual curation. In total 1,337 abstracts from these papers were reviewed to create a manually curated database of fungal interactions. Medical Subject Headings (MeSH) were used for annotations whenever applicable. In addition, FUNGuild was integrated to provide information about the trophic modes in an ecological context (Nguyen et al., 2016).

Diversity is calculated using the R packages vegan (Dixon, 2003) and phyloseq (McMurdie and Holmes, 2013). Various methods, including principal coordinates analysis (PCoA) and non-metric multidimensional scaling (NMDS), can be used to generate ordination plots. FastSpar implementation of the SparCC algorithm can be used to create correlation networks of co-abundant taxa (Friedman and Alm, 2012; Watts et al., 2019). Alternatively, BAnOCC can be chosen to account for the compositionality of NGS abundance data (Schwager et al., 2017). The correlation analysis can be executed for each sample group individually, e.g., to compare networks of “case” and “control” samples. If a metadata table is provided, group-wise statistics are performed using Mann–Whitney U test for binary response variables and Kruskal–Wallis one-way analysis of variance in combination with Dunn’s post hoc test (Dunn, 1964) for other nominal responses. Spearman’s rank correlation is used for continuous responses instead. Features significant in any of these tests are annotated with our manually curated database of fungal interactions and clinical samples. Machine learning is applied to categorical response variables using the R package caret (Kuhn, 2008). Both random forest (RF) and support vector machines (SVMs) are used in combination with ANOVA filter and recursive feature selection. Best performing models according to the area under the receiver operating curve (AUC) in 5-fold cross-validation and feature importance scores are reported.

The overall pipeline of the DAnIEL web server consists of two parts: A front-end the user is interacting with to upload and visualise the data and a back-end workflow responsible for processing the uploaded data. The front-end of DAnIEL web server is implemented as an R shiny app. For visualisation ggplot2 is used (Wickham et al., 2019) and Rmarkdown to create summary reports. The back-end is built as a Snakemake workflow (Koster and Rahmann, 2012). This allows running the workflow separately on any Linux system including computing clusters. Conda is used to create reproducible environments for installing and running scripts and individual tools. A unique identifier will be assigned to each project to access the results later on. This also acts as a token for authentication. The tutorial consisting of 38 samples usually takes approximately half an hour wall time using 10 threads to be fully processed. Reports and visualisations can be accessed at the front-end once the corresponding step in the workflow has finished. This includes interactive plots and a summary consisting of findings, annotations, methods and references in a single HTML file.

An overview on related software packages for analysing fungal amplicon sequencing data is given in Table 1. QIIME2 is a command-line focused tool; therefore, it is not ideal for researchers without programming skills (Bolyen et al., 2019). ITScan covers profiling of operational taxonomic units (OTUs), however it does not cover quality control of raw reads (Ferro et al., 2014). CloVR-ITS was designed for pyrosequencing data; whereas DAnIEL is built for illumina paired-end data (White et al., 2013). Most tools are lacking the ability to calculate correlation networks especially those aware of the compositional nature of taxon counts, which is crucial in most analyses (Gloor et al., 2017). Tools like PipeCraft and LotuS focus on calculating the OTU table (Hildebrand et al., 2014; Anslan et al., 2017). Many existing tools are general and do not account for properties of a typical ITS dataset by default. For example, the length of ITS1 can range from 9 to 1,181bp (Yang et al., 2018). We chose 50bp as the default minimal QC read length as a trade-off to be able to detect fungi with a short ITS region while still have enough bases left for an accurate taxonomic classification. To the best of our knowledge, DAnIEL is the only web server available covering the whole workflow of ITS analysis beginning from raw reads to publication ready figures and tables, as well as, integration with a relational database for biological evaluation of statistical findings and comparative analysis with public available mycobiome data sets.

We demonstrated the functionality of the DAnIEL web server by running it on public available cohorts of soil and human mycobiomes. The first cohort investigated the effects of wildfire on the soil fungi in northwestern Canadian boreal forest (N=300, NCBI Accession PRJNA564811; Whitman et al., 2019). The cohort was selected from the integrated database of fungal projects. DAnIEL was run with default parameters. The results are shown in Figures 2, 3. All sub figures were directly generated by the web server. We confirmed that fungal communities were strongly dissimilar between burned and unburned sites. Burned sites showed significantly decreased Shannon and Chao1 alpha diversity metrics (Wilcoxon rank sum test, p<0.01). Fungal Bray-Curtis dissimilarities differ also significantly (Adonis PERMANOVA, p<10⁻⁴). Furthermore, a disrupted co-abundance pattern was observed in burned sites using SparCC correlation networks. Genus node degree and betweenness centrality are significantly decreased (Wilcoxon rank sum test, p<0.01). We increased the minimal absolute correlation coefficient to rmin,SparCC=0.3, in the interactive network exploration of DAnIEL to emphasise this coabundance fragmentation. A Random Forest was picked to be the best model in predicting the state (burned vs. unburned) based on the fungal abundance profile (AUC=98% in 5-fold cross validation).

Secondly, we performed a meta-analysis of two publicly available human skin mycobiomes about dandruff (N=966, ITS1, PRJNA415710; Saxena et al., 2021) and chronic wounds (N=384, ITS1, PRJNA324668; Kalan et al., 2018). Results of the human case study are shown in Supplementary Figures 1, 2. Wound sites showed significantly decreased Shannon and Chao1 alpha diversity metrics (Wilcoxon, p<0.001). Fungal Bray-Curtis dissimilarities differ also significantly (Adonis PERMANOVA, p<10⁻⁴). Wound samples showed increased abundances in Malassezia and Saccharomyces. A fungal co-abundance network was only possible to be constructed for dandruff samples after filtering with default parameters. The minimum absolute SparCC correlation coefficient was lowered to 0.1 to check the robustness of the networks. This revealed three co-abundant genera pairs in the wound samples, but this network was still much sparser than the network obtained from dandruff samples. A random forest model was picked to be the best one in predicting the skin type based on the fungal abundance profile (AUC=99% in 5-fold cross validation) with Saccharomyces and Ascomycota spp. Showing high Gini feature importance. The higher AUC value compared to the soil example indicated a stronger non-linear biological signal discriminating the sample groups using the fungal abundance profile.

Taxa found significant in either of differential abundance or co-abundance analysis were annotated with FUNGuild (Nguyen et al., 2016) and our integrated relational databases. The genera found significantly different in the human cohort were assigned with 867 interactions to other bacteria or cytokines and 141 infection related samples from our manually curated database (Supplementary Figure 3). The number of annotations varied across the clades according to their coverage of published literature. For the soil cohort, however, the genera found significant showed only 214 interactions reported in the literature and 25 infection related samples being reported confirming the human focus of our database.

Processing durations were benchmarked on a Docker container provided with 10 cores and 100GB of RAM. It took 2.9h to process the soil samples and 62.9h to process the human samples (see Supplementary Figure 4). The time limiting step in big cohorts of the pipeline is the denoising process, which can be parallelised much better for ASV profiling compared to OTU profiling. Downloading and quality control takes usually less than a minute per sample.

The performance of taxonomic profiling was evaluated using simulated reads. Grinder was used to simulate 10 samples in which each of them consists of 500×103 150bp paired-end (PE) reads (Angly et al., 2012). Primers ITS1 and ITS2 targeting the ITS1 sub region were utilised to simulate abundances from 100 different reference sequences following an exponential distribution. The same database UNITE 8.2 dynamic was used for simulation and training the classifiers to enable a fair comparison (Nilsson et al., 2019). To simulate biological variability, a uniform mutation rate of 1% was incorporated, in which substitutions were four times more likely than insertions or deletions. DAnIEL was run on the data with different methods for denoising (DADA2 and PIPITS) and taxonomic classification (BLAST consensus and Naive Bayes). Benchmarking performance of profiling abundance was based on (Ye et al., 2019) using counts pooled at genus rank. Briefly, the Euclidean distance (L2 norm) between the measured and the true abundance profile was calculated for each sample. Furthermore, differences of abundances were calculated for each sample and taxon separately. Pure taxon occurrence was evaluated by counting samples, in which a taxon was both measured and simulated. Precision, Sensitivity, Specificity and F1 score were calculated based on this contingency table. Taxon occurrences were very specifically but less sensitively profiled (98 and 33% on average, respectively, see Supplementary Figure 5). DADA2 outperformed PIPITS in all metrics. The outperformance of using ASV compared to OTU is consistent with the literature (Callahan et al., 2017; Caruso et al., 2019). ASVs are more accurate and allow a more detailed analysis of the mycobiome studied. However, this can be also disadvantageous. The detailed taxonomic profile of ASVs in taxonomically diverse studies like environmental samples can make a manual curation of sequence alignments and statistical downstream analyses more difficult. Furthermore, intragenomic variation can result in multiple ASVs originating from the same fungal cell overestimating the true diversity (Schoch et al., 2012). On the other hand, OTU profiling is less sensitive to potentially unwanted details, and we still considered it as a necessary function in DAnIEL to make the results comparable with older studies.

Naive Bayes classification outperformed the BLAST consensus approach in terms of specificity and precision but not in sensitivity. Most accurate abundance profiles were generated using DADA2 (see Supplementary Figure 6). PIPITS underestimated many abundances, which lead to increased distances in some samples especially in combination with the Naive Bayes classifier. The values were very similar for different phyla. Therefore, taxonomy seemed to have only little influence on the classification performance.

Furthermore, we compared the results in the case studies using both ASV and OTU profiling. The results are shown in Supplementary Figures 7, 8. Most of the high abundant taxa were found using any denoising method with similar abundance values and correlation networks. The alpha diversity was significantly higher using OTU profiling (Wilcoxon rank sum test, p<0.001).

We evaluated the generalisability of the DAnIEL interactions database on 30 other ITS studies from various habitats. The results are shown in Supplementary Figure 9. Only genera prevalent in at least 10% of the samples in any habitat were considered for this analysis. On average, 29% of the prevalent genera in host habitats, 28% in aquatic and 15% in soil samples were already found in our manually curated database.