dataesr/scientific-paragraphs-categorization huggingface dataset:
We noted also that Si-O-Si angles exhibit a much larger variation of values compared to the distances. This is certainly one of the reason why the metrics on angles are much more important than metrics on distances. In the case of SOAP descriptors, distributions seem to be shifted compared to the PCOD2 database, suggesting that the hypothetical environments may be too diverse compared to feasible zeolites -and also showing that even though using SOAP descriptors resulted in the best predictions throughout this work, they are not suited for a prediction on the PCOD2 database. The geometrical descriptors were then used to predict the thermal expansion and bulk modulus of the PCOD2 database, which highlighted a trend where zeolites with a higher bulk modulus tends to have larger negative thermal expansion.
We introduce a homemade pipeline allowing to improve the quality of the metagenomic annotations carried out when using shallow depth metagenomic datasets. The main motivation being to be able to quantify more precisely, with greater certainty, the genes involved in bacterial interactions. The limitation in our experimental design is that we use a sequencing technique with a low throughput (miSeq) compared to the metagenomic standard (hiSeq) because we carry out a fairly large sampling (almost a hundred samples) in time series. This methodological constraint from our study means that the assembly of the sequences is not very exhaustive (less than 50% of the sequences manage to be assembled). In this chapter, we will therefore present a new pipeline designed to specifically deal with such kind of data.
We used co-assembly and a sequence annotation strategy in order to recover the sequences that could not be mapped on the assembled contigs. In addition, in order to avoid adding too much noise, when rescuing reads, we have built an algorithm to define a threshold of e-value based on the noise of the sequence annotation learned from sequences mapped in the assembly. We have selected several recent tools known to be effective for assembling, mapping and annotating these data. In addition, this pipeline was also built in order to be very user-friendly in terms of installation. In this idea of reproducibility, accessibility and transparency, we have designed an installation script to allow each user to install each tool required for the pipeline in a simple and reproducible way. Regarding the performances of this pipeline, we were able to show that the expected error rate (False discovery rate) for the annotation was close to 5%.
Introduction Metagenomic approaches are useful for investigating both the diversity and functioning of environmental microbial communities. Over the last decade, several tools and workflows have been released to assemble and analyse these datasets (e. g. Li et al. 2016; Bankevich et al. 2012; Wood and Salzberg 2014; Buchfink, Xie, and Huson 2015; Menzel, Ng, and Krogh 2016). Metagenomic data (derived from shotgun sequencing of total extracted DNA) has recently been used to assemble putative genomes, called metagenomic assembled genomes (MAGs), to determine taxonomy and to identify metabolic functions of microbial communities. Sequence reads are assembled into contigs, which improves the accuracy of annotation, and binned into MAGs. Both the contigs and the bins (putative genomes -MAGs) can provide more accurate taxonomical annotations.
The accuracy of the processes suggested for assembling the reads in contigs and for binning them are dependent on sequencing depth; if it is too shallow, the coverage of the contigs will be low and only a small fraction of reads will be recruited into the assembly. As a consequence, the results might be less reliable and other strategies must be used. There are several methods for improving the assembly into contigs and for the subsequent binning into MAGs. Co-assembly is one strategy that has been developed to improve the assembly quality when sequencing depth is limited. If the metagenomic dataset is composed of several samples that are known to contain similar bacterial communities, they can be pooled during the assembly step to increase the sequencing depth and improve the quality of the assembly.
In this case, reads from different samples can be used to create a contig. The subsequent binning of contigs has often in the past relied on nucleotide frequency discrimination, although the binning could be improved by using differential coverage of the reads in the contigs in the different samples. Co-assembly of the contigs from multiple samples could, therefore, improve the quality of MAGs retrieved by increasing the discrimination among taxonomically related species or strains. For example,  Delmont and Eren (2018)  retrieved more MAGs from the Tara ocean dataset, including several strains of Prochlorococcus, than found in the original publication  (Sunagawa et al. 2015 ) by using both co-assembly of the contigs and differential coverage for the binning. Co-assembly can also be used to track changes in taxonomy and functions in time-series datasets.
If the genes of interest are abundant, shallower sequencing technologies other than hiseq (e. g. miSeq) can be used and co-assembly can increase the assembly quality. However, pipelines to deal with these specificities (time series, shallow depth sequencing, low assembly) are currently missing. In this paper, we introduce a new workflow ("Eggvio") designed to annotate and analyse metagenomic datasets of any size (but optimized for shallow sequencing datasets). Coassembly was coupled to a read annotation strategy to rescue reads that were not easily mapped back to the assembled contigs. In addition, an algorithm was built to define an evalue threshold based on the noise of read annotation derived from the reads used in the assembly. Several recent tools to carry out the assembly, mapping, binning and annotation were selected to build a user-friendly pipeline.
Pipelines for metagenomics generally involve the installation of many tools or dependencies and require some previous informatics knowledge. To improve the reproducibility, accessibility and transparency, we designed a script that allows every user to install all the tools needed for the pipeline. EggVio is a flexible pipeline optimized for co-assembly, binning and annotation of low depth MiSeq samples (<10 6 reads per sample). Due to the computational needs for read annotation, it was designed to be used on a linux server with a batch slurm queue for job submission. It is coded mainly into bash scripts to process fastq data starting from quality filtering to functional annotations. R (R Development Core Team 2011) programming is used for merging annotations from contigs and reads that did not map back onto the assembly and also to determine the e-value threshold.
The programs used to carry out every step of the pipeline are shortly summarized in the flowchart shown in Figure  1. In addition, an installation script ( EggVio_install_tools. sh) is provided for most of the tools that do not require any root privileges for installation (this excludes the assembler, megahit). The user needs to download the annotation database for eggnog mapper and the taxonomic database for contig taxonomic annotation by kaiju separately. The databases were not included in the installation, because multiple annotation databases are available and the user can choose which one to download themselves. As a guide, we included both steps in the wiki pipeline (https://gitlab. com/R_addict/eggvio/wikis/Download-and-installation-of-EggVio). The scripts for the pipeline are freely available on GitLab (https://gitlab.
These data can then be imported into anvio for binning of MAGs and assembly visualization. For functional annotation of the metagenomes, several other steps are needed and are described in the next section with a focus on the algorithm used to learn the e-value threshold for read annotation. If the annotations of reads and genes are different, the read annotation is considered as spurious (False Positive = FP) and if they are identical, then they are considered as true positives (TP) F. The algorithm to learn the threshold is a greedy algorithm based on successive refining steps. It works on a (decreasingly) ordered vector of e-values from read annotation with a corresponding vector returning the information whether the corresponding e-value returned a TP or a FP annotation. The algorithm will compute the percentage of FP starting from the whole dataset and then by removing iteratively the highest 1000 e-value annotations.
When it finds the local optima, it will start computing the same statistics, but starting from the interval identified as optimal and refine it by removing iteratively 100 annotations and then for the last refinement step only 1. This e-value returned will then be used to annotate the reads that did not map onto genes with an expected probability of FP<0. 05. After co-assembly of the different samples of the dataset (Figure  2  A), gene detection on contigs (Figure  2 B ) is carried out using prodigal  (Hyatt et al. 2010 ). The coverage is then computed using Bowtie2 (Figure  2 C ) and results are converted into counts using a custom bash script relying on functions from samtools (H. Li et al. 2009 ) and bedtools  (Quinlan and Hall 2010). The genes and the reads which mapped onto them are then functionally annotated ( ).
These gene annotations and their respective mapped reads are compared to identify false positive (FP) annotations. If the read annotations differ from the genes they are mapped onto (considered as the 'gold standard' annotation), this is considered as a false positive (Figure  2 E ). Based on these data, the e-value threshold (E. T. ) learning algorithm is used to define a suitable threshold for read annotation where the percentage of expected FP in the annotations considered as significant would be p-value <= 0. 05. This algorithm is written in R and is presented below. #Variables: #dataF= a dataframe with zero (FP) or 1(TP) to inform about how good is the corresponding annotation #stepN= a numeric value to define the size of the interval to use to compute the FP #modeU= defines if the number of iterations has to be defined by the size of the dataset (nrow(dataF)) or fixed to a specific number of steps (10).
#startPos= the starting position (a row number in dataF) from where to start computing the FP for the first iteration falsePosFind = function(dataF, stepN, modeU, startPos){ if(modeU == 'raw'){ roundNumb <-nrow(dataF) -(nrow(dataF) %% stepN) stepMax <-roundNumb / stepN }else{ stepMax <-10 } res <-data. frame(rowN=rep(NA,stepMax+1), Fpos= rep(NA,stepMax+1)) for(iter in c(1:stepMax)){ if(iter == 1){ res$rowN[1] <-startPos res$Fpos[1] <- length(which(dataF[c(startPos:nrow(dataF)),1]==0))/(nrow(dataF)-startPos+1) } res$rowN[iter+1] <-(iter*stepN)+startPos fpos <- length(which(dataF[c(((iter*stepN)+startPos):nrow(dataF)),1]==0))/(nrow(dataF)- ((iter*stepN)+startPos)) if(length(fpos)==0){ res$Fpos[iter+1] <-0 }else{ res$Fpos[iter+1] <-fpos } } return(res) } #Core of the script #steps defines the size of the intervals for each iteration of the learning algorithm steps <-c(1000,100,10,1) #the dataframe from where to learn(learningData)is composed of two columns: V1= zero (FP) or 1(TP) to inform about how good is the corresponding annotation and V2= e-value associated to the annotation.
#The learning dataframe is first ordered by decreasing e-value (V2). dataForLearn <-([order(learningData$V2,decreasing = TRUE),] #pValThresh=0. 05 = fraction of FP tolerated for(stepSize in steps){ if(stepSize == 1000){ threshRes <-falsePosFind(dataForLearn,stepSize,'raw',1) }else{ threshRes <-falsePosFind(dataForLearn,stepSize,'default',whereToCut) } if(pValThresh < min(threshRes$Fpos)){ errorM <-paste('The p-value threshold set is smaller than the best result possible to reach for this dataset. Please set it above the following minimum possible:', min(threshRes$Fpos)) stop(errorM) } whereToCut <-threshRes$rowN[which(threshRes$Fpos < pValThresh)[1]] if (stepSize > 1){ whereToCut<-whereToCut -stepSize } #the e-value threshold final result is here threshold_seed_ortholog_evalue <-dataForLearn[whereToCut,2] exactPval <-round(threshRes$Fpos[which(threshRes$Fpos < pValThresh)[1]],5) } This function is designed to find an e-value threshold such that after filtering, the false discovery rate FDR =  + = 0.
05 (where FP = False Positive and TP = True Positive) At the same time, we would like to minimize the rejection of correct annotations. To meet both criteria, the function orders the dataset by decreasing e-value (the least significant evalue has the row name 1 in the ordered dataset). Then, the script calls the function "FalsePosFind" to compute the FDR with a threshold set every 1000 annotations. Since the data are ordered, the first threshold that meets the criteria FDR < 0. 05 will be the best solution, as it will preserve the highest amount of TP. Once this position is found, the loop will iterate the function "FalsePosFind" a second time, but starting from the position located 1000 annotations higher than that initially found. The optimization serves to minimize the amount of TP rejected as false negative (FN) and to refine the threshold such that it is less restrictive by computing FDR every 100 annotations starting 1000 annotations away from the local optima found on the previous iteration.
Benchmarking the EggVio pipeline Presentation of the dataset used in the study: To test how our pipeline could enhance the annotation of shallow sequenced datasets, we used one of our in-house datasets from a polluted site bioremediation project (MISS). It consisted of a time series of 30 samples tracking chlorinated compounds in a polluted ground water site. After the injection of organic carbon to induce the biodegradation of chlorinated solvents in the groundwater (C1), three other samples were taken one month apart (C2, C3 and C4). Six replicates were collected for each time point. The biodegradation of chlorinated compounds was evaluated at each time point and qPCR analysis was carried out in order to evaluate the abundance of the pceA gene coding for an enzyme involved in reductive dechlorination of tetrachloroethene (PCE) to trichloroethene (TCE).
Evaluation of the threshold learning algorithm: We first used this dataset to evaluate how our threshold learning for read annotation based on assembled data would perform. We computed the threshold for every annotation (gene id, gene name, Kegg orthologs = KO and Gene Onthology = GO) to compare the thresholds returned. We then focused on the KO annotation to test how learning on a subset of the dataset (since we could not assemble the whole dataset) would affect the threshold and FDR estimate. To test this, we randomly subsampled the learning dataset to predict the threshold and then observed how the FRD was affected when computed on the whole dataset.
Comparison of the sequencing results and the qPCR results on the tracking of the gene pceA We also used this dataset to compare how the hybrid annotation (genes and reads annotations) would affect the results as compared to gene only annotation using the intermediate results from EggVio. Both annotations (hybrid and genes only) were then normalized by the RPKM method using the R package GenomEnvironR (github: https://gitlab. com/R_addict/genomenvironr). The percentage of sample being annotated as genes and as reads was then investigated for every sample. A NMDS of the samples was carried out for both types of annotation at the KO (Kegg orthologs) level. These sample representations were compared to an NMDS plot of the samples that included chemical data of several chlorinated compound concentrations measured in the water.
Evaluation of the read annotation threshold learning algorithm In order to evaluate how our learning threshold performed, we first ran the algorithm on the fully assembled data. Then, we assessed how random subsampling of the training dataset would impact our learning threshold and its ability to keep the amount of false positive below 0. 05. To do this, we subsampled the mapped reads to smaller fractions (respectively 80%, 50 % and 10% of the original dataset) 1000 times and ran the threshold analysis. The threshold results were then compared to the amount of True Positive (TP) and False positive (FP) in the original dataset at different annotation levels (gene ID, gene name, KO and GO, Table  1 ). An FDR = 0. 05 was obtained for all the annotation levels except for Gene ID, where the error was too high with a threshold at FDR = 0.
The dashed vertical lines represent the median e-value of their respective distributions. We then evaluated the impact of dataset subsampling on the e-value threshold estimate for read annotation and our FDR estimate for the whole dataset. After plotting the results of 1000 subsampling at different percentages of the dataset, we observed that the e-value estimate became more variable at smaller subsampling sizes (Figure  5 ) and the boxplots tended to deviate from the true estimate toward smaller e-values as they were moved toward the bottom of the y axis (Figure  5 ). This trend was confirmed by the amount of TP observed for each of the e-value thresholds returned for the different subsampling. The second quartile of the boxplots was generally above 0. 95 (Figure  6 ) for the subsampling replicates, indicating that for more than 50% of the subsamples, the threshold returned was enriched with correct annotations and thus showed an FRD <0.
05. In addition, the replicates with the lowest amount of TP were always above 0. 945, and FRD <0. 055 was observed in the worst case during the subsampling experiments (Figure  6 ). 4. 2 The effect of the read annotation on the dataset representability and the detection of pceA genes. We investigated how the NMDS representations of the MISS dataset based on the metagenomic data agreed with the NMDS representation based on the chemistry (Figure  7 ). The NMDS from both types of annotation (assembly and hybrid) clustered more together at a given sampling time (e. g. , C1 with samples C1) than between samples from different sampling times, but their grouping at a finer scale showed differences. For example, the sample C1pz-20 was far from all the other samples from the dataset when considering only the assembly annotation (Figure  7  A) but did not display any chemical features on the chemistry NMDS that could explain its position (Figure  7 C ).
Discussion 5. 1 Assembly missed meaningful information in shallow depth datasets, but can be complemented by rescuing reads The percentage of reads recruited in the assembly varied widely across samples (Figure  3 ), with an average of 11% of reads recruited and the amount of gene ids detected in the assembly were 22 times lower than in the read annotations. This is lower than what was observed in a study where full datasets were assembled and 10% to 30% more annotations were retrieved using an assembly free method  (Anwar et al. 2019). The undersampling of the dataset was further supported by the NMDS representation that showed that sample C1pz-20 was the most unique sample based on the assembled data (Figure  7 A ). However, this sample had the lowest coverage (<1%) in the assembly (Figure  3 ).
After applying the hybrid method, the NMDS could be corrected (Figure  7 B ) and showed a representation of the samples that was much more in accordance with the chemical dataset (Figure  7 C ). In addition, the hybrid method increased the sensitivity of correlation detection between the relative abundance of pceA genes in the metagenomes and the concentrations of PCE over TCE (Figure  8  E versus F). This negative correlation was not detected in the qPCR data (p-val>0. 05). We interpreted this correlation as PCE degradation occurring prior to sampling time rather than an instantaneous measure of the PCE degradation potential of the bacterial community (rather detected by metatranscritomics). Thus, at lower PCE/TCE ratios, a higher abundance of pceA genes could be interpreted as a selection for organisms able to degrade PCE into TCE in the microbial community.
EggVio E. T. algorithm can accurately quantify the noise added when rescuing reads for annotation By randomly subsampling the learning dataset to evaluate how the e-value threshold estimate would affect the estimation of the FDR at the KO annotation level, we showed that in the worst case, the FDR would only be impacted by 5% leading to an accuracy (= TP/(FP+TP)) of 0. 945 (Figure  6 ). We also observed that the recall (= TP/(TP+FN) = 1 -FRACTION OF CORRECT REJECTED) was above 69% if we exclude the gene id (seed eggnog ortholog). These estimated performances can be compared to other read annotation tools such as miFaser  (Zhu et al. 2018) , recently released and based on a custom high quality database using a custom score modeled after the HSSP metric for function transfer between full-length proteins  (Schneider, de Daruvar, and Sander 1997).
Based on reads generated from their database, they assessed that the accuracy of miFaser could reach 90% and the recall 50%  (Zhu et al. 2018 ). However, since the databases are different and do not use the same criteria, comparisons are difficult to perform. Concerning the lowest level of annotation accessible by eggnog-mapper (gene id = seed eggnog ortholog), the noise is high and we would, therefore, not recommend using this annotation level for any downstream analyses. The high noise at this annotation level compared to others is likely related to the inclusion of taxonomical information. Gene id annotation is composed of a tax id at the species level and the gene id from its original genome, making it challenging to accurately retrieve annotations from non-assembled reads. Although not shown in this work, EggVio also integrates kaiju taxonomical annotation of the contigs, which is more accurate since it is carried out on longer reads, thus making the resolution much higher (as for functional annotation) and enables taxonomic analysis of the dominant strains of the population.
Conclusion We showed that EggVio is a flexible pipeline for processing shallow depth sequencing datasets. In addition to assembling reads into contigs, the pipeline can rescue unrecruited reads while adding a predicted amount of noise (FDR = 5%) to the annotation. EggVio improved the correlation between metagenomic and environmental chemistry data. In addition, it removed the artifacts observed in the NMDS when calculated using the assembly alone and showed a more reliable community composition and allow a more robust sample comparison in the actual dataset. Nonetheless, further validation, using a mock community for example or benchmarking it against other possible tools such as miFaser is required to validate the full performance of the pipeline.
The algorithm to learn the threshold is a greedy algorithm based on successive refining steps. It works on a (decreasingly) ordered vector of e-values from read annotation with a corresponding vector returning the information whether the corresponding e-value returned a TP or a FP annotation. The algorithm will compute the percentage of FP starting from the whole dataset and then by removing iteratively the highest 1000 e-value annotations. When it finds the local optima, it will start computing the same statistics, but starting from the interval identified as optimal and refine it by removing iteratively 100 annotations and then for the last refinement step only 1. This e-value returned will then be used to annotate the reads that did not map onto genes with an expected probability of FP<0.