dataesr/scientific-paragraphs-categorization huggingface dataset:
The 45 resequencing datasets were mapped on the complete Chinese Spring (CS) genome sequence (IWGSC RefSeq_v1. 0). On average, 98. 6% of reads were mapped, demonstrating the absence of a major contamination in our samples (Table  1 ) and that sequence divergence between different Triticum genomes did not prevent from mapping reads accurately on the CS reference. On average, 81% of mapped reads were mapped onto the chr3B pseudomolecule. The remaining 19% are mapped onto unanchored scaffolds (chrUn; 4% of mapped reads) but also to other chromosomes, suggesting the presence of DNA from other chromosomes in the sorted DNA samples (especially chr2B whose size is close to that of chr3B; Additional file 1: Supplementary Table  S1 ). The distribution of mapped reads along the 20 other chromosomes did not reveal any case of homeologous exchange with 3A or 3D, neither translocated from other chromosomes.
They represent on average 145 genes (2%) per accession, ranging from 42 to 217 (0. 7-3. 5%; Table  2 ; Additional file 1: Supplementary Table  S3 ). As expected, East Asian T. aestivum accessions exhibit the lower rate of variations with 75 genes (1%) showing PAVs, while T. carthlicum exhibit the highest: 184 (3%) genes. Only one gene is absent in all tested accessions, meaning it is specific to CS: a cupredoxin-encoding gene (TraesCS3B01G482300). A de novo assembly of the unmapped reads was conducted for each sample to identify genes that are absent from the reference. We retrieved between 695 and 11,173 contigs (>500 bps) per dataset representing between 620 kb and 7. 96 Mb (<1% of chr3B size). We searched for similarity with known proteins from related plants (see Materials and Methods) and applied a clustering step to detect genes present in several accessions.
aestivum (574 TEs) represents less than 1% of the TE content of the chromosome. Thus, it is obvious that TE transposition is not the main factor of this variability. This was confirmed by the strong conservation of the TE copy number estimated for each family. Indeed, only a few families are present in different proportions between accessions. This supports previous findings suggesting that the TE composition is highly controlled since A-B-D divergence, with no major TE burst  (Wicker et al. , 2018). Most of SVs affecting TEs originate from deletions, agreeing with previous conclusions made based on MITEs  (Keidar-Friedman et al. , 2018). We used SV data to study the phylogenetic relationships between Triticeae. The tree topology suggests that hexaploid species originated from independent hybridization events and do not share a common tetraploid parent.
Phylogenetic Analyses To get robust markers, we first pruned SNPs using plink (parameters: -blocks no-pheno-req -allow-no-sex -blocks-maxkb 500000 -indep-pairwise 500000 5 0. 99; v1. 90b6. 2; Purcell et al. , 2007)  while keeping undetermined values (NA) in our genotyping matrix. We then used R package pvclust on both SNP and PAV matrix to perform a hierarchical clustering (function: pvclust(), distance: simple matching, clustering: Ward. D2, Bootstrap: 10000). Dendrograms were then compared and aligned using R package dendextend  (Galili, 2015). Entanglement score estimated the quality of the dendrogram alignment, zero being perfectly untangle (SNP tree vs PAV tree = 0. 01). To evaluate the similarity between the two dendrograms, Baker's gamma index  (Baker, 1974)  was calculated using dendextend function cor_bakers_gamma.
Forty-seven genes not covered after mapping CS resequencing dataset were filtered out, leaving 6,082 genes that could be analyzed. Genes were considered absent (PAV) when covered over less than 10% of their length, while considering only coding exons. For the detection of CNVs affecting genes, we calculated a normalized ratio between the read count for each gene in accession X and the read count in Chinese Spring. First, we normalized the read counts based on the proportion of the length of the gene covered by reads. Second, for each accession, we normalized read counts based on the total number of reads mapped on all chromosome 3B genes. Finally, we computed a log2ratio of the normalized read counts for each gene between accession X and Chinese Spring. For each accession, a boxplot of the 6,082 gene log2ratios was computed and we used the R native function boxplot upper and lower whiskers as thresholds to identify potential genes with up-or downCNVs, i.
Detection of TE inDels and CNVs We developed a Perl script (getISBPfromBed. pl; available at 1 ) in order to extract 150 bps encompassing the junction 1 https://github. com/fchoulet/biotools between every TE predicted across CS chromosome 3B and its insertion site, each corresponding to a potential 5 and 3 ISBP marker. In case of major overlap between two neighbor ISBPs (>99 pbs), we kept only one ISBP to avoid redundancies. ISBPs containing one or more Ns were filtered out. In addition, we discarded non-unique ISBPs (detected by mapping ISBP sequences using BWAmem). These later represent typical cases where a TE was predicted as two or more adjacent fragments, leading to design ISBPs that are fully included into an element and, thus, are repeated. TE inDels were assessed by mapping reads on the set of ISBPs while keeping only 100-bp reads fully mapped with at maximum 2 mismatches, meaning that the junction between the TE extremity of the insertion site is covered.
Detection of PAVs and CNVs Affecting the Non-coding Non-repeated Fraction of Chromosome 3B We built a BED file as a template to call PAVs only in genes, non-coding low-copy DNA, and TE extremities (ISBPs). We retrieved the positions of low-copy regions (including genes) as the complement to the TE positions. These segments were split into 100-bp windows, while discarding all segments smaller than 51 bps. We then added the positions of the complete list of 150 bp ISBPs to the BED file. Altogether, these regions account for 134 Mb split into 1,593,666 loci of 100/150 bps. Presence/absence of each locus was computed by analyzing read mapping data as described above. This BED file was also used to detect deletion blocks. Blocks were defined using a sliding window of 50 neighbor loci (with step = 5).
Assembly of Unmapped Reads and Gene Prediction We used samtools to select unmapped reads (-F4) and to create 2 paired-end and 1 single end fastq files (-f 64, -f 128). We assembled these reads independently for each accession using MIRA (parameters: genome, de novo, accurate; v4. 0. 2, Chevreux, 2004)  and kept only contigs larger than 500 bps. Potential contaminations were searched by aligning all contigs against the NCBI nt database (August 2018 2 ) using BLASTn (parameters: -F F -e 1e-5,  Altschul, 1997). Contigs with best hits outside the plant kingdom (id 80%, length overlap cutoff 200 bps) were considered contamination. Then, we searched for the presence of coding genes using Exonerate (options -maxintron 20000 -showtargetgff TRUE -model protein2genome; v2. 2. 0  Slater and Birney, 2005)  and the predicted proteins originating from 9 plant genomes: Sorghum bicolor, Zea mays, Oryza sativa, Brachypodium distachyon, Aegilops tauschii, Hordeum vulgare, Triticum urartu, Populus trichocarpa, and Arabidopsis thaliana.
ACKNOWLEDGMENTS We thank  Dr. Marie Kubaláková  and  Bc. Romana Šperková  for technical assistance with chromosome flow sorting and chromosome DNA amplification. We thank the  INRAE URGI platform (Versailles, France)  and the  Mésocentre UCA  for their support and computing facilities. We also thank the  Biological Resource Center on small grain cereals  ( INRAE, Clermont-Ferrand ) for providing seed samples.
Such aberrant transcription and translation events lead to the generation of peptide sequences that are missing in conventional protein sequence repositories  8, 9. If such translation events lead to the presentation of tumor-specific and immunogenic HLA ligands, these occurrences could substantially expand the repertoire of targetable epitopes for cancer immunotherapy  [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]. Currently, ~1% of the entire genome is annotated as protein-coding regions, yet 75% of the genome can be transcribed and theoretically translated, potentially offering a pool of previously unexplored peptide targets  20. To date, mass spectrometry (MS) is the only analytical methodology that allows the direct identification of the HLAp repertoire in vivo  21. Often, MS-based immunopeptidomic discoveries are limited to the standard, available protein sequence database, usually containing only annotated proteome-derived sequences.
All three forward open reading frames (ORFs) (stop-to-stop) in the above transcripts were subsequently in silico translated into polypeptide sequences. For every sample, the polypeptide sequences were concatenated to personalized canonical proteome references containing allelic variant information from patient tumor exome data. Finally, we searched the MS immunopeptidomics data against these personalized reference databases. Database size affects false positives in noncHLAp detection. In silico translation of transcripts in three forward reading frames results in a large number of potential polypeptide sequences. In proteogenomics, searching MS data against such inflated protein reference databases may propagate false positives  28, 29. Hence, our first investigative step was to understand the impact of database size on the level of false positives in immunopeptidomics datasets.
We searched reference databases containing canonical (i. e. , UniProt) and our non-canonical polypeptide sequences with a single search tool (MaxQuant) and at a global 1% FDR. The accuracy was assessed by assigning HLA binding prediction scores to the MS-identified peptides with MixMHCpred  30. We reasoned that non-canonical HLA class-I bound peptides (non-cHLAIp) should follow the same binding rules as protHLAIp  31. First, we compared a generic non-canonical protein sequence database derived from the three forward frame (three-frame) translation of all non-coding transcripts from ENCODE  32  with a sample-specific protein sequence database derived from the threeframe translation of lncRNAs and pseudogenes from the RNA-Seq data using an expression cutoff value of > 0 fragments per kilobase of transcript per million mapped reads (FPKM).
Additional databases of decreasing size were assembled by retaining only the sequences that originated from more highly expressed genes (FPKM > 2, > 5 or > 10). Reducing the size of the database by personalizing and focusing on highly expressed genes led to an increase in the percentage of noncHLAIp that were predicted to bind to their respective HLA alleles (MixMHCpred p-value ≤ 0. 05) (Fig. 1b ). Restricting the database to polypeptide sequences originating from highly expressed genes should, on the one hand, improve the accuracy of MS-based non-canonical peptide identification, while on the other hand, lead to the potential loss of peptides encoded by lower-expressing transcripts. Hence, in this study, we included all non-coding transcripts with FPKM > 0 to circumvent the need to exclude polypeptide sequences based on low-expressing transcripts.
NewAnce improves the accurate identification of noncHLAp. We developed the computational module NewAnce, which combines the MS search tools MaxQuant  33  and Comet  34  , with the implementation of a group-specific strategy for the FDR calculation (see Methods section for more information and Supplementary Fig. 1a -g for performance evaluation). All HLAp identified by either of the search tools were consequently matched against an up-to-date UniProt/TrEMBL sequence database (95,106 protein sequences of the human reference proteome (up000005640), with isoforms) to extract noncHLAp that do not map back to known human proteins in UniProt. For every sample, FDRs were calculated separately for protHLAp and noncHLAp identified by Comet (Fig. 1c  and Supplementary Fig. 1a ). Only consensus (intersection) peptide-spectrum matches (PSMs) from Comet and MaxQuant were retained for further downstream analyses.
a A schematic of the entire workflow is shown, where tissue samples or tumor cell lines were obtained from patients, and exome, RNA-and Ribo-Seq were performed to provide a framework to assess the non-canonical antigen repertoire. HLAp were immunoaffinity-purified from cancer cell lines and matched tumor/healthy lung tissues and then analyzed by MS. Immunopeptidomics spectra were then searched against RNA-and Ribo-Seq-based personalized protein sequence databases that contain non-canonical polypeptide sequences. MS-identified noncHLAIp were validated by targeted MS-based PRM and tested for immunogenicity using autologous T cells or PBMCs. b The percentage of predicted HLA binders of length 8-14 mer peptides with a MixMHCpred p-value ≤ 0. 05 was used to evaluate the accuracy of the identified HLAIp by MaxQuant at 1% FDR as a function of database size (blue line).
The percentage of predicted binders obtained for each condition is shown for each bar for the melanoma cell line 0D5P. c Different protein sequence databases combining whole-exome sequencing and inferences from RNA-Seq and Ribo-Seq data were utilized. NewAnce was implemented by retaining the PSM intersection of the two MS search tools MaxQuant and Comet, and applying group-specific FDR calculations for protHLAp and noncHLAp. Source data are provided as a Source Data file. Low RNA expression limits noncHLAIp presentation. We then characterized the expression levels of source RNAs encoding HLAIp in more depth. For this purpose, we compared all identified source genes of protHLAIp to source genes of noncHLAIp in the 0D5P sample. The protein-coding source genes had a median FPKM value of 9. 3, whereas the presumed non-coding source genes showed lower expression overall, with a median FPKM of 2.
Nonetheless, we identified and confirmed the presence of a unique DAC-induced noncHLAIp derived from a lncRNA (Supplementary Fig. 4r ). Ribo-Seq improves the coverage of protHLAIp and non-cHLAIp. Next, we hypothesized that immunopeptidomes would correlate more closely with translatome than transcriptome data. To build the translatome-based database for the MS search, all ORFs showing periodic RPF distribution were extracted for the 0D5P sample, and translated in silico. This technique reduced the size of the search space, and we used this independent discovery method in our study to identify additional noncHLAIp, including those derived from previously unexplored ORFs in coding genes. We investigated the extent by which a protein sequence database inferred by Ribo-Seq could replace or complement the search performed with our personalized references comprising canonical protein sequences concatenated with polypeptidecsequences from the three-frame translation of expressed non-coding transcripts.
Using 0D5P as a representative immunopeptidomic dataset, we  observed a positive correlation between RNA expression and HLAIp-sampling (see Methods section) searched against the personalized protein sequence database (r = 0. 392) (Fig. 4e ). Then, we searched the same immunopeptidomics MS data against the de novo-assembled Ribo-Seq inferred database, and we correlated this HLAIp-sampling with RNA abundance (Fig. 4f ) or with translation rates based on the spectral coefficient of the 3periodic signal in the Ribo-Seq data (see Methods section) (Fig. 4g ). This approach resulted in a significantly higher positive correlation between HLAIp-sampling searched against the Ribo-Seq inferred database and the translation rate (r = 0. 574) than the overall RNA abundance (r = 0. 431, two-sided p-value < 10e-16).
Thus, evidence exists that the immunopeptidome, at least for the 0D5P sample, is better captured by the translatome than the transcriptome. Notably, restricting the database to actual translation products detected by Ribo-Seq provided a deeper coverage of the immunopeptidome than a canonical protein sequence database (Fig. 4h ). This enhanced coverage led to the identification of additional noncHLAIp derived mainly from ORFs that are not included in canonical annotation but still showing periodic footprint of translation, such as those originating in 5ʹ or 3ʹ UTRs, presumed non-coding RNAs, retained introns, and pseudogenes. The majority of those identified were derived from either upstream ORFs or other unannotated ORFs (Fig. 4i-j ). Many of these additional non-cHLAIp were missed using the RNA-Seq inferred database.
Identification of tumor-specific noncHLAIp. As our initial MS search space was not restricted to polypeptide sequences derived from tumor-specific transcripts, we retrospectively investigated the potential of identified noncHLAIp to be classified as tumorspecific. A public database of RNA sequencing data from 30 different healthy tissues (Genotype-Tissue Expression, GTEx  46  ) was assessed at a strict 90 th percentile, which sets the expression of a gene at the top 10% of its expression across all samples. We identified 335 noncHLAIp from 280 lncRNA genes in the seven melanoma samples, of which 23% were expressed in any of our tumor samples and not in the healthy tissues (excluding testis due to its immunoprivileged nature) (Fig. 6 ). Among these genes was the tumor-specific LINC00518 gene, which has been proposed as a two-gene classifier for melanoma detection, together with the TAA PRAME  47.
Using an in-house curated inventory of human TE-derived polypeptide sequences (from three-frame translations) as a reference, we found 88 unique TE-HLAIp in our whole dataset. Some were derived from autonomous TEs, such as long-tandem repeat (LTR) retrotransposons and long interspersed nuclear Fig. 3  MS and ribosome footprint-based evidence of non-canonical peptide generation. A set of proteome-derived tumor-associated antigens, and noncHLAIp (lncRNAs and TEs), from melanoma 0D5P were synthesized in their heavy-labeled form and spiked back into replicates of HLAIp eluted from 0D5P cells to confirm the presence of endogenous HLAIp. The proportions of confirmed and non-confirmed HLAIp as determined by a PRM and b Ribo-Seq-targeted validation are shown for each of the antigen classes. c An example of the co-elution profiles of the transitions of heavy-labeled and endogenous noncHLAIp (from lncRNA; SYLRRHLDF) from 0D5P (left) is shown.
6a ). Importantly, 60 of the 88 TE-HLAIp were found in presumed non-coding TE regions and therefore represent previously unknown HLA peptides. For example, peptides derived from AluSq2 SINE/Alu and L1PA16 LINE/L1 elements were expressed in only skin and testis. These TE-HLAIp would have been overlooked in canonical MS-based searches. We next examined whether our approach could identify tumor-specific non-canonical targets in the ideal case in which normal and tumor biopsies are available, i. e. , from the two lung cancer patient samples included in the present dataset. For the C3N-02671 lung tumor sample, 21 noncHLAp were detected by MS; however, none of the peptides were tumor-specific. In the C3N-02289 sample, we identified 45 noncHLAIp by MS (Fig. 7a ), among which ten peptides were identified uniquely in the tumor tissue.
Four of these source genes were also entirely absent at the RNA level in the adjacent lung (Fig. 7b ). Interestingly, the noncHLAIp from RP11-566H8. 3 was also testis-specific in the GTEx database (90 th percentile transcripts per million (TPM) ≤ 1) (Fig. 7b ), thus qualifying as a non-canonical cancer/ testis antigen. The same analyses of TEs in lung tumor sample C3N-02289 resulted in the identification of one LTR7B LTR/ERV1 TE-HLAIp that was present in the tumor tissue; however, this gene is also expressed in healthy brain. For sample C3N-02671, no TEderived HLAIp were detected. To comparatively assess the expression of canonical tumor antigens in the same samples, we investigated select TAAs using the same methodology (Supplementary Fig. 6b, c ). We identified six TAA protHLAIp that were exclusively detected in the tumor tissue of C3N-02289 (BIRC5, TERT, FAP, SPAG4, MAGEA9, and BCL2L1).
7a, b ). Further, the shared presentation of the noncHLAIp AAF-DRAVHF, derived from the family of LINEs (LINE/L2) on chromosome 6, was confirmed in two melanoma samples (Supplementary Fig. 7c, d ). Interestingly, the corresponding source RNA expression is restricted to the skin and testis. Next, we assessed a large collection of immunopeptidomic datasets (ipMSDB  48  , 91 biological cancer tissue/cell line sources, 35 biological healthy tissues/cell line sources; 1,102 MS raw files in total) and obtained the first large-scale signature of noncHLAIp presentation (Supplementary Data 8). In total, 220,293 peptides were obtained from healthy samples versus 280,385 peptides from cancer samples. We re-identified in ipMSDB 92 tumor-specific noncHLAIp (source genes described above and were expressed at 90th percentile TPM ≤ 1 in a maximum of three tissues) (Fig.
Discussion Our proteogenomics approach led to the stringent identification of hundreds of noncHLAIp derived from presumed non-coding genes, TEs and alternative ORFs. This feat was achieved with NewAnce, a computational module that overcomes the challenge of reduced sensitivity and specificity when searching against large MS search spaces, and it can be applied to any (non-canonical) protein sequence database of interest  28, 50. We rigorously tested the validity of noncHLAIp identifications with HLA binding predictions, sequence-specific hydrophobicity characteristics, targeted MS analyses, and provided evidence of translation in peptide-encoding ORFs by Ribo-Seq. Using all of these strategies together, we confirmed that NewAnce was superior to MaxQuant and Comet alone, across all the investigated samples.
Determination of the correlation between gene expression and HLA-I sampling was assessed by fitting a polynomial curve of degree 3 to each dataset. Pearson correlation values were calculated to assess the correlation between the fitted curve and the corresponding dataset. h With data derived from 0D5P, a comparison of the overall overlap in unique HLAIp identified with RNA-Seq-based and Ribo-Seq-based assembled databases for MS search is shown. i Overlap of noncHLAIp identified by RNA-Seq-and Ribo-Seq-based searches. j The total number of noncHLAIp identified by Ribo-Seq is depicted for each of the respective ORF types. Source data are provided as a Source Data file. detection by MS. Interestingly, the expression and translation of microproteins derived from presumed non-coding RNAs in the heart were recently discovered using a Ribo-Seq directed proteogenomics approach.
Similar to our results, evidence of translation was confirmed for 22. 5% of the lncRNAs, while 55. 4% of the micropeptides were validated by PRM MS  51. Importantly, our results additionally demonstrate that the correct identification of noncHLAIp in proteogenomic workflows requires proper FDR control and validation using multiple independent methods. Combining immunopeptidomics with RNA-Seq and Ribo-Seq datasets enables the comprehensive assessment of how transcription, translation and HLA presentation are correlated. Despite the different methodological challenges  [52] [53] [54]  , we observed the expected correlations between HLA presentation level and expression, especially by Ribo-Seq, presumably because translation is biologically closer to antigen processing and presentation than transcription is.