dataesr/scientific-paragraphs-categorization huggingface dataset:
In particular, many problems in genomics map to a supervised learning framework with a goal of determining functional sets of genes from partial annotations and feature data. A correspondingly high number of methods and assessments report comparative evaluation using traditional machine learning statistics, such as the area under the receiver-operator characteristic curve (AUROC). However, genomics poses unique challenges and opportunities relating to the extreme scalability of data collection and analyses, both across novel contexts, such as conditions, tissues or species, and the ability to collect high-throughput data in consistent assays. The shared ancestry of organisms forms the basis of many ways we extend results from one system to another. Across species, this shared ancestry is the basis for functional annotation using homology  (Altenhoff et al.
ROC curves from the literature We systematically sampled 35 ROC-curve-containing research articles from the PLoS One journal during one calendar year (genomicsrelated Subject Areas) and selected 15 high-profile research articles (see Supplementary Appendix S1 for a detailed list of papers and figures extracted). We used the Engauge Digitizer (https://doi. org/10. 5281/zenodo. 3941227) software to extract curves from the selected figures, following the standard procedure (selection of 3 axis points for scale, automatic segment detection). In instances where the figure contained too many overlapping curves and individual curves proved too difficult to extract, we removed the figure from the analysis. After the extraction process, Engauge Digitizer generated CSV files with data points evenly distributed along the curve.
ROC curves from PPI, co-expression and co-domain data We downloaded the mouse PPI data BIOGRID-ALL version 4. 4. 197 from the BIOGRID website  (Oughtred et al. , 2019; Stark et al. , 2006). We filtered the BIOGRID data for mouse (taxonomy ID 10090) and physical interaction ('Experimental. System. Type' =¼ 'physical'). This initial network contained 57 337 interactions across 10 172 genes. To take into account indirect connections  (Gillis and Pavlidis, 2011b) , we propagated the existing interactions, setting the weight for each pair of proteins as the inverse of the shortest path between the two proteins. We downloaded the mouse co-expression network from the CoCoCoNet  (Lee et al. , 2020)  website (last updated on April 20, 2021). We converted Ensembl identifiers gene to symbols using the mapIds function from the AnnotationDbi and the org.
ROC curves from single-cell RNA sequencing data We downloaded the Tabula Muris  (Schaum et al. , 2018)  single-cell RNA sequencing (scRNAseq) dataset from FigShare, specifically Version 2 of the 10Â [Single-cell RNA-seq data from microfluidic emulsion (v2), 2018] and Smart-Seq2 [Single-cell RNA-seq data from Smart-seq2 sequencing of FACS sorted cells (v2), 2018] data, along with metadata and annotations, keeping all annotated cells (100 605 cells). We applied CP10K (counts per 10k) normalization for the 10Â data and CPM (counts per million) normalization for the SmartSeq data. We computed markers for each mouse by tissue combination using the MetaMarkers package  (Fischer and Gillis, 2021). We subset the datasets to a given mouse using the 'mouse_id' metadata, then ran the compute_markers function on the normalized counts, using the 'cell_ontology_class' as cell type labels, and 'tissue' metadata as group labels (stratifying marker search by tissue).
We removed genes with low detection rate (<10% in all mice by tissue by cell type combinations) and only kept markers inferred for cell types containing at least 20 cells. To compute ROC curves, we asked if a set of reference markers were the top markers in other mouse by tissue combinations. For each mouse by tissue combination, we ordered genes according to the effect size of the ROC test, 'auroc' column in the MetaMarkers table, then used this list as a predictor for the reference markers. We used the 'prediction' and 'performance' functions from the R ROCR package to compute the ROC curve ('tpr' and 'fpr' statistics) and the AUROC ('auc' statistic). We obtained a first set of reference markers by selecting the top 20 markers for the '3_10_M' mouse from the Smart-Seq dataset, for the 'B cell' cell type in the 'Fat' tissue.
Straight lines in ROC curves are pervasive in the genomics literature While reviewing the genomics literature, we were struck by the recurrence of straight lines in published ROC curves. To confirm our intuition that straight segments are surprisingly common, we extracted ROC curves from 50 research articles, composed of an unbiased selection of 35 articles from the PLoS One journal and 15 manually curated high-profile articles (Supplementary Appendix S1). In total, we extracted 77 ROC curves using the Engauge Digitizer (see Section 2) software (Fig. 1a ). The predictions summarized by the ROC curves were either gene-centric predictions (74/77 curves) or interaction predictions (3/77 curves). From the methods description and the figures, we estimated that the number of positives ranged from 24 to 6000, while the total number of objects ranged from 121 to 119 149.
Straight lines in ROC space suggest the presence of discrete functional modules Why are straight lines so common in genomics ROC curves? To provide an initial intuition, we consider a toy example of protein function prediction (Fig. 2 ). A machine learning classifier is applied to a high-throughput dataset measuring the likelihood that two proteins interact (Fig. 2a ). The classifier is trained on sparse annotations obtained by a low-throughput assay that labeled a subset of genes with 50% True Positive Rate (TPR) and 10% False Negative Rate (FNR) (Fig. 1a ). In this ideal scenario, the classifier identifies a functional module containing an even mix of previously annotated genes and unannotated functional genes. The presence of this functional module in the data is immediately visible on the ROC curve through the presence of straight lines.
We built a PPI network by aggregating all interactions from the BIOGRID  (Oughtred et al. , 2019; Stark et 2006)  database annotated as 'Mouse' and 'Physical Interaction', resulting in a network containing 10 172 proteins and 57 337 interactions. As PPI networks are typically sparse, we used a propagation algorithm to obtain a dense network, which accounts for indirect interactions between proteins (see Section 2). We downloaded the mouse coexpression network from the CoCoCoNet  (Lee et al. , 2020)  database. The network was obtained by aggregating 3359 samples over 85 experiments, resulting in a dense network containing 17 834 genes. To allow comparisons between the two modalities, we restricted the two networks to 9058 common genes. To assess whether a function is supported by a network's topology, we used the guilt-by-association framework implemented by the EGAD  (Ballouz et al.
Discussion In this study, we showed that the shape of ROC curves offers a visual and data-driven interpretation of the extent of biological functions. The presence of straight lines in the ROC curve suggests that the data are compatible with the extension of a functional gene set to unannotated genes. We call these straight lines FECs, because they define discrete classes of genes that are equivalent with respect to the functional property investigated. Our examples show that the extensibility of gene sets is context specific: we found that a subset of B-cell markers was conserved across tissues, while secondary markers varied from tissue to tissue. One of the strengths of FECs is that the generalizability and extensibility a gene set can be probed with one look at the ROC curve. Either the gene set works perfectly well in the new context (AUROC ¼ 1), or performance is suboptimal and FECs suggest how the gene set can be reorganized in the new context.
Our focus on large aggregate databases (PPI, co-expression, drug-target interactions) confirms the presence of widespread modular structure in each data type. However, recent algorithms combine increasingly broad data resources  (Gligorijevi cet al. , 2021; Kulmanov and Hoehndorf, 2020; Rifaioglu et al. , 2018; You et al. , 2021) , which may smooth out the strong modularity suggested by each data type alone, thus reducing the frequency of FECs. Another interesting avenue is the integration of high-throughput functional assessments generated from CRISPR screens  (Bock et al. , 2022) , in particular the elucidation of gene function at the single-cell resolution  (Replogle et al. , 2022) , to probe the generalizability of functional modules across cell types, tissues and conditions. In summary, FECs define a formal framework to visualize and probe the context-specificity of functional gene sets.
Data availability The datasets were derived from sources in the public domain: co-expression data were obtained from the CoCoCoNet website (http://labshare. cshl. edu/ shares/gillislab/resource/CoCoCoNet, 20 April 2021, date last accessed), protein-protein interaction data from BIOGRID (https://downloads. thebiogrid. org/Download/BioGRID, version 4. 4. 197), drug-target interaction from STITCH (http://stitch. embl. de/download, version 5. 0), domain information from UniprotKB (https://rest. uniprot. org/uniprotkb, 12 August 2022, date last accessed), gene ontology information from the gene ontology website (http:// current. geneontology. org/ontology/go-basic. obo, 19 December 2018, date last accessed) and the JAX website (http://www. informatics. jax. org/downloads/ reports/gene_association.
mgi. gz, 12 December 2018, date last accessed), data from the Tabula Muris consortium from FigShare (https://figshare. com/articles/ dataset/Single-cell_RNA-seq_data_from_microfluidic_emulsion_v2_/5968960/ 3; https://figshare. com/articles/dataset/Single-cell_RNA-seq_data_from_Smart-seq2_sequencing_of_FACS_sorted_cells_v2_/5829687/8).
haemophilum and M. tuberculosis was compared with the program MAUVE (23) (Fig. 3 ). Many differences were observed between the two leprosy bacilli, by one way, and between M. haemophilum and M. tuberculosis, by the other, indicating events of chromosomal rearrangement. To infer the minimal number of chromosomal rearrangements among them, a pairwise inversion distance matrix was constructed with the program GRIMM  (26) , which predicted 18 inversion events between M. leprae and M. lepromatosis and 11 between M. haemophilum and M. tuberculosis (Fig. 3 ). In addition, the MRCA of leprosy bacilli incurred ;50 chromosomal rearrangements after divergence from M. haemophilum. Comparisons of the consecutive CDS and pseudogenes in M. leprae and M. lepromatosis further rendered 25 syntenic blocks as the likely outcome of major rearrangement events (Table  S4 ).
lepromatosis. These chromosomal rearrangements may also produce gene duplications and deletions and affect gene expression. Similar effects were also previously observed for the genome of M. leprae  (25). Orthologous clusters of proteins were estimated with the program OrthoFinder v2. 3. 3  (27)  for proteomes of the four mycobacteria (Fig. 4 ). A core proteome of 1,332 protein clusters, with at least one member encoded in each genome, was detected. Of them, a set of 1,227 clusters included a single protein in each proteome (single copy orthogroups). The number of clusters with only proteins from one leprosy proteome was relatively high (82 in M. leprae and 300 in M. lepromatosis). However, most of them were annotated as hypothetical proteins, only a few with an annotated function (13 in M.
1  (71). Initial pseudogene prediction was based on Prokka annotation of adjacent CDS presenting the same annotation and a custom python script. In brief, the script uses LAST  (72)  alignments to detect CDS which are 75% shorter than a reference set of clustered proteomes at 95% amino acid identity (Table  S8 ). Then, it mixes adjacent fragmented genes when they overlap with the same reference hit. Pseudogenes annotations were revised exhaustively by using several types of BLAST analyses. The main one was based on BLASTX using the complete M. lepromatosis FJ924 genome sequence as a query against a protein database of M. haemophilum (E value 0. 01). This revision produced length changes of several pseudogenes, fusion of others, and conversion of some initially annotated CDS to pseudogenes.
GenoPlotR was used to plot Mauve results  (79). The permutation file was converted with the program GRIMM in a distant matrix, which contained the minimal number of inversions required between a pair of genomes to explain differences in gene order  (26). Distance matrix was charged in MEGA7 (76) and a neighbor joining algorithm was used to infer the rearrangement phylogeny  (80). Orthology between CDS and pseudogenes from M. leprae Br4923 and M. lepromatosis FJ924 was obtained with a reciprocal BLASTN best hit strategy (E value = 1. 0E-05). Nonreciprocal hits and those hits in disagreement with gene order were revised upon visual inspection. Synteny was the first criterion to assign orthology when more than one similar hit was detected allowing up to 5 hits in new BLASTN searches. Identification and analysis of repetitive sequences.
3. 3 using the msa (mafft aligner) and IQ-TREE options  (27, 73, 74). Data on the clusters were used for several analyses. Estimation of the ages of pseudogenes. The relative ages of M. leprae and M. lepromatosis pseudogenes were estimated according to the method of Gomez-Valero et al. (7)  with inclusion of the recent M. haemophilum genome data  (12)  for closer phylogeny. Briefly, this method used M. tuberculosis as outgroup to root the tree and compared M. haemophilum with either M. leprae or M. lepromatosis. The divergence between the two species determines the relative age of pseudogenization, ranging from 1, when it took place at the start of the divergence, to 0 at this time. The method is based on the estimation of the dN and dS (the number of nonsynonymous and synonymous substitutions per site, respectively) for each cluster of orthologous gene-pseudogene in the branches leading to M.
haemophilum and a leprosy bacillus. Estimation of the dN and dS values required use of clusters of orthologous genes. To produce clusters of orthology that incorporated M. leprae and/or M. lepromatosis pseudogenes and corresponding genes from both M. haemophilum and M. tuberculosis, a reciprocal best hit strategy was used with E value 1. 0E-05. The set of pseudogenes/genes of M. lepromatosis and M. leprae were then compared, and they were further aligned to the sets of M. haemophilum and M. tuberculosis. Alignments of genes and pseudogenes to maintain the codon structure, despite the indel mutations of the pseudogenes, were performed with MACSE  (83). This program was also used to align clusters of functional genes in the four species. The program codeML  (84)  was used to obtain the dN and dS values as previously described  (58).
A custom python script was used to run those analyses. To estimate relative pseudogene ages (Rpa), the following formulas were used, for M. lepromatosis (up) and for M. leprae (down): The i (intercept) corresponds to the node of divergence between the leprosy bacilli and M. haemophilum lineages. The values dN ilp_ps and dN ile_ps are the numbers of nonsynonymous substitutions in the branches of a pseudogene from M. lepromatosis and M. leprae, respectively. To be pointed out, the term dN only shows a strict sense for the period of evolution as a gene. Rpa ¼ ðdN The parameters f(dN ih ) and f'(dN ih ) are estimates, for orthologous genes, of the expected dN value for the branch of a leprosy bacillus as a function of the dN value for the branch of M. haemophilum considering that this gene had been evolving as a functional gene for the complete period of time.
2813 for M. lepromatosis and M. leprae, respectively). These are the expected dN values for pseudogenes with Rpa values of 1. The time of divergence in Mya between the lineage of leprosy bacilli and the M. haemophilum lineage (Rpa = 1) was estimated with the RelTime method  (28). All mentioned scripts are available at https:// github. com/DiegoSantos-Garcia/Scripts.
ACKNOWLEDGMENTS This work was supported in part by a  University Cancer Foundation  grant from the  University of Texas MD Anderson Cancer Center  (to X. Y. H. ), by the  National Institutes of Health (NIH)  grant  P30 CA016672  for  MD Anderson's core Sequencing and Microarray Facility and core Bioinformatics Shared Resource , by the  NIH  grant  CA143883  for  MD Anderson 's TCGA Genome Data Analysis Center, and by the  Mary K. Chapman Foundation. D. S. G. was recipient of a contract from the project  ANR Hmicmac 16-CE02-0014. F. J. S. was funded by  Ministerio de Ciencia, Innovación y Universidades (Spain)  ( PGC2018-099344-B-I00 ) and  Generalitat Valenciana  ( Prometeo/2018/A/133 ) and cofinanced by the  European Regional Development Fund (ERDF). We thank  Nipun Mistry  for participation in initial work.
20 ), and a state-of-the-art land surface model, ORCHIDEE-Hillslope (Organizing Carbon and Hydrology In Dynamic Ecosystems-Hillslope)  21, 22. This model partitions the water in each model grid cell, accounting for the  wetter lowlands which result from the convergence of hillslope water flows (see details in Methods). Based on the SM output from ORCHIDEE-Hillslope, the flooded area dynamics of each grid cell is then calculated using a subgrid hydrological model  [23] [24] [25]  calibrated to match satellite observations of inundated areas (Methods). We first perform historical simulations for 2000-2016 with ORCHIDEE-Hillslope forced by climate data from GSWP3-W5E5 (refs. 26, 27 ) and annual land-cover maps based on forest change from China's NFI. The contribution of forest change to wetland change is isolated with two factorial simulations: S0 (simulation with climate change (hereafter called CC), elevated CO 2 concentration (hereafter eCO 2 ), but without forest change) and S1 (simulation with CC, eCO 2 , and forest change).
3d ). On the contrary, in dry regions where more than 80% of annual precipitation is lost to the atmosphere as ET (Supplementary Fig. 2a ), the substantial increase of ET due to afforestation leads to a more negative δSM δA forest (~-0. 04 mm km -2 ) and therefore a higher δA wet δA forest. In the mesic climate zone with a PET/P of 1-2, δA wet δSM is larger than in wet or dry regions (Fig. 3d ), which thus raises the δA wet δA forest in this region. According to the List of Protected Wetlands in China  37  , we identify the basins, at level 6 as classified by the global HydroBASINS database  38  containing National, Provincial, as well as Municipal and County-level wetland conservations (hereafter called BAS N , BAS P , and BAS MC , respectively; Fig. 4a ) and then investigate the potential impacts of continuous tree planting on those protected basins in China.
Methods Forest inventory data. To investigate the forest change in China during the last few decades, we collected data from the second to ninth National Forest Inventory (NFI) released by China's State Forestry Administration (data for the first NFI are not available)  53, 54. At the province scale, the forest inventory in China is carried out every 4-5 years, covering 1973-2018. Forests in the NFI are defined as lands with more than 20% tree cover, including arbor forests, bamboo forests, and shrubs. Bamboo forests and shrubs are under-represented in ORCHIDEE. We, therefore, only focus on arbor forests in this study; they account for ~90% of all forest area and ~80% of the annual increase in forest area in China. To generate a forest map for each year, we combined the spatial information from the 1:1,000,000 Chinese Vegetation Map  55  with the annual forest area change linearly interpolated from the five-year-interval forest inventory data.
We used two observationbased wetland data sources to conduct wetland-related analysis. First, to estimate historical wetland change, we used the satellite-based global inundation product GIEMS-2 (Global Inundation Estimate from Multiple Satellites version 2; ref. 20  ) to estimate the wetland change in China from 2000-2015. By combining passive and active microwaves, along with visible and near-infrared observations, this 0. 25°× 0. 25°product gives monthly estimates of surface water extent, including wetlands, open water, and rice paddies, but potentially excluding large lakes, rivers, and reservoirs. Some small water bodies that cover less than 10% of the grid could be missed in GIEMS-2 due to the relatively coarse resolution and the dense vegetation. The spatial distribution and temporal variation of the global inundation area from GIEMS-2 have been validated with existing independent products such as precipitation and altimeter river height  20.
Due to the small landscape structure of most WNRs, we extracted the basins at level 6 as classified by the global HydroBASINS database  38  intersecting with these WNRs, to investigate the effects of forest change on protected wetlands at a basin scale. Budyko conceptual model. The Budyko framework was proposed five decades ago. It is based on the empirical relationship between annual mean ET/P and PET/P (ref. 31 ). By assuming changes in soil and groundwater storage are negligible over annual time scales, it indicates that most P in a catchment is allocated to Q in wet and energy-limited areas, while most P goes to ET in dry and water-limited areas (Supplementary Fig. 2 ). Despite the simple and "lumped parameter" structure, the intuitive framework has been successfully applied to explain and predict how the terrestrial hydrological cycle has changed up until now  33, [62] [63] [64].
6-2. 2, where the decrease in Q due to 20% forest gain is equivalent to more than 2% of P (Fig. 2 ). ORCHIDEE-Hillslope Simulation. The ORCHIDEE land surface model 21 simulates the terrestrial carbon and hydrological processes and has been widely used for the detection and attribution of the global or regional carbon and hydrological cycles  [65] [66] [67]. The carbon module simulates photosynthesis, litterfall, and soil carbon dynamics, while the hydrological module describes the partitioning of P into ET and Q, and the water redistribution in the 2-m soil column is based on the Richards equation  66, 68, 69. Both the carbon and water processes in ORCHIDEE are very dependent on the vegetation cover, which is described as a mosaic of up to 13 plant functional types (PFT) including bare soil, nine forest types, C3 and C4 grasslands, and croplands for each grid cell.
Acknowledgements The study was supported by the  National Natural Science Foundation of China  (grant numbers  41830643  and  41722101 ), and the  National Key Research and Development Program of China  ( 2016YFC0500203 ). P. C. acknowledges support from the  ANR CLAND convergence institute. We are grateful for the computational resources provided by the  High-performance Computing Platform of Peking University 's supercomputing facility.
Data availability All observation and model data that support the findings of this study are available as follows. The National Forest Inventory data were available from China's State Forestry Administration (http://www. forestry. gov. cn/). The GIEMS-2 data set analyzed during the current study over China from 2000 to 2015 has been deposited on Zenodo (https://doi. org/10. 5281/zenodo. 5750962)  76. The RFW data sets are available at https://doi. pangaea. de/10. 1594/PANGAEA. 892657. The HYDE v3. 2 data set are available at https://easy. dans. knaw. nl/ui/datasets/id/easy-dataset:74467. The historical and future climate data from GSWP3-W5E5 and ISIMIP3b are obtained from https://esg. pik-potsdam. de/search/isimip/. The protected wetland locations in China are obtained from http://www. zrbhq.