dataesr/scientific-paragraphs-categorization huggingface dataset:
The summed proportion for the nonlocal (EHG + Pınarbaşı) ancestry component in Sicily LM HGs ranges between 17±2% and 33±4% (Fig. 2B , Fig. 5A , Data file 1). Interestingly, whereas a Near Eastern-related ancestry component is discerned in many of the Mesolithic HGs from the Iron Gates and other areas in the Balkan, Baltic or Scandinavia, this is not the case for any of the Villabruna cluster individuals. In contrast to previous statements that the Villabruna cluster individuals form a genetically homogenous group  (16, 25) , here the individuals show rather diverse ancestry profiles with various combinations of Sicily EM, EHG and GoyetQ2-related ancestry. Since the Sicily LM HG gene pool contains the distinct Near Eastern-related ancestry but not the GoyetQ2-related ancestry, it is unlikely that the diversity of this group originated solely from genetic drift from these Villabruna cluster individuals.
Expanding early farmers replace local HGs in Sicily Recent ancient human DNA studies have shown that the genetic variation in Early Neolithic groups from central and southwestern Europe is a subset from that found in early farmers from Barcin in northwestern Anatolia and Revenia in northern Greece  (30, (34) (35) (36) (37) 39). Interestingly, early farmers from Diros in Peloponnese Greece might harbour an ancestry component that places them outside of the genetic diversity represented by those from Anatolia Barcin  (17). In PCA (Fig. 1C , Data file 1) the Sicilian early farmers (Sicily EN) plot closest to Early Neolithic groups from the Balkan, Serbia, Hungary, Greece and Anatolia, but not from Iberia. With the f3-outgroup statistic f3(Mbuti; Sicily EN, X) we aimed to determine the best genetic proxies for the overall gene pool related to Sicily EN.
The excess genetic attraction for southeastern European HGs might hence be driven by the farmer, rather than HG, ancestry component in Sicily EN. To further investigate what combination of early farmer and HG ancestry fits the Sicily EN gene pool best, and their respective admixture proportions, we used qpWave (PWave) and qpAdm-based (PAdm) ancestry models. We required a test result to be more extreme (larger) than P = 0. 1 in order not to reject the null-hypothesis of a full ancestry fit (full rank). We found that the ancestry for approximately half of the Sicilian early farmers can be fitted as entirely early farmer ancestry as found in Anatolia EN Barcin, Greece EN Peloponnese, Croatia EN Cardial and Impressa, Hungary EN Koros or Germany LBK (Fig. 4 , Data file 7). However, these early farmer sources were rejected as sole ancestry sources to the combined Sicily EN gene pool (max.
PWave ≤ 0. 011 for Croatia EN Cardial, Data file 7). We could improve the model fit to the Sicily EN gene pool by adding ~4-9% ancestry from Sicily LM or Sicily EM HGs as a second source to early farmer ancestry from either Anatolia EN Barcin (PAdm ≥ 0. 016), Greece EN Peloponnese (PAdm ≥ 0. 059), Croatia EN Cardial (PAdm ≥ 0. 056) or Hungary EN Koros (PAdm ≥ 0. 300) (Table  1 , Data file 7). Notably, three-way mixture models, in which the ancestry from Anatolia EN Barcin was combined with that from early farmers from Ganj Dareh as a proxy for the early farmer ancestry, improved the fit to the Sicily EN gene pool even further (PAdm ≥ 0. 335, Table  1 ). We can therefore not exclude the possibility that the early farmer ancestry in Sicily EN falls partly outside of the broader genetic diversity of the early farmers from Europe and Anatolia Barcin (Supplementary Section S6).
Materials and Methods aDNA analysis. All pre-amplification laboratory work was performed in dedicated clean rooms (82) at the Max Planck Institute (MPI) for the Science of Human history (SHH) in Jena and MPI for Evolutionary Anthropology (EVA) in Leipzig, Germany. At the MPI-SHH the individuals were sampled for bone or tooth powder, originating from various skeletal elements (e. g. petrous, molars, teeth, humerus, phalange, tibia, see Data file 1). The outer layer of the skeletal elements was removed with high-pressured powdered aluminium oxide in a sandblasting instrument, and the element was irradiated with ultraviolet (UV) light for 15 minutes on all sides. The elements were then sampled using various strategies, including grinding with mortar and pestle or cutting and followed by drilling into denser regions (Data file 1).
7. 12) aln and samse programs (90) with a lenient stringency parameter of '-n 0. 01' that allows more mismatches, and '-l 16500' to disable seeding. We excluded reads with Phred-scaled mapping quality (MAPQ) <25. Duplicate reads, identified by having identical strand orientation, start and end positions, were removed using DeDup v. 0. 12. 1  (88). aDNA authentication and quality control. We assessed the authenticity and contamination levels in our ancient DNA libraries (unmerged and merged per-individual) in several ways. First, we checked the cytosine deamination rates at the end of the reads (91) using DamageProfiler v0. 3 (https://github. com/apeltzer/DamageProfiler). After merging the libraries for each individual, we observed 21-52% C>T mismatch rates at the first base in the terminal nucleotide at the 5'-end, an observation that is compatible with the presence of authentic ancient DNA molecules.
Y-chromosome haplogroup determination. To determine the Y chromosome haplogroup for genetic males we used the yHaplo program ((101), available via: https://github. com/23andMe/yhaplo). We based our haplogroup assignment on 13,581 strand-unambiguous ancestry informative SNPs from the ISOGG (International Society of Genetic Genealogy) data set. We called genotypes for these SNP sites by randomly choosing one allele with 'pileupCaller' (https://github. com/stschiff/sequenceTools). Using an in-house script (Choongwon Jeong) the genotypes were converted to an input format for yHaplo. In ancient genomes, missing data for key diagnostic sites may halt the automated search before the most derived haplogroup is reached. We therefore manually checked the coverage for informative SNPs further downstream for the haplogroup that was assigned to each individual.
Acknowledgements. For helpful comments we thank Stephanie Eisenmann ,  Maite Rivollat ,  Thiseas Lamnidis  and other members of the  Department of Archaeogenetics , and  Barbara Pavlek  of the Minds & Traditions research group of the  Max Planck Institute for the Science of Human History. We are grateful for comments on the manuscript from  Didier Binder  from the  French National Center for Scientific Research and Detlef Gronenborn  from the  Leibniz Research Institute for Archaeology. We thank  David Reich ,  Shop Mallick  and  Ian Mathieson  for access to unpublished data. We thank  Sven Steinbrenner  of the  Department of Human Evolution from the Max Planck Institute for Evolutionary Anthropology  for undertaking the stable isotope analysis. Funding. The  Max Planck Society  financed the genetic, isotopic and radiocarbon analyses.
Author contributions. M. Ma, W. H, and J. K conceived the study. M. Ma ,  A. T, M. P , S. T, C. C, V. S and R. dS provided the ancient human remains and input for the archaeological interpretation. M. vdL, C. F, S. N and L. K performed laboratory work with the help of F. A, M. B, R. R, R. S, A. W, G. B and M. Me. M. vdL conducted the population genetic analyses with the help of K. P, C. J and W. H. S. T performed the AMS radiocarbon dating analysis, and  M. Ma  the isotope analysis. M. vdL,  M. Ma , W. H, V. V-M, C. P, S. S, C. J, K. P and J. K. wrote the paper with input from all co-authors. Competing interests. The authors declare no competing interests. Corresponding authors. Correspondence to M. vdL,  M. Ma  or J. K. Data and materials availability. Genomic data (BAM format) are available through the  Sequence Read Archive  (accession number X) and consensus mitogenome sequences (FASTA format) in GenBank (accession numbers X to X).
In addition, UZZ79 and UZZ81 showed an elevated PMMR indicative of a kinship relation (Data file 4). 10 (19) and Schmutzi v1. 0 (94) (Data file 1). Before running Schmutzi, we realigned the reads to the rCRS using CircularMapper v1. 93. 4 filtering with MAPQ < 30. Dataset. For genotyping we extracted reads with high mapping quality (MAPQ ≥ 37) to the autosomes using samtools v1. 3. The DNA damage plots indicated that misincorporations could extend up to 10 bp from the read termini in non-UDG treated and up to 3bp in UDG-half treated libraries. We hence clipped the reads accordingly, thereby removing G>A transitions from the terminal read ends in ds-libraries and C>T transitions in both ss-and ds-libraries. For each individual, we randomly chose a single base per SNP site as a pseudo-haploid genotype with our custom program 'pileupCaller' (https://github.
HP1 proteins are involved in the regulation of specific gene expression programs To investigate the role of HP1 in the regulation of gene expression in liver, we performed an unbiased RNA-seq transcriptomic analysis of libraries prepared from 7 weekold control and HP1-TKO liver RNA. We found that 1215 genes were differentially expressed (730 up-regulated and 485 down-regulated) between control and HP1-TKO livers (with a 1. 5fold threshold difference and an adjusted P ≤0. 05) (Fig. 4A  and supplementary Table  1 ). Analysis of differentially expressed genes (thereafter called HP1-dependent genes) using David Gene Ontology (https://david. ncifcrf. gov/) and Gene Set Enrichment Analysis (GSEA; http://software. broadinstitute. org/gsea/index. jsp) softwares revealed that several biological processes were significantly affected in HP1-TKO livers.
To identify genes commonly deregulated in HP1-TKO and HP1αβ-liverKO livers that might explain the tumor protective role of HP1, we performed a RNAseq analysis in HP1αβ-liverKO (GSE84734). Only 18 genes (13-up and 5 down-regulated) were similarly deregulated in HP1αβ-liverKO and HP1-TKO livers (supplementary table  7 ). Interestingly, on the 13 up-regulated genes, 5 belonged to the family of transcriptional repressors KRAB-Zinc Finger Proteins (KRAB-ZFP). We thus analyzed the expression of several of these krabzfp genes by RT-qPCR in HP1αβ-liverKO, HP1αγ-liverKO and control livers. As illustrated in Figure  4D , zfp951 and zfp992 were up-regulated exclusively in HP1αβ-liverKO livers, zfp345, exclusively in HP1αγ-liverKO livers whereas zfp984 and 5730507C01Rik were up-regulated in both HP1αβ-liverKO and HP1αγ-liverKO livers.
HP1 loss lead to reactivation of specific retrotransposons and over-expression of neighboring genes KRAB-ZFP are transcriptional repressors known to repress themselves and retrotransposons of the endogenous retroviruses (ERV) family  33, 34. To determine whether the HP1-dependent deregulation of KRAB-ZFP was associated with altered ERV expression in liver, we investigated the expression of DNA repeats in our RNA-seq dataset. Coordinates of all annotated DNA repeats of the RepeatMasker database (mm10 assembly) were aligned against the RNA-seq reads and only those that could be assigned unambiguously to a specific genomic locus were analyzed. In total, 846 repeats were deregulated in HP1-TKO compared with control livers with 71. 3% being up-regulated and 28. 7% down-regulated (Fig. 5A  & Supplementary Table  8 ).
elegans  42. However, it is important to note that liver has very specific properties, being mostly quiescent throughout life but able to regenerate upon stress essentially through the re-entry of quiescent and fully differentiated hepatocytes into cell cycle rather than via stem cell proliferation, as it is the case in other tissues  43, 44. These specific properties of hepatocytes could rely on a peculiar loose chromatin organization that might be less sensitive to the loss of HP1 as compared to other cell types for cell viability but more prone to cell transformation as we observed in BMEL cells. In support of this hypothesis, it is noteworthy that, among the 17 human cancers tested in the database from Human Protein Atlas (https://www. proteinatlas. org and Fig. S3 ), liver cancer is the one expressing the lowest levels of cbx1, 3 and 5.
Statistics and reproducibility. The Microsoft Excel or biostatgv (https://biostatgv. sentiweb. fr) softwares were used for statistical analyses; statistical tests, number of independent experiments, and P-values are listed in the individual figure legends. All experiments were repeated at least twice unless otherwise stated. images of paraffin-embedded liver tissue sections from 7-week-old control (Ctl) and HP1-TKO (TKO) mice stained with DAPI (63x magnification). To select mostly hepatocytes, only the largest nuclei with a size comprised between 70 and 150 µm 2 and with a circular shape were selected for this analysis. 2D sections of nuclei were divided in four concentric areas (1 to 4) and DAPI staining intensity was quantified using the cell profiler software. The mean fractional intensity at a given radius was calculated as the fraction of the total intensity normalized to the fraction of pixels at a given radius in n=584 control and n=762 HP1-TKO (TKO) nuclei.
Data are the mean ± SEM. ***p value <0. 001. (E) ChIP analysis of H3K9me3 and H3K27me3 enrichment on major satellite repeats within control (Ctl, C3 and C5) and HP1-TKO (TKO, KO1 and KO3) BMEL cells. Graphs represent the average of 5 independent experiments and are shown as the mean ± SEM. ***p value <0. 001, ns: not significant (Student's t-test). (F) Loss of the three HP1 proteins did not affect the expression of major satellites in neither liver nor BMEL cells. qPCR assays were performed using total RNA from livers of 7-week-old control (n=4) and HP1-TKO mice (n=4) and on control (Ctl, C3 and C5) and HP1-TKO (TKO, KO1 and KO3) BMEL cells. data were normalized to Hprt expression and are shown as the mean ± SEM. ns: not significant (Student's t-test). (G) Satellite repeats were quantified by qPCR on genomic DNA from the same animals as those used for (F).
ACKNOWLEDGMENTS: We thank  P. Chambon ,  C. Sardet ,  T. Forné ,  D. Fisher  and  C. Grimaud  for helpful discussions and critical reading of the manuscript. We thank  F. Bernex  and  L. LeCam ,  C. Keime  and B. Jost for fruitful discussions. We thank  L. Papon ,  H. Fontaine  and  C. Bonhomme  for technical assistance and  M. Oulad-Abdelghani  and the  IGBMC  for the anti-HP1 and TRIM28 antibodies. We also thank the RHEM technical facility and particularly  J. Simony  for histological analysis and the IGBMC/ICS transgenic and animal facility for the initial establishment of the HP1 and TRIM28 mouse models. We thank  C. Vincent  and the  IRCM animal core facility  for the day to day care of the animal models. Finally, we thank S. Chamroeun for counting positive cells on TMA. We acknowledge the imaging facility MRI, member of the national infrastructure  France-BioImaging  and supported by the  French National Research Agency  ( ANR-10-INBS-04 ,  «Investments  for the future»).
This work was supported by funds from the  Centre National de la Recherche Scientifique (CNRS) , the  Institut National de la Santé et de la Recherche Médicale (INSERM) , the  University of Montpellier  and the  Institut regional de Cancérologie de Montpellier (ICM) , ) and  SIRIC Montpellier Cancer , Grant INCa_Inserm_DGOS_12553. SH was funded by an  Erasmus PhD  fellowship. We also thank the  Ministry of Education, Science and Technology of the Republic of Kosovo  for a scholarships to support SH. FC was supported by grants from  ANR  ( ANR 2009 BLAN 021 91 ;  ANR-16CE15-0018-03 ),  INCa  ( PLBIO13-146 ),  ARC  ( PJA20131200357 ), and  La ligue Régionale contre le Cancer  ( 128-R13021FF-RAB13006FFA ). Sequencing was performed by the MGX facility. Montpellier, France.
To alleviate the requirement of a reference, recent methodologies use k-mer representation; they directly compare the counts of nucleotide sequences of length k between samples  [6]. These k-mer based approaches have been core to the field of metagenomics, where they are used to discover unique k-mers or k-mer signatures to classify organisms  [7, 8]. However, when translated to mammalian genomes, k-mer representation results in a k-mer count matrix with as many columns as there are samples and as many rows as there are k-mers, generally billions. Exploring such large matrices to find biologically relevant k-mers is intractable unless the analysis focuses only on a very small subset of the sequencing data  [5]  or by using metaheuristics that provide partial solutions  [9]. Here we present iMOKA (interactive multi-objective k-mer analysis), a novel approach and software that allows non-specialists to make use of k-mers to explore large amounts of mammalian sequencing data.
This approach is agnostic of the type of sequencing data used, is not biased towards annotated genetic elements, and can analyze transcript levels and single nucleotide variations in one pass. Importantly, iMOKA is interactive; it allows the user to import and merge samples from different studies and tailor their exploration of k-mers to specific genomic elements of interest such as splicing events, mutations, or global gene expression. We tested iMOKA on four clinical datasets: the classification of breast cancer subtypes and response to chemotherapy of breast, ovarian cancer, and diffuse large B cell lymphoma (DLBCL). We find that iMOKA found features that are more accurate than classical bioinformatics approaches, takes up less space, uses less memory, has faster runtimes, and can be run on a computer cluster or on a laptop.
Results iMOKA design iMOKA imports sequencing files in FASTQ, FASTA, BAM format, or SRR identifiers via its user interface. It then counts the occurrences of all sequences of given length k (default 31)  [9]  using the KMC3 software  [10]  in each sample (Fig. 1 ). It then extracts labels from the sequencing metadata so that the user can define groups they wish to compare. Importantly, each sample is stored as a sorted vector of k-mer counts in a dedicated binary file using a custom prefix-suffix structure that drastically reduces the disk space requirements ("Methods" section). For each sample, a JSON file is created that contains metadata and a rescaling factor for k-mer count normalization that allows the user to remove or add samples without having to recalculate an entire k-mer matrix.
Benchmarking datasets and algorithms iMOKA uses a k-mer based analysis to detect sequence features and create classification models from large cohorts of mammalian RNA sequencing data. To test its performance, we selected four studies that were distinct in their data structures, classification objectives, and sizes. The first was a non-binary classification of 1038 patients aiming to define 4 subtypes of breast cancer which were luminal A (LumA), luminal B (LumB), HER2-enriched (HER2), and basal-like. The second was a cohort of 240 ovarian cancer patients where the objective was to predict response to chemotherapy. The third was a smaller cohort of 118 breast cancer patients where the objective was also to predict response to chemotherapy. The last was an even smaller cohort of 17 DLBCL patients divided according to their responsiveness to the chemotherapy.
In our benchmark, we included methods based on four different types of features which were k-mer counts, percentage-spliced-in (PSI), transcripts per kilobase million (TPM), and sequencing counts. The two latter were measured and tested across annotated genes and transcripts separately. The algorithms we benchmarked were DESeq2  [12] , edgeR  [13] , and limmaVoom  [14]  for TPM and sequencing counts; iMOKA for kmer counts; and Whippet  [15]  for alternative splice site usage. We excluded four other k-mer based methods HAWK  [16] , KOVER  [17] , Kissplice  [11] , and GECKO  [9]  because they were respectively impossible to run on such big datasets due to segmentation fault errors, were unable to find k-mers that could classify the input samples or, for the last two methods, were killed after 2 weeks of runtime on our computer cluster.
In our benchmark, we compared the list of features output by each algorithm by using them in a random forest classifier and determining their out of bag scores (OOB score). The out of bag score tests how well each classifier performs without having to set aside a portion of the data specifically as a test set. It is as reliable as using a test set  [18, 19]  without having to set aside part of the data. We chose the random forest classifier because it is a non-parametric approach and because the importance of each input feature is easy to evaluate. Finally, for the largest dataset, the molecular classification of breast cancer, we performed a 5-fold cross validation of the entire iMOKA procedure and all other benchmarked algorithms, using 4/5 of the dataset for data reduction and creation of a random forest model and 1/5 of the dataset as the test set.
Classification of breast cancer subtypes Breast cancer is a transcriptionally heterogeneous disease with multiple subtypes that determine prognosis, treatment, and patient outcome. Although breast cancer classification is constantly being updated, a broadly accepted stratification defines four groups which are luminal A (LumA), luminal B (LumB), HER2-enriched (HER2), and basallike  [20]. We benchmarked iMOKA on a dataset of 1038 mRNA-Seq breast cancer samples from the Cancer Genome Atlas (TCGA) Pan-Gyn cohort  [21]  (patients per class: basal 190, Her2 82, LumA 559, LumB 207) and tested how well the outputs of each approach could accurately predict the four classes. We found that the list of kmers output by iMOKA (Additional file 1, Fig. S5 ) was above all other methods in their ability to classify the four types of breast cancer (Fig.
2d  and Additional file 4). iMOKA identifies events associated with the response to treatment in ovarian cancer patients Our second benchmark was performed on a dataset of high-grade serous ovarian cancers taken from the TCGA_OV cohort  [29]. We included patients having an annotated  [30]  response to a first-line treatment to the combination platinum and taxane chemotherapy (patients per class: 174 responsive, 66 non-responsive). iMOKA identified 138 k-mers with individual accuracy between 65 and 75% (Table  S1  and Additional file 2). Again, the k-mers found by iMOKA gave the most accurate oob scores for response to chemotherapy (Fig. 3a ). The gain compared to other methods is much higher than for the previous breast cancer classification. This can be explained by the fact that most of the methods we benchmark against only make use of gene or transcript expression or splicing sites.
In agreement with this, the gene ontology (Fig. 3d ) analysis shows a functional enrichment of small molecule binding proteins. iMOKA identifies events associated with the response to neoadjuvant chemotherapy in breast cancer patients The third test dataset was taken from the Breast Cancer Genome Guided Therapy (BEAUTY) study  [31]  and consisted of patients with all 4 types of breast cancer for which we tested the response to neoadjuvant chemotherapy with paclitaxel and anthracycline. This allowed us to test the binary classification of more heterogeneous cell populations on smaller sample sizes: 36 patients that had a complete response to chemotherapy and 82 that did not. It is worth noting that this dataset presented a significant batch effect, detected using the R package DASC  [32] , associated with the load date of the samples (Fig.