dataesr/scientific-paragraphs-categorization huggingface dataset:
This led to a high accuracy from the initial training epochs on both Training and Validation sets. Different learning rates were also tested for each dataset to try to obtain high accuracy on the Validation set with a reduced number of epochs. Tables  2, 3 , 4 and 5 report the results of the experiments using 10 epochs, the accuracies on the test set are highlighted in bold. From these tables we can see that using datasets of different sizes and different learning rates led to very similar training accuracy. Dataset ds4 (i. e. 5k for Neutral class and 5k for Selection class) is the one that produces the lowest accuracy, but still in line with the other results, since images in the datasets are still very similar. It can also be observed that dataset ds3 provides a good balance between number of samples, accuracy and training time.
Note that, the value of the learning rate that leads to higher accuracy is 0. 001. In order to increase the accuracy, we first increased the number of training epochs. Further experiments with 50 epochs were performed and led to a higher accuracy. As depicted in Fig. 8 , the training accuracy and the training loss of the model trained on the various datasets during 50 epochs are very similar. Observe that dataset ds4 initially deviates from the others but converges similarly when the maximum number of epochs is reached. 6  and 7  respectively, the similarity of the training, validation and test accuracies can be observed. This clearly highlights the fact that the implemented model is generalizing correctly and does not exhibit overfitting. Further experiments were performed but they did not lead to remarkable improvements.
Application to real data As a proof of principle, we tested our CNN using a real dataset. We chose a clear example of a selective sweep in the SLC24A5 gene (variant rs1426654, Chromosome 15:48134287; Ensembl. org), which was identified as responsible for skin pigmentation loss in human populations of European ancestry, see  [32, 33]. Using the variant data from the 1000 Genome Project, we selected the genomic region from position 48000000 to 48500000 on chromosome 15 for 24 individuals sampled in Tuscany (Italy). As the signature left by the selective sweep in this gene spans more than 100 kb, we adjusted our CNN by training it with a new set of simulated windows of 10,000 base pairs generated as in the experiment before but setting an effective population size of 10,000 and a mutation rate of 2E10-8.
The training set included 3910 images (1990 as selection and 1920 as neutral), whereas both the training and the validation set included 500 images of each type (dataset ds5). We trained the CNN in Fig. 7  using the Adam optimizer with a learning rate of 0. 001. The training was done with early stopping, i. e. , the training stops if there is no improvement in the validation accuracy after three epochs. The trend of the accuracy and loss on both training and validation is shown in Figs. 12 and 13 respectively. Note that, contrary to the training described in the Experiments section in which after the first few batches of the first epoch the model is already able to distinguish neutral and selection, the model struggles in the first 7 epochs. Only from the 8th epochs onwards the model starts extracting relevant features necessary to distinguish neutral from selected images.
This is probably due to fact that images of 10000x48 contain more information than images of 1000x48. Table  9  recaps the accuracies and losses on the training, validation, and test sets for simulated genomic images. We converted the real genomic data into 50 images of 10000 rows and 48 columns using a custom python script (vcf2raster. py) and categorized the windows between position 48100000 and 48220000 as selection (12 windows) while the remaining 38 windows before and after those positions were labeled as neutral (38 windows). Besides the window including the causative variant rs1426654, the other genomic windows were labeled as under selection if nucleotide diversity was below 0. 00003 substitutions per site. This labeling implies that the identification of the regions flanking the causative variant and characterized by a similar genetic variation pattern are also of interest.
A denote such matrix: in the first row, A 11 (called True Positives TP) counts the number of Neutral examples correctly classified while A 12 (called False Negatives FN) denotes the number of Neutral examples wrongly classified as Selection. Similarly, in the second row A 21 (False Positives FP) and A 22 (True Negatives TN) denote the number of Selection examples wrongly and correctly classified respectively. Figure 11 presents the confusion matrix for dataset ds3 using 50 epochs. The results are clearly in line with the values shown in Tables.
Let A denote such matrix: in the first row, A 11 (called True Positives TP) counts the number of Neutral examples correctly classified while A 12 (called False Negatives FN) denotes the number of Neutral examples wrongly classified as Selection. Similarly, in the second row A 21 (False Positives FP) and A 22 (True Negatives TN) denote the number of Selection examples wrongly and correctly classified respectively. Figure 11 presents the confusion matrix for dataset ds3 using 50 epochs. The results are clearly in line with the values shown in Tables.
1. Allele size was quantified using an Applied Biosystems 3130 Genetic Analyzer. To score alleles, we used the package Fragman  (Covarrubias-Pazaran et al. , 2016)  in R 3. 5. 1 software (R Development Core Team. R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2018. http://www. r-project. org). To reduce genotyping errors, we used a comparative multiple-tube approach to determine multilocus genotypes as recommended by  Taberlet et al. (1996). We used the allelematch package  (Galpern et al. , 2012)  in R software to perform multiple comparisons between genetic profiles.
Genetic clusters and genetic diversity We used a Bayesian approach implemented in STRUCTURE 2. 3 software  (Pritchard et al. 2000)  to determine if the Gehuaqing/Xiangguqing subpopulation was sub-structured into K genetic clusters, based on the multilocus genotypes of the 74 individuals sampled (see Results). Twenty independent runs for each value of K (1e10) were performed with 1,000,000 Markov chain Monte Carlo (MCMC) repetitions and a 100,000 burn-in period using no prior information and assuming correlated allele frequencies and admixture. K was identified using two methodologies: the maximum values of ln P(D) (the log probability of the data for a given K) returned by the software and DK based on the rate of change in the likelihood of K  (Evanno et al. 2005). With this model, STRUCTURE estimates the posterior probability q i , which is the proportion of each individual genotype that is derived from each of the K genetic clusters  (Pritchard et al.
We constructed a minimum spanning network of mtDNA haplotypes using an infinite site model of DNA sequences  (Templeton et al. 1992). Statistical analyses of mtDNA were performed using the pegas R package  (Paradis 2010). The mitochondrial DNA sequences generated in this study were aligned with sequences from  Liu et al. (2007)  using the Decipher algorithm as implemented in the Decipher R package  (Wright et al. 2012). A phylogenetic tree was constructed by testing which nucleotide substitution model best fit to our sequence data by using the phangorn package in R software  (Schliep, 2011). We estimated the allelic richness per subgroup using a rarefaction procedure corrected by the number of individuals sampled per subgroup. Genetic variation was estimated over all loci within each subgroup by assessing the observed (H O ) and expected (H E ) heterozygosities.
Deviations from Hardy-Weinberg equilibrium (HWE) and from linkage equilibrium for each locus within each subgroup were estimated using an MCMC-based test. All these calculations were performed using the R packages adegenet  (Jombart 2008) , hierfstat  (Goudet 2005) , and pegas. The genetic differentiation between subgroups was quantified by computing the  Weir and Cockerham (1984)  estimator F ST with the R package hierfstat. Individual heterozygosity was assessed as the proportion of the 10 microsatellite loci for which a given individual had two distinct alleles. Individual heterozygosity was then compared between the fed and wild subgroups using Fisher's test after analysis of variance (ANOVA).
Relatedness We compared relatedness between fed individuals and between wild individuals using the R package related  (Pew et al. , 2015). First, we compared the performance of six relatedness estimators using simulations to select the best estimator given our genotype data. We used the moment estimators of  Li et al. (1993) ,  Queller and Goodnight (1989) ,  Lynch and Ritland (1999) , and  Wang (2002) , and the likelihood estimators of  Milligan (2003)  and  Wang (2007). We simulated data sets of 100 pairs for four relatedness categories (parenteoffspring [PO], full-sibling [FS], half-sibling [HS], and unrelated [UR]), with allele frequencies obtained from values observed in our field dataset. We then assessed relatedness using the six estimators listed above and compared the relatedness estimates from simulations with the expected relatedness values (PO ¼ 0.
5, FS ¼ 0. 5, HS ¼ 0. 25, and UR ¼ 0. 0) by assessing Pearson's correlation coefficients. We selected the estimator with the highest correlation coefficient for the subsequent analyses. We then assessed relatedness for all individuals in our dataset, and we tested whether individuals within the wild and fed subgroups were more related than expected if these subgroups simply represented random clusters of individuals  (Pew et al. , 2015). We performed a permutation test using 100 iterations, to compare relatedness between the fed and wild subgroups.
Subgroup sizes Estimating population or group size based on non-invasive DNA capture-mark-recapture data is increasingly used in wildlife monitoring because of its high level of reliability  (Petit and Valiere 2006). However, all the methods developed to assess population size are based on key assumptions that are not fully met in the Gehuaqing/Xiangguqing Yunnan snubnosed monkey subpopulation, especially because the demography of the fed subgroup is artificially driven by human intervention. Reserve officers are very familiar with the composition of the fed subgroup and can recognize individuals based on their appearance. We thus estimated subgroup size for the fed and wild subgroups however we considered that the estimate of the fed subgroup size made by reserve officers is likely the most reliable.
The DNA-based mark-recapture capwire package for R software developed by  Pennell et al. (2013) , based on the models developed by  Miller et al. (2005) , was used to estimate the subgroup sizes. Designed for the estimation of population size with non-invasive sampling of wild populations, the methods developed in the capwire package handle data inferred from multiple observations of the same individuals within a sampling session. The even capture probability model (ECM) assumes that all individuals have an equal probability of being captured, whereas the two innate rates model (TIRM) assumes that the population contains a mixture of individuals with distinct capture probabilities. We first computed maximum likelihood estimates of group size using ECM and TIRM models, and we selected the most suitable model by using a likelihood ratio test (LRT).
Simulation study We explored the predicted effects of gene flow from the wild to fed subgroup on the genetic diversity of the fed subgroup over a period of ten years, under the scenario of human intervention by individual translocations. We implemented a simulation routine in R software, based on yearly estimates of genetic diversity descriptors (F ST ; expected and observed heterozygosities H E and H O , census subgroup size N and effective subgroup size N e ; and mean relatedness). Little data are available on the demographic parameters of snub-nosed monkeys. Consequently, we used the few known parameters  (Cui et al. , 2008; Kirkpatrick et al. , 1998; Xia et al. , 2019; Liu et al. 2007; Ren et al. , 2011)  and complementarily a life table for populations of blue monkeys (Cercopithecus mitis)  (Bronikowski et al.
, 1998; Xia et al. , 2019). At time 0, a given number of wild mature females chosen arbitrarily (m) were virtually introduced into the fed subgroup and considered fed individuals for the rest of the simulation. For each simulation, breeding males were randomly selected among mature fed individuals based on the value of one male for four females in OMUs, as reported in the literature  (Cui et al. , 2008). These males conserved their status as breeding males and had a maximal survival rate for the entire simulation. Each virtual year of the simulation (between 1 and 10), we applied the mortality rate based on survival data as described above. Breeding females were selected among females at least five years of age with the probability defined above. We then simulated mating between breeding males and mature females in the fed subgroup, supposing (i) random mating, and (ii) the production of one infant per female.
Genetic clusters and genetic diversity The 20 independent runs of the Bayesian analysis implemented in STRUCTURE software showed maximized average ln P(D) and DK values when K ¼ 2 (Fig. 3 ). Thus, we concluded that K ¼ 2 was the most likely number of genetic clusters within the Gehuaqing/Xiangguqing subpopulation. Then, we compared our assignment of each individual based on field observations with the assignment based on the STRUCTURE posterior probability q (Table  1 ). The mean posterior probability was 0. 187 for cluster 1 and 0. 813 for cluster 2 in fed individuals and 0. 743 for cluster 1 and 0. 257 for cluster 2 in wild individuals. Depending on the threshold chosen, the assignment based on field observations coincided variably with the assignment based on the STRUCTURE posterior probability q: 59.
Fig. 4b  shows that this aggregation was mainly due to the low genetic distance between fed individuals. Wild individuals for which samples were collected geographically close to the fed subgroup were genetically similar to fed individuals, and genetic distance increased with geographical distance. All except one of the 74 individuals sequenced for the CR HVI were characterized by mtDNA haplotypes previously published  (Liu et al. , 2007). Nine of the ten mtDNA haplotypes found in individuals sampled in the current study were described by  Liu et al. (2007)  in the central and southern parts of the distribution range of the Yunnan snub-nosed monkey (M12-M24). However, we noticed a mismatch between the name of mtDNA haplotypes in this reference  (Liu et al. , 2007)  and information available on the GenBank website (see Supplementary information S3).
One haplotype (TCW01) has never been recorded in the literature and was deposited on GenBank (accession number: MN171525). Using a minimum spanning network and a maximum likelihood tree, we observed the two mtDNA haplogroups previously described for the Yunnan snub-nosed monkey (A and B; Figs. 5  and 6 ), although clade B was infrequent. The most abundant haplotypes (M16-M18 and M20) coincided with  Liu et al. 's observations, except for M24 (named M26 in Liu et al. , (2007) ) which has not yet been observed in the Gehuaqing/XianGuQing subpopulation (G09 in  Liu et al. , 2007) , but rather in southwestern groups. Haplotype diversity was similar between fed (h ¼ 0. 82 ± 0. 01) and wild (h ¼ 0. 82 ± 0. 01) subgroups, in accordance with values observed by  Liu et al. , (2007). Nucleotide diversity was low for both the fed (p ¼ 0.
Acknowledgments The authors sincerely thank all the guards of the  National Nature Reserve of BaiMa XueShan  for their invaluable help and advice. Sequencing analyses were performed thanks to the genomic platform of the  University of Franche-Comt e and to the firm Biomnigene (Besançon, France). This work was carried out with the support of GDRI EHEDE (http://gdri-ehede. univ-.
Deep learning methods are widely used for medical applications to assist medical doctors in their daily routine. While performances reach expert's level, interpretability (highlighting how and what a trained model learned and why it makes a specific decision) is the next important challenge that deep learning methods need to answer to be fully integrated in the medical field. In this paper, we address the question of interpretability in the context of whole slide images (WSI) classification with the formalization of the design of WSI classification architectures and propose a piece-wise interpretability approach, relying on gradient-based methods, feature visualization and multiple instance learning context. After training two WSI classification architectures on Camelyon-16 WSI dataset, highlighting discriminative features learned, and validating our approach with pathologists, we propose a novel manner of computing interpretability slide-level heat-maps, based on the extracted features, that improves tile-level classification performances.
Introduction Since their successful application for image classification  [1]  on ImageNet  [2] , deep learning methods, especially Convolutional Neural Networks (CNN), have been extensively used and adapted to tackle efficiently a wide range of health issues  [3, 4]. Along with these new methods, the recent emergence of Whole Slide Imaging (WSI), microscopy slides digitized at a high resolution, represents a real opportunity for the development of efficient Computer-Aided Diagnosis (CAD) tools to assist pathologists in their work. Indeed, over the last three years, notably due to the WSI publicly available datasets, such as Camelyon-16  [5]  and TCGA  [6] , and in spite of the very large size of these images (generally around 10 giga pixels per slide), deep learning architectures for WSI classification have been developed and proved to be really efficient.
It can be performed using Otsu thresholding  [8] , color space transformation (e. g. , thresholding on the saturation channel of HSV color space in  [9] ), or semantic segmentation using a U-Net  [10]  (e. g. , in  [11] ). Once tissue is detected, another step that consists in stain normalization is generally carried out. The motivation behind this step is to improve the transferability to other datasets that come from a different hospital that might use a different scanner or staining to create their WSIs  [12]. Color deconvolution  [13]  is the most popular approach  [14] [15] [16]  but other methods such as color transfer  [17]  or more recently cycle-GAN [18] can be applied as well. Finally, tiling and feature extraction (i. e. , encoding) are performed to enter the MIL context and reduce the input size.
This approach reaches an AUC of 0. 858 on Camelyon-16 and 0. 915 on TCGA-lung (subset of TCGA dataset related to lung cancer). In  [26] , an attention module  [33]  is used instead of a min-max layer. AUC values of 0. 775 for a breast cancer dataset and 0. 968 for a colon cancer dataset were reported. Recently, more works on large datasets proposed architectures that follow the same design  [22, 24, 29, 31]. Heat-maps based on intermediate scores computed in these architectures are what we call architectural explainability that results from prior knowledge on WSI problems that is introduced by design in the architecture. They are of great interest and have proved to be really efficient to the point of being able to spot cancerous lesions that had been missed by experts (in  [22] ). However explanations are relying on a single "medical" score which might limit the interpretability regarding complex tissue structures that can be found on these slides.
Architectures We validate our approach on two WSI classification trained architectures: CHOWDER and Attention-based classification. CHOWDER  [21]  uses a 1 × 1 convolution layer to turn each tile descriptor into a single tile score. These scores are then aggregated using a min-max layer, that keeps the top-R and bottom-R scores (e. g. , empirically R = 5 gives the best results), to give a slide descriptor (M = 2 × R). Attention-based architecture  [26]  uses an attention module (two 1 × 1 convolution layers with respectively 128 and 1 channels and a softmax layer) to compute competitive and normalized (sum to 1) tile scores from tile descriptors. Then, the slide descriptor is computed as the weighted (by tile scores) sum of tile descriptors (M = N). Note that in our experiments the feature extractor is a ResNet-50  [19]  (N = 2048) trained on ImageNet and is part of the preprocessing thus is not fine-tuned.
Datasets We validate our approach using Camelyon-16 dataset that contains 345 WSI divided into 209 "normal" cases and 136 "tumor" cases. This dataset contains slides digitized at 40× magnification from which we perform sample detection using Otsu thresholding  [8]  on a thumbnail of the slide downscaled by a factor 32 and keeping tiles that contain at least 50% of foreground pixels w. r. t. Otsu segmentation. Then, we extract, with regard to a nonoverlapping grid, 224 × 224 pixels tiles at 20× magnification without stain normalization. Then, we pre-compute, for each tile, 2048-tile descriptors using a ResNet-50 model trained on ImageNet as it is done in  [9, [20] [21] [22] [23] [24]. 216 slides are used to train our models while 129 slides form the test set to evaluate performances of the different trained models.
Slide-Level Explanations: Heat-Maps Coherence between patches exposed for a better interpretability led us to think about another way to present features to pathologists. Indeed, since tissues have a coherent and somehow organized structure, a relevant feature for histological problems would be activated in a coherent and somehow organized way over slides. Thus, along with patch-based visualization, we propose to access features activation heat-maps H c,i over slides, as presented in Section 3. Figure  7  illustrates qualitative results. Quantitatively, we report a tile-level localization AUC of 0. 884 for CHOWDER model and 0. 739 for Attention-based model, using featurebased heat-map values (that are the average normalized feature activation over all features identified for the"tumor" class, see H c,i,j computation in Section 3) as a "tumor" prediction score and using lesion annotations provided by Camelyon-16 dataset to get the groundtruth label per tile.
Both AUCs are significantly high, which validates our approach of identifying features that are relevant and of computing heat-maps for interpretation and explanation. Note that the AUC computed using tile scores is 0. 684 for CHOWDER model and 0. 421 for Attention-based model (see Table  1 ). We can also note that there is a gap in interpretability between CHOWDER model and Attention-based model while classification performances are similar (AUC of 0. 82 for the CHOWDER model and 0. 83 for the Attention-based model). This gap can be explained by the fact that, in the context of Camelyon-16, identifying one tumorous tile is enough to label a slide as "tumor", so implicit tile classification does not need to be exhaustive to provide meaningful information to the slide level decision module, however if so interpretability will decrease.
Approach This study about the impact of the number of selected features on heat-map quality shows that three kinds of features stand out for Camelyon-16 dataset among Imagenet features: features of interest that activate homogeneously mostly in tumorous regions, adversarial features that activate homogeneously not only in tumorous regions, and unrelated features that either activate non homogeneously or almost do not activate over slides. Figure  11  illustrates the difference between features of interest and adversarial features. It also gives another way to think about the third question we put forward ("Are these features of tile descriptors relevant medically and representative of histopathological information?") by introducing a manner to measure the potential transfer of each individual feature to a given histopathological problem.