dataesr/scientific-paragraphs-categorization huggingface dataset:
ACKNOWLEDGEMENTS We thank the anonymous referee for their suggestions that greatly improved the quality of the manuscript. We also thank  Alice Quillen ,  Friedrich Anders , and  Kyle Oman  for their comments and enlightening discussions. Funding for this work and for RAVE has been provided by: the  Australian Astronomical Observatory ; the  Leibniz-Institut für Astrophysik Potsdam (AIP) ; the  Australian National University ; the  Australian Research Council ; the  European Research Council  under the European Union's  Seventh Framework Programme  (Grant Agreement 240271 and  321067 ); the  French National Research Agency ; the  German Research Foundation  ( SPP 1177  and  SFB 881 ); the  Istituto Nazionale di Astrofisica at Padova ; The  Johns Hopkins University ; the  National Science Foundation of the USA  ( AST-0908326 ); the W.
Conclusion: Using high throughput strand-oriented sequencing for whole ticks in different stages and feeding conditions, we obtained a de novo assembly that significantly increased the genomic resources available for I. ricinus. Our study illustrates the importance of analysing the transcriptome at the level of the whole body to gain additional insights into how gene expression changes over the life-cycle of an organism. Our comparison of several RNAseq datasets shows the power of transcriptomic data to accurately characterize genetic polymorphism and for comparing different populations or sources of sequencing material.
ricinus, which was expected to show low heterozygosity. We then compared the results from our study to the results from previously published RNA-Seq studies that used different sources of I. ricinus tick material: (i) wild ticks; (ii) F1 offspring obtained from a mating between two wild ticks; and (iii) a tick cell line. Specifically, we compared levels of polymorphism and heterozygosity between four different sources of I. ricinus tick material that were expected to differ with respect to their genetic variability. This information provides a large catalogue of polymorphic sites in expressed regions and provides an important database for future population genetic studies.
Quality and assembly of reads To produce the dataset, the raw paired-end reads (2 × 125 bp) were first cleaned. Adapters were clipped and low-quality regions were filtered using Trimmomatic (release 0. 36)  [21] ; only reads with a minimum of 36 high quality-scored contiguous bases were kept. Summary statistics of the sequence quality were checked for each library by visualizing the FastQC report (release 0. 11. 5) [22].
Assessment of transcriptome completeness A common test used to assess the "coverage" (or information completeness) of a given sequence dataset is to analyse random samples of reads, and to create a saturation curve describing the relationship between different metrics (numbers of contigs of a specified length, number of matches to a known set of genes, etc. ) and the read sample size. Complete datasets have saturation curves that plateau more quickly compared to incomplete datasets. To determine the coverage or completeness of our dataset, we randomly sampled 1, 2, 5, 10, 20, 50, 80, 100, 140 and 160 million reads from our cleaned libraries, as detailed below. Completeness was also estimated by determining the presence of homologs of conserved arthropod genes using the BUSCO approach  [25].
BUSCO v1 uses a reference database of 2675 conserved arthropod genes (BUSCO genes) and searches for potential homologs in the database of interest by running BLAST  [26]  and HMMER  [27]. Conserved BUSCO genes are assigned to four classes of genes: (i) missing; (ii) fragmented; (iii) duplicated; and (iv) complete. To determine if the open reading frames (ORFs) were correctly predicted, we checked the strands of the predicted genes within the contigs matching the BUSCOs. To produce the final assembly, we reduced the potential redundancy resulting from the presence of alternative transcripts in the contigs. We clustered similar sequences using cd-hit-est  [28]  with 98% of identity, retaining the longest transcript of each cluster. Identity parameters were chosen to cluster nearly identical sequences resulting from alternative splicing.
Gene prediction and annotation The prediction of coding sequences (> 100 amino acids) was performed using TransDecoder, which is part of the Trinity software  [24]. The TransDecoder options were set to account for the strand orientation of the sequencing (i. e. ORFs were searched only on the forward strand). Finally, annotations from comparison with public databases were used to filter multiple ORF predictions by transcripts (see below). Following the pipeline recommendation of Trinotate (release 2. 0. 2)  [29] , both contigs and predicted peptides were compared by blastx+ and blastp+ (release 2. 2. 29)  [26]  to releases of Swissprot and Uniref90 (available at https://data. broadinstitute. org/Trinity/__depreca-ted_trinotate_resources/Trinotate_v2. 0_RESOURCES/ v2. 0 RESOURCES/). Protein domains were identified using HMMER (release 3.
0 from March 2010)  [27]  with PFAM-A  [30] , signal peptides with SignalP  [31] , and transmembrane domains with TMHMM  [32]. We tagged ribosomal RNAs using RNAmmer  [33]. All these layers of annotation were combined by Trinotate to assign gene ontology (GO) information to each contig. In the case of multiple ORF predictions for a contig, if one ORF was similar to a known protein while the others was not, only the former ORF was retained. In other cases (several ORFs with similarity to known proteins, or several ORFs with no similarity to known proteins), the different ORFs were retained. As we were particularly interested in cuticular proteins, all peptides that contained the chitin-binding domain (PF00379) were classified using the CutProtFam-Pred webserver  [34].
Comparison of the completeness of our assembly relative to other assembled transcriptomes In recent years, several research groups have produced RNAseq datasets for I. ricinus. For most of these projects, the de novo assemblies (or sets of predicted genes derived from these assemblies) have been published in the Transcriptomes Shotgun Assembly (TSA) division of GenBank. Using statistics provided by BUSCO  [25] , we compared the completeness of our final assembly to that of six different assemblies/gene sets, which were obtained from six different RNAseq projects (see Additional file 1: Table  S2 ; TSA accessions: GADI01, GANP01, GBIH01, GCJO01, GEFM01 and GEGO01). In addition, to determine the relative contribution of each of the different datasets (TSA and or own CDS prediction) to the complete gene collection, we analysed the clustering of all CDSs using cd-hit-est with default parameters  [28]  -for GCJO01, which corresponded to contigs, we predicted CDSs by using Transdecoder  [28].
Differential expression and GO enrichment In addition to the reconstruction of transcript sequences, RNA sequencing also allows the user to quantify transcript expression by counting the number of sequenced reads that map to a given transcript. Paired reads for each library were pseudo-aligned on the Transcriptome de Bruijn Graph (T-DBG) using Kallisto  [35]. We chose this method, based on a k-mer approach, because it is much faster while providing the same accuracy as the best mapping approaches  [35]. This method produced raw counts and normalized count statistics (TPM, or transcripts per million reads) for each assembled contig. These counts allowed us to test for differential expression between conditions. To assess feeding-related changes in gene expression, we compared partially fed ticks versus unfed ticks.
For this, we performed a comparison between two ensembles of libraries, respectively D/E/F/M/N/O and A/B/C/G/H/J/K/ L. Differential expression analyses were performed with the R package DESeq2  [36]  using the raw counts from Kallisto, each library been taken as an independent replicate. To describe the relevant biological changes between conditions, we used predictions produced by gene ontology (GO) term annotation, which included: "molecular function", "biological process" and "cellular localization"  [37]. GO terms were compared between transcripts that were not differentially expressed between conditions (unbiased transcripts) versus transcripts that were differentially expressed between conditions. We defined unbiased transcripts as contigs with no significant change in expression between conditions (fold change less than 2 and adjusted P-value higher than 0.
05). Enrichment analysis was performed using the elim method using Kolmogorov-Smirnov tests developed and implemented by Alexa  [38]  in the R package TopGO. As suggested by this author, multiple testing was taken into consideration by using the false discovery rate (FDR) on the enrichment test P-values. The resulting GO enrichments were analyzed using the R package GOprofiles  [39] , which provides visualization tools. GO enrichment comparisons are by definition limited to contigs with assigned GOs, whereas many more contigs can show significant changes in expression between conditions. A substantial number of contigs had no assigned GOs but did have other annotations, such as domains identified through the PFAM analysis. We therefore used a text mining analysis of the PFAM domains to further compare changes in gene expression between the different conditions (unfed/fed, nymph/adult, male/female).
This approach allowed us to distinguish the most common terms associated within a text. We therefore extracted the PFAM terms associated with contigs over-expressed in each of the conditions defined above (fed or unfed) and treated them as a single text for each category. To prevent over-representation of transcripts with many PFAM domains, only the first ten PFAM terms with the lowest e-value were retained for a transcript and over-expressed contigs were defined as contigs with a fold change larger than 4 and a significant P-value (P < 0. 05) in the DESeq2 analysis. PFAM descriptions were edited to remove the less informative terms (e. g. "protein", "domain", "motif"). We then used an in-house R script to draw clouds of words with word sizes proportional to their frequency in the text. This approach is complementary to the GO-enrichment tests, as it helps to visualize major shifts in expression among conditions.
Polymorphisms Polymorphism was surveyed in the reads produced through our project (corresponding to an inbred line, here after referred to as NEU) but also for data from five other RNAseq projects of I. ricinus publically available in GenBank. Three published datasets used wild tick populations from Sénart in France (SEN)  [16] , from the Czech Republic (CZ-W)  [13] , and a mixture of wild tick populations (LUX) that was provided by Charles River Laboratory  [17]. The other two datasets were based on F1 full sibs from a cross between wild ticks from the Czech Republic (CZ-F1)  [19] , and a tick cell line (CL) deposited by the Broad Institute under BioProject accession numbers PRJNA238785-88. More details on those 5 datasets are given in Additional file 1: Table  S3. Reads sequences are available at the NCBI Sequence Reads Archive (SRA) and organized by BioProject.
After downloading reads from the SRA archives, the reads were cleaned using Trimmomatic (with the same parameters as above). As estimates of polymorphism and heterozygosity depend on sample size, we standardized the sample size for each dataset by randomly sampling 30 million reads from each SRA. The combined datasets were analysed to detect single nucleotide polymorphisms (SNPs). SNPs were predicted with a "direct from the reads" approach, using KisSplice (release 2. 4. 0) [40] , a software that identifies variations by detecting "bubbles" in the De Bruijn graph. SNPs were mapped on our final set of contigs using Blat (version 36)  [41]. A report assessing various parameters for each SNP (location, reliability, etc. ) was provided by Kiss2refTranscriptome  [42]. To minimize false positives, we retained only SNPs that respected the following criteria: (i) SNPs had to be covered by at least ten reads in each dataset; and (ii) SNPs needed to be uniquely mapped (e.
Taxonomic assignation We used taxonomic information from the code name of the best hit on the Uniref90 proteins cluster. We found that 6. 4% of assembled transcripts (n = 12,368) were assigned to fungi (considering a minimum of 50% of protein identity and an E-value lower than 10e-5). Fungal contamination of individual ticks has been observed in the I. ricinus colony at the University of Neuchâtel. To facilitate moulting, blood-engorged larvae are placed in tubes with moistened filter paper, which also facilitates the growth of opportunistic fungi, which are the most likely source of the contamination observed in the present study. We counted 887,767 reads that mapped to the fungi-like contigs, which represents only 1. 2% of the total number of counting events by Kallisto, suggesting that contamination was modest.
For both species, the highest peak of genes had a very high percentage of identity (> 95%). This peak probably corresponds to the orthologous genes found in both I. scapularis and I. ricinus. However, the distribution of identities for I. ricinus shows a secondary peak of hits with a low identity (∼45%). These low identity hits cannot correspond to the same gene, but probably represent distant paralogs of genes, or genes having a similar protein domain. The presence of this secondary peak decreases the mean identity of hits to I. ricinus and explains why the mean identity of our dataset is lower for I. ricinus than I. scapularis.
Completeness Of the 2675 conserved arthropod genes in the BUSCO database, our assembly contained 2033 complete genes (completeness of 76%) and 250 partial genes (extended completeness of 85. 3%). Previously published assemblies or collections of predicted genes (TSA archives) for I. ricinus all produced lower percentages of complete conserved genes (Fig. 2a ). A comparison of the overlap between our new dataset and the combined TSA datasets found that our assembly contains 295 BUSCO genes (11%) not present in any of these TSA datasets, whereas the combined TSA datasets contained 231 BUSCO genes not present in our assembly (Fig. 2b ). We obtained Fig. 1  Taxonomic assignment of contigs using Uniref90 best blast hits. The bars correspond to the number of contigs assigned to different taxonomic groups (the percentage of contigs assigned to each taxonomic group is shown in brackets).
Best hits to "Acari", which include ticks and mites are hits to Acari species other than I. ricinus and I. scapularis. Best hits to "Arachnida", which also include spiders and scorpions, are hits to Arachnida species other than Acari. The mean percentage of identity at the amino acid level and the standard deviation are shown after the taxon name a similar percentage of BUSCO contigs with correctly predicted strands (98. 7%) compared to the combined TSA datasets (mean of 98. 6%, minimum of 95. 8% for GEGO, maximum of 100% for GANP). We then performed a direct comparison between all the CDSs from our assembly and from the TSA datasets, with a clustering approach. We showed that our assembly permitted to identify 36,652 new CDSs which were not present in any of the published dataset (Fig. 2d ).
3 ); specifically, the expression of cuticle-associated genes was significantly higher in partially fed ticks than unfed ticks. Out of 39 transcripts containing a chitin-binding PFAM domain (PF00379), one contained two PF00379 domains and 22 contained a signal peptide. All these transcripts were classified as members of group 2 using the CutProFam-Pred webserver. Expression levels in each library of cuticle-related transcripts are shown in a heatmap (Fig. 4 ). Of the three most expressed cuticular transcripts, two transcripts shared a high identity with the proteins described by Andersen & Roepstorff  [43]. The c252234_g1_i1 transcript in our study had 97. 03% identity with the Ir-ACP10. 9 protein (belonging to the CU109 cluster of Uniprot)  [43]. This transcript was found to be highly over-expressed in partially fed ticks (588-fold change between unfed and partially fed ticks; P-value < 0.
Polymorphism Starting from the 368,367 SNPs initially detected by KisSplice using the default parameters, we applied more stringent criteria to increase the specificity: only 8955 SNPs were supported by at least 10 reads in each strain. Fig. 3  Abundance of significantly enriched GO terms in partially fed ticks. For comparison, the abundance of the same terms for differentially expressed transcripts in both conditions (unfed and partially fed ticks) are also represented. GO enrichment for the comparison between unfed and partially fed ticks. The X-axis shows the top enriched GO terms for the three categories: Biological Process, Molecular Function and Cellular Component. The Y-axis shows the number of GO terms found in differentially expressed transcripts for unfed and partially fed ticks (see Additional file 1: Tables S4-S9 for detailed GO results).
Discussion Several recent RNAseq studies on I. ricinus have focused on the transcriptomes of specific tissues such as the midgut, salivary glands, haemocytes and ovaries  [13] [14] [15] [16] [17] [18] [19]. In contrast, we tried to obtain a broader picture of the transcriptome of ricinus by sequencing and assembling transcripts from whole ticks in different conditions defined by developmental stage, sex and feeding status. First, our de novo assembly showed high overall completeness and significantly enlarged the catalogue of known coding sequences for I. ricinus. Secondly, we detected blood-feeding-induced changes in gene expression at the level of the whole body. Thirdly, our analysis of polymorphism in the transcriptome sequence data of our study and previously published studies allowed us to compare the levels of genetic diversity between outbred wild strains of I.
We found again that large numbers of CDSs were predicted only either in our study or in published datasets. Therefore, each analysis had a relatively large level of specificity, a finding thus similar to that of the BUSCO approach. This suggests that it is important to combine the different sources of data to obtain the most complete description of the transcriptome of the species. When analysing shared CDSs in the different analyses (by a clustering method), we found that our predicted CDSs often provided the longest sequence in each cluster, i. e. the most complete gene sequence. Therefore, our predicted collection of CDSs significantly enriched the extant knowledge for I. ricinus genes. To understand the differences in completeness between datasets (and their complementary aspects), we stress that tissue-specific transcriptomes (previous studies) will detect transcripts that may be rare at the level of the whole body (our study), whereas our whole-body transcriptome may have helped to identify transcripts present mostly or only in tissues not already covered by previous studies.
scapularis than to I. ricinus. One explanation for this counter-intuitive result is that the available genomic resources are more complete for I. scapularis than I. ricinus. A complete genome has been published for I. scapularis  [46] , but only a genome survey is available for I. ricinus  [17, 47]. For I. ricinus, most protein sequences in the databases are derived from de novo assembled transcriptomes that are still rather incomplete. Thus, one explanation for the relatively low mean identity of the hits to I. ricinus could be the absence of the same gene in the data banks (matches would often correspond to paralogous gene copies or shared protein domains). Gene expression in response to blood-feeding Differential expression (DE) analyses were focused on the comparison between two feeding conditions, the unfed and partially fed ticks.
Genetic diversity We therefore established a catalogue of robust SNPs which could be useful for future population genetics analyses, as it could be used to design a genotyping study of multiple individuals, allowing a shift in scale compared to traditional population genetics approaches based on small numbers of markers. To our knowledge, only one similar analysis with RNAseq data has been conducted to date for I. ricinus  [14] , but the numbers of identified SNPs between this study and our work are difficult to compare given the different strategies (de novo SNP identification from the reads versus mapping to a set of CDSs) and because the reference sets used for localization of the SNPs are not the same. We however stress that the de novo (direct from the reads) approach for RNAseq datasets has been successfully applied in recent studies for different organisms  [50].
Even though the different datasets analysed in this work comprise pools of individuals or materials of very different origins, our approach illustrates the power of these analyses to determine the genetic characteristics of the different materials. Our study of single nucleotide polymorphisms found indeed substantial differences in heterozygosity between the different sources of I. ricinus used in RNAseq studies. The factorial component analysis indicated that three datasets corresponding to wild ticks (SEN, CZ-W and LUX) clustered together, whereas each of the three laboratory strains (NEU, CZ-F1 and CL) differed markedly from that cluster and from each other. The clustering of the wild tick populations indicates that they have similar levels of heterozygosity (and should not be interpreted as the absence of genetic diversity between the three wild tick populations).
Conclusions Our de novo assembly showed high overall completeness and significantly enlarged the catalogue of known coding sequences for I. ricinus. Our study investigated the transcriptome of whole ticks that differed with respect to their developmental stage, sex and feeding condition. This approach allowed us to detect changes in gene expression at the level of the whole body instead of specific tissues. We found that blood-feeding induced a strong upregulation of transcripts associated with cuticle production. Finally, our analysis of polymorphism in the transcriptome sequence data of our study and previously published studies allowed us to identify 3866 robust SNPs from expressed-regions, and to compare the levels of genetic diversity between outbred wild strains of I. ricinus versus inbred laboratory strains and tick cell lines.
Acknowledgements We are grateful to the  Genotoul bioinformatics platform Toulouse Midi-Pyrénées (Bioinfo Genotoul)  for providing support and computing resources. Part of the computations (multiple transcriptome assembly for subsamples of reads) were performed at the  Bordeaux Bioinformatics Center (CbiB). We thank two anonymous reviewers and the editors for their constructive criticism and comments which helped improve this manuscript.
Availability of data and materials Reads from sequencing the 15 libraries were deposited in the SRA section of NCBI under the BioProject ID: PRJNA395009. The contigs are available under accession GFVZ00000000 (NCBI, TSA). The annotation table was deposited on the Zenodo repository with a DOI (10. 5281/zenodo. 1137702), available on https://zenodo. org/record/1137702.
By applying a less stringent gene definition, which includes 50 kb upstream of the TSS, we identified a pool of 16,406 targeted genes enriched in 35 gene ontology functional categories (Table  S2 ) and 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Table  S3 ) by Database for Annotation, Visualization and Integrated Discovery (DAVID) analysis (p < 0. 05 after Bonferroni correction for false discovery rate). Enriched pathways included metabolism and transport of nucleic acids and proteins, chromatin modifications, regulation of gene expression and mitosis, ubiquitination, splicing, and endocytosis.