dataesr/scientific-paragraphs-categorization huggingface dataset:
2016). As explained in Diesing et al. (  2014 ) and  Li et al. (2016) , the main underlying assumption of this method is that the predictive power of multiple classification trees (the elemental unit of machine learning methods) is higher than that of a single tree. Bootstrapped samples from the training data are used to construct the individual trees in the forest introducing the first element of randomness. In turn, a random subset of the predictor features is used at the node splits throughout the construction of the model. The result is the construction of unique trees. Decisions about the class allocation (labelling) are made on the basis of majority votes of individual trees. After a feature selection procedure, the RF was run growing 501 trees and leaving the parameters as default. The routine was implemented in R (R Development Core Team 2008) using the RandomForest package  (Liaw and Wiener 2002).
To validate the application of an unsupervised classifier, paired-pixel metrics of map agreements were computed after  Foody (2004) ,  Pontius and Millones (2011) , and  Pontius and Santacruz (2014). The R package diffeR was used  (Pontius and Santacruz 2015). Components of allocation are used to derive the agreement between maps at the level of the entire landscape and per category. Quantity and Allocation describe the amount of change that is respectively due to the proportion of categories between reference and test instances and to the amount of spatial mismatch between categories.
Change detection Three types of analysis were performed on the backscatter time series in order to extract trends and patterns of change in substrate classes: pre-, post-and ensemble-methods classification. Ensemble approaches combine supervised and unsupervised classifiers, whereas a pre-classification method focuses on the unclassified data values (similarly to directly relying on spectral bands in satellite imagery). The aim of a post-classification approach is to allocate classlabels to the data values to produce thematic maps.
Ensemble approach classification An ensemble method, combining supervised and unsupervised classifications was also applied. K-means classes (dB ranges) identified in T 7 were used to reclassify the complete dataset for which ground-truth data were not available. From this classified dataset, proportion counts were extracted to observe temporal trends. Prior to transforming the successional backscatter mosaics into classified data, the Within Group Sum of Squared Distances plot was computed independently for each dataset. This ensured that the number of classes in each time series was maintained, and allowed testing. This also serves to test the class discrimination potential of data gathered at 100 and 300 kHz from the EM1002S and EM3002D, respectively. This technique is similar to computing a silhouette plot where the optimal number and size of classes in a dataset becomes visible  (Eleftherakis et al.
Supervised map of the study area Figure  5  shows the visual validation of the sample datasets. This showed an overall good representation of the BS variation in the study area (Fig. 5a ). Mean backscatter values, extracted within 10 m circular buffers at the sample's positions, indicate good separation of the classes (Fig. 5b ) where coarse-hard and fine-soft classes exhibit the highest and lower backscatter values respectively. Similarly, the separation using the bathymetry evidences the distribution of substrate types within different depth zones (i. e. fS on the top of the sandbank, mS + b transiting to the deepest area, and cS + G in the gully; Fig. 5c, d ). The predicted substrate classes' distribution by the Random Forest supervised classification is shown in Fig. 6b. The most important variables selected by the feature selection tool were BS, BS 3 × 3 mean filter, BS Local Moran and bathymetry.
Ensemble approach classification Class proportion counts per survey were extracted from the classified EM3002D dataset (ensemble approach) and are shown in Fig. 9. Temporal trends' and classes' relationships are shown for the entire study area, as well as for the three selected gravel refugia. The fS class appears to be relatively stable across all instances and survey. An inversely correlated relationship is evident between cS + G and mS + b classes. This is also shown in Fig. 10  where the proportion counts per survey are plotted against each other. At the level of the entire study area (Fig. 9d ), and similarly to the pre-classification analysis, this method indicates that the class proportions return to their original state. On the contrary, within the gravel refugia zones (Fig. 9a-c ), the cS + G class experiences a net loss in favour of finer substrate types with no indication of re-establishment to a previous state.
Multibeam backscatter in a monitoring context The basic premise in using MBES backscatter data for seafloor change detection is that changes in substrate cover must result in changes in backscatter values and changes in backscatter due to seafloor cover change must be large with respect to changes caused by other factors (readapted from Singh 1989) such as sea conditions, sensor's intrinsic characteristics, changes in on-board acquisition parameters, vessel speed and direction of survey  (Rattray et al. 2013; Lurton and Lamarche 2015). As such, verification of MBES backscatter stability is critical and should be controlled  (Anderson et al. 2008). In this study, these limitations were mostly overcome as the dataset used was acquired by maintaining rigorous standards of acquisition and processing, including careful attention on environmentally dependent transmission losses (i.
e. by regular control of absorption coefficient). To verify instrumental drift on the medium to long term, the trend in backscatter levels was compared against a time series in backscatter levels at a known reference area (KWGS reference area; Blue Polygon in Fig. 1 ;  Roche et al. 2016 ). As such, average backscatter levels of the RV Belgica EM3002D could be compared from one survey to another during a similar period and allowed obtaining a dataset with temporally consistent dB ranges (yet relative). However, changing environmental factors and seabed conditions may affect the backscatter values also. The effect of biological activity, which is seasonally driven and linked to the spawning and recruitment period of benthic species, is probably the most prominent factor. From literature, it is known that megabenthic zoo-and/or phytobenthic structuring species can be responsible for significant changes in the acoustic signal (e.
Change detection methods The pre-classification analysis of the backscatter time series indicated that within the selected regions of interest, no significant changes in seabed substrate could be detected across the timespan analysed. Since the first dataset was recorded with a former-generation echosounder, which was not cross-calibrated with the EM3002D and using a different frequency range, the values derived could not be directly compared in terms of the range in insonification values. The only evident behaviour in the data was in the barchanoid and gravel gully regions where locally, the mean backscatter level fluctuated around the 1 σ deviation (Fig. 8a, c ). Since the comparison is rather focused on the spatial delineation and areal extent of the substrate classes rather than the intrinsic, physical characteristics of a circumscribed area, the post-classification methods, as adopted similarly to  Rattray et al.
(2013) , did allow comparing data from different echosounders, and acquisition parameters. This was also possible due to the agreement in the number and size of classes discriminated by the two echosounders. The approach revealed information on the behavioural tendencies of certain substrates to undergo a certain change such as the negative trend of the hard substrate class and gain of the finer substrates. In this study, there was a high agreement between supervised and unsupervised models, using quantity and allocation agreement/disagreement metrics, which allowed extending the analysis to the entire time-series dataset. As such, the initially ground-validated information could be fully exploited and extended to the full backscatter dataset. Substrate class proportions over time could be extracted so that global changes could also be accounted for.
This is unlike the pre-classification approach that is limited to selected sub-areas where backscatter levels were extracted from. Therefore, the ensemble approach combined supervised and unsupervised classification algorithms (similarly to  Ierodiaconou et al. 2005 ) which allowed using one ground-truth dataset to train a classification that was subsequently applied to the whole time series. This is a big advantage since sampling of each time series is most often not realistic given survey time and cost restrictions. Here, application of consistent data acquisition and processing allowed the comparability of instant statistical analysis results at various times. Accordingly to recent research  (Li et al. 2016) , it was shown that the Random Forest classifier is a highly valuable tool for seafloor applications, producing accurate models and providing information on the most important feature layers used in the classification.
Conclusion This study highlights the importance of researching approaches and testing tools usable for local-and regionalscale environmental assessments (i. e. for MSFD implementation). A selection of useful methodologies was presented to detect changes in seafloor substrate types. The investigation showed how under specific standardised multibeam backscatter acquisition procedures, the confidence of repeated acoustic observations could be enhanced significantly and how the valuable, but expensive groundtruth information could be propagated from one survey to a time-series dataset via the application of supervised and unsupervised classification routines. The serial backscatter dataset was analysed using techniques developed in the remote sensing terrestrial realm showing that the methodologies are applicable for marine environmental monitoring.
Acknowledgements This paper is a contribution to the projects  MOZ4  ( Flemish Authorities ,  Agency Maritime Services  and Coast, Coast, contract  211. 177 ),INDI67 (BELSPO, contract  BR/143/A2/ INDI67 ), ZAGRI, a federal Belgian programme for the continuous monitoring of sand and gravel extraction, paid from private revenues, and  TILES (BELSPO , contract  BR/121/A2/TILES ). Ship time RV Belgica was provided by  BELSPO  and  RBINS-OD Nature. The  Renard Centre of Marine Geology of Ghent University (RCMG)  is acknowledged for the use of the sediment analysis facilities and instruments. Flanders Marine Institute (VLIZ)  is acknowledged for the use of the video frame. Flanders Hydrography  is acknowledged for the use of background bathymetric data. Two anonymous reviewers are acknowledged for their useful comments and improvements on the initial manuscript.
Microarchitectural parameter measurements Image treatment and parameters measurements were performed by a single operator blinded for the assessment using MicroView software ABA 2. 2 (General Electric, Fairfield, USA). An anisotropic filter was applied to the images. Individual semi-automatic segmentation was performed based on Otsu method  [35]  to distinguish subchondral bone plate from trabecular bone. Subsequently, based on Hounsfield numbers, the region of interest (ROI) was defined manually using a contour-based tool on the weight-bearing area of the medial tibial condyle, as follows and illustrated in Fig. 1. Firstly, the X (medio-lateral) and Y (cranio-caudal) axes were set from the intercondylar area to the lateral edge of the cortical bone and from the intercondylar area to the caudal aspect to the medial condyle, respectively.
Acknowledgements The authors wish to thank  BIOVIVO  from the  Claude Bourgelat Institute (IClB), Lyon, France, and Centre Lago, Vonnas, France  for the animal care and logistical assistance. We would also like to thank the IVTV  ANR-10-EQPX-06-01  for their financial support and  Robert Gurny  for giving us the opportunity to further characterize the tissues from his previous work. Finally, the authors would like to thank  Dr. Joseph Hamm  for his kind review of the paper for English quality.
Exon counting was obtained using the featureCounts function of the Rsubread R/ Bioconductor package, which implements a flexible and powerful counting algorithm  [31]. Between-sample normalization at gene level was performed by the R/ Bioconductor package edgeR  [32]  using the weighted trimmed mean method (TMM)  [33]. Relative scaling factors for the libraries were calculated, and normalized counts per million (CPM) values in edgeR were obtained. Additionally, we used fragments per kilobase of transcript per million mapped reads (FPKM) as a measure that should be invariant to the length of genomic features in RNA-seq. FPKM were calculated using the Cuffdiff 2 algorithm  [34]. As this measure is not recommended for differential expression analysis  [35] , we used it only for correlation analysis.
Data transformations for exploratory analysis Exploratory data analysis was performed on log 2 -transformed values of signal intensity. Two measures for gene expression were used. The first measure was log 2 expression (intensity, counts or CPM) of a genomic feature. Calculated as the mean of log 2 expression of corresponding probe sets, it is an absolute measure that characterizes each sample independently and allows, to some extent, estimating the quantity of the corresponding mRNA. To avoid infinity after log-transformation of RNA-seq data, a small constant offset was added to the measurements. For CPM, we used 0. 5, as proposed by the developers of the edgeR package  [38] , therefore logCPM = log 2 (CPM + 0. 5). For FPKM, which can have very small positive values, the constant offset was calculated as 1% percentile of all non-zero values (0.
005 in our dataset); therefore logFPKM = log 2 (FPKM + 0. 005). The second measure was log 2 fold change (logFC). It provides a relative expression of a gene or exon in tumour versus near-by normal tissue. When calculated for a pair of samples, it removes the part of the transcription signal common to both tissues for the same patient and allows concentrating on the differences between paired tumour and normal adjacent tissues.
Differential expression analysis (DEA) Microarray data were analysed by linear models with empirical Bayes statistics from the limma package of R/ Bioconductor in order to detect differentially expressed genes  [40]. For RNA-seq we selected one of the most used methods -edgeR  [32]  which is based on negative binomial models for counting data. The method was recently reported as adequate for analyses with low number of replicates  [41]. We tried both paired analysis, which accounts for patient effect, and unpaired analysis, which assumes no linkage between tumour and normal samples coming from the same patient. As no improvement was seen with paired DEA, we used unpaired analysis for the genes. Benjamini-Hochberg correction was used to control false discovery rate (FDR) among selected significant features (mRNAs or lncRNAs).
We performed DEA on all available features for the two considered types of molecules separately: protein-coding mRNAs and lncRNAs. Significant protein-coding genes identified in this analysis were used afterwards for the functional enrichment analysis. To measure the similarity between the gene lists we used the Jaccard index  [42] : (N HTA ∩ N RNA-seq ) / (N HTA ∪ N RNA-seq ), where N HTA , N RNA-seq are the numbers of the selected features identified by HTA and RNA-seq platforms. Next, public data for lung squamous cell carcinoma (LUSC RNA-seq v2 dataset) from The Cancer Genome Atlas (TCGA) Research Network (http://cancergenome. nih. gov) was used as a reference dataset. This dataset, generated by highly standardized procedures with thorough quality control, contains RNA-seq data from 51 normal and 502 tumour samples.
Functional analysis Enrichment analysis on significantly differentially expressed genes was performed using Gene Ontology (GO) terms from all the available domains (biological processes, molecular functions and cellular components)  [43]. The main analysis was performed for genes with FDR < 10 -4 in order to get a manageable number of significant genes, compatible with regular practice of functional annotation, and to level the size of the selected gene lists. However, in order to avoid stochastic effects and artefacts selecting a fixed threshold for the platforms, we verified our findings in parallel with 4 thresholds for FDR (10 -2 , 10 -3 , 10 -4 and 10 -5 ) and 4 thresholds for the number of top significant genes (500, 1000, 2000, 4000). The enrichment analysis was repeated on these lists using the topGO package of R/ Bioconductor.
Enrichment of gene ontology categories with significant genes was quantified by Fisher's exact test and resulting p-values adjusted by Benjamini-Hochberg's FDR procedure. Next, we summarized the extensive lists of enriched ontology terms by removing redundant terms using REVIGO tools (http://revigo. irb. hr) [44]  with the semantic similarity measure "Resnik" and dispensability scores of the categories <0. 4. To avoid a tool-related bias, we confirmed our findings using the ReactomePA package of R/Bioconductor  [45]. In addition, in order to reduce the transcript length effect on the results of the enrichment analysis, we applied the goseq package of R/Bioconductor. This package was developed specifically to address the problem of length-related bias in RNA-seq data when performing functional annotation of significant genes  [10].
Analysis of differential exon usage The analysis of the differential exon usage was performed by R/Bioconductor tools widely used in the field: diffSplice method from limma for HTA arrays and DEX-Seq package aimed at RNA-seq data analysis  [36]. Both diffSplice and DEXSeq are based on the splicing index (SI)the difference between exon and gene logFCand provide the significance of SI. Both algorithms were used twice for each platform: for expression of exons and for expression of exon-exon junctions. Following the requirements of DEXSeq, paired analysis was used for exon level. Resulting p-values were corrected by the Benjamini-Hochberg's method. Note that SI based methods can detect exons with zero logFC when the absolute gene logFC is high. To prevent this scenario, we combined FDR and exon logFC selecting significant exons.
Exploratory analysis of the gene expression data Sequencing shows higher variability in the expression data We started the gene expression analysis by mapping Affymetrix probe-sets on Ensembl-defined exons. This increased the compatibility between the platforms in terms of the lists of characterized genes, especially for protein coding mRNAs: almost 100% of these genes were found in common between both platforms (Additional file 1: Figure  S3 ). Regarding lncRNAs, 5855 (93%) were also found in HTA data. Overlap of exon IDs was high as well: 92% for protein coding mRNAs and 90% for lncRNAs. Therefore, the list of common genes was considered representative enough and further exploratory analyses, such as principal component, variance and correlation analysis, were performed on this list. Initially, principal component analysis (PCA) of the expression data showed a strong platform effect (Fig.
The part of variability that could not be explained by these factors, or "within-group" variability, was represented by residuals. For both types of RNA, the highest fraction of variability detected by HTA arrays was associated with the tissue state (Fig. 1d, e ). The patient effect was less pronounced than the tissue effect and comparable between the platforms. Importantly, the highest residuals were observed in RNA-seq data, suggesting higher levels of stochastic noise in this dataset. To prove this, we investigated the variability between biological replicates in normal and tumour tissues with respect to the average expression level (Additional file 1: Figure  S4 ). For both tissues, higher variability was observed in RNA-seq data, especially for low expressed genes (Additional file 1: Figure  S4a, d ).
Limma, which uses a more general normal model for the signal, tends to    3  Differentially expressed genes identified by the platforms. Evolution of the number of significant genes identified with variable FDR thresholds (a), using edgeR and limma with voom correction for analysis of RNA-seq data, and using limma for HTA data. Solid lines show unpaired analyses, while dotted lines show analyses paired by patient. Differentially expressed protein coding mRNAs (b) and lncRNAs (c) were obtained by unpaired differential expression analysis using edgeR for RNA-seq and limma for HTA (FDR < 0. 01) and represented as proportional Euler-Venn diagrams. The lists of differentially expressed genes were confirmed by the top 25% significant genes detected in the LUSC-TCGA dataset: 4569 protein coding genes (FDR < 10 -18 ) and 111 lncRNAs (FDR < 10 -8 ) were used.
Evolution of Jaccard index for coding mRNAs with variable FDR thresholds (d) between the two platforms (violet) shows a monotonic behaviour. Similarity between the TCGA validation gene list and each of the platforms-RNA-seq (red) and HTA (blue) showed a slight outperformance of HTA (marked by an arrow) penalize additional factors strongly. Taking these facts into account we decided to continue the study using the unpaired analysis. Interestingly, RNA-seq data analysed by edgeR reported fewer genes than HTA data analysed by limma for moderate FDR values (0. 05-10 -4 ). For the standard threshold (FDR < 0. 01), 6173 protein-coding mRNAs were detected from HTA data while only 4777 were obtained from RNA-seq data, with 3683 genes in common (Jaccard index of 0. 507), as shown in Fig. 3b. Likewise, 1219 lncRNAs were found to be significant on the HTA platform against 892 in the RNA-seq experiment (Fig.
3c ), but with much lower similarity (0. 217). However, genes with a more stringent threshold (FDR < 10 -4 ) were observed mainly for edgeR with its negative binomial model. When the analysis assumed the normal model, the reported number of significant genes fell even faster with the decrease of FDR for RNA-seq than for HTA arrays (Fig. 3a ). As the two curves representing the number of significant genes in HTA and RNA-seq crossed around FDR < 10 -4 (Fig. 3a ), this threshold was used to select the similar-size sets of significant genes for the functional annotation. In order to verify the obtained lists of differentially expressed genes, we compared them with external results obtained on squamous cell carcinoma samples from TCGA (LUSC dataset). This dataset can be used as a reference due to its high quality and the large size of the patient cohort used.
As mentioned in Methods, 25% of the most significant genes in the TCGA dataset were selected creating a validation gene set of similar size to the set of significant genes in our experiments. For protein-coding mRNAs the validation set was composed of 4569 significant genes with an FDR < 10 -18. The number of genes identified in our analysis and confirmed by the TCGA dataset are summarized in Fig. 3b, c. HTA showed a higher number of confirmed mRNAs, but proportionally the confirmation level varied only slightly: 37% for HTA, and 39% for RNA-seq. Common mRNAs showed the highest confirmation level (43%), while the uniquely identified genes had lower confirmation rates of 23% and 26% for sequencing and microarrays, respectively. We also investigated the evolution of overlap between platforms and the confirmation rate within each platform (Fig.