dataesr/scientific-paragraphs-categorization huggingface dataset:
Quantitative real-time qRT-PCR was performed using LightCycler 480 SYBR Green Master (Roche) and different sets of primers on LightCycler 480 Real-Time PCR System (Roche). Samples were normalized on GAPDH and analyzed according to the ΔΔCt method. Primer sequences are provided in Supplementary Table  1. Single-cell RNA-Seq analysis. Publicly available single-cell data were used to infer a transcriptome signature of ITGB8-expressing T cells and Treg cells across  different tumor types. scRNAseq counts were downloaded from the GEO repository for melanoma (SKCM: GSE115978), colorectal cancer (CCA: GSE108989), liver cancer (HCC: GSE98638), and non-small cell lung cancer (NSCLC: GSE99254). All single-cell data were analyzed with the "Seurat" package v. 3. 1. 0  40. After creating a seurat object, the sctransform wrapper was used for normalization, scaling, and identification of variable features within each dataset.
UMAP projections were used for the visualization of expression and co-expression patterns. The "FindMarkers" function was used with standard settings to identify genes differentially expressed in ITGB8 pos Tregs (CCA, HCC, and NSCLC). For each set of ITGB8 single-cell markers, the singscore gene signature scoring method was used to score ITGB8 pos Treg activity (Foroutan et al. 41  ). The singscore method uses rank-based statistics to analyze the sample's gene-expression profile and scores the expression activities of gene sets at a single-sample level. TCGA data exploitation. Melanoma (SKCM), colon cancer (COAD), liver cancer (LIHC), and lung cancer (LUAD) gene expression and clinical data were downloaded from the The Cancer Genome Atlas (TCGA) repository [https:// www. cancer. gov/tcga] using the R packages "TCGAbiolinks" v.
2. 9. 4  42  and "RTCGA. clinical" v. 20151101. 8. 0 (Kosinski  43  ), respectively. Briefly, each dataset was queried for Illumina HiSeq RNA-Seq results. Downloaded data were preprocessed, normalized and filtered using "TCGAbiolinks" and edgeR (Robinson et al. 44  ), and a Itgb8 pos score was calculated for each sample, as described above. Tumor samples were classified into two groups, Low and High, based on the median of the Itgb8 pos score. Survival data (i. e. , time to last follow-up and overall survival status), was used to fit a Cox proportional hazards regression model using the "survival" package v. 2. 44-1. 1 [https://CRAN. R-project. org/package=survival]. Survival Kaplan-Meier curves were plotted with "survminer" v. 0. 4. 5 [https:// CRAN. R-project. org/package=survminer]. Additional validations and score correlations were also performed in the melanoma expression dataset  45.
Acknowledgements The expert assistance of  V. Bernet ,  S. Rodriguez ,  E. Guillemot , and  H. Tarayre  is acknowledged as well as the advices and analysis from V. Fontanier. We also thank  Drs. S. Soudja  and  V. Dardalhon  for critical comments and the  Marie lab members  for their helpful discussions. We thank  Dr. Bartholin  for sharing TGFβRI CA mice and  Prof. Ribas  for sharing data on melanoma patients. The help of the core facilities of the  CRCL: PIC (Microscopy) , Flow cytometry, Anican and  B. Manship  for editing is also acknowledged. This work was supported by grants from  LabEx DEVweCAN ANR investissement d'Avenir   ANR-10-LABX-61  (J. C. M. ), the  Helmholtz-DKFZ-Inserm program  (J. C. M. ),  BMS foundation Grant , and Foncer contre le cancer, ARC grant #  2019-1047 , the labelisation ligue Nationale contre cancer (J.
C. M. EL-2016 and  EL-2021 ). Agence Nationale de la Recherche   ANR-13-PDOC-0019  (H. P. ) and  People Program  ( Marie Sklodowska-Curie Actions ) of the  European Union  ( PIIF-GA-2012-330432 ) (H. P. ). A. L. was supported by a  LabEx DEVweCAN  grant and  Ligue Nationale Contre le Cancer. O. L. was supported by the  Ligue Nationale Contre le Cancer  and S. T. by a PhD fellowship from the  French Ministry of Higher Education. The  SIRIC LYriCAN INCa-DGOS-Inserm_12563  is acknowledged.
Data availability Published Single Cell RNA-seq data used in this study was obtained from the Gene Expression Omnibus (GEO) repository under accession codes: GSE115978 (https:// doi. org/10. 1016/j. cell. 2018. 09. 006), GSE108989 (https://doi. org/10. 1038/s41586-018-0694-x), GSE98638 (https://doi. org/10. 1016/j. cell. 2017. 05. 035), and GSE99254 (https:// doi. org/10. 1038/s41591-018-0045-3). Gene expression and survival data were obtained from The Cancer Genome Atlas (TCGA) repository (https://www. cancer. gov/tcga, datasets: SKCM, COAD, LIHC, and LUAD), and the National Center for Biotechnology Information database of Genotypes and Phenotypes (https://www. ncbi. nlm. nih. gov/gap/) with accession number phs001919 (https://doi. org/10. 1038/s43018-019-0003-0). All codes used in this study were based on the following R packages: Seurat v.
3. 1. 0 (doi: Article file 070921. doc), singscore v1. 10. 0 (https://doi. org/10. 1186/s12859-018-2435-4), TCGAbiolinks v. 2. 9. 4 (https://doi. org/10. 1093/nar/gkv1507), RTCGA. clinical v. 20151101. 8. 0, edgeR v. 3. 31. 1 (https://doi. org/10. 1093/bioinformatics/btp616), survival v. 2. 44-1. 1 (https://github. com/therneau/survival), and survminer v. 0. 4. 5 (https:// rpkgs. datanovia. com/survminer/index. html). Data supporting the findings of this study are provided either in supplementary figures or in the Source data file. Source data are provided with this paper.
microRNAs (miRNAs) associate with Ago proteins to post-transcriptionally silence gene expression by targeting mRNAs. To characterize the modes of miRNA-binding, we developed a novel computational framework, called optiCLIP, which considers the reproducibility of the identified peaks among replicates based on the peak overlap. We identified 98 999 binding sites for mouse and human miRNAs, from eleven Ago2 CLIP-seq datasets. Clustering the binding preferences, we found heterogeneity of the mode of binding for different miRNAs. Finally, we set up a quantitative model, named miRgame, based on an adaptation of the game theory. We have developed a new algorithm to translate the miRgame into a score that corresponds to a miRNA degree of occupancy for each Ago2 peak. The degree of occupancy summarizes the number of miRNA-binding sites and miRNAs targeting each binding site, and binding energy of each miRNA::RNA heteroduplex in each peak.
To identify miRNA-binding sites for endogenously expressed miRNAs, we analyzed eleven publicly available Ago2 CLIP-seq datasets from raw data using a computational framework called optiCLIP (optimized pipeline of CLIP-seq data analysis) based on a novel benchmarked workflow for pre-processing and peak calling  (12) , peaks reproducibility assessment among replicates by a newly developed strategy based on Jaccard index calculation, and finally employing miRBShunter prediction program for miRNA-binding sites identification  (4). To evaluate the performance of optiCLIP, we repeated the analysis using other prediction programs, namely, Tar-getScan  (13) , miRanda  (14) , TarPmiR  (15)  and RNAhybrid  (16). We have also set up a quantitative model to study the interplay of miRNA-binding sites, by proposing a degree of occupancy that stratifies Ago2 peaks.
optiCLIP: an integrated and optimized framework for Ago2 CLIP-seq data analysis In order to analyze Ago2 CLIP-seq experiments from raw data files and identify miRNA-binding sites, a computational framework, called optiCLIP, was developed thanks to our recent computational benchmarks of several software used in CLIP-seq analysis workflow  (4, 12). optiCLIP starts with raw data preprocessing that includes adapter removal with cutadapt  (17) , quality filters of reads to identify possible sequencing errors or biases requiring a minimal quality score of 15 and a minimal length of 10 nucleotides, and duplicate collapsing with Prinseq  (18)  with the option '-derep 123' to collapse duplicated reads, '-min qual mean 15 -min len 10' to finally select reads based on quality and length. The alignment of sequence reads to a reference genome is done using Novoalign software (http://www.
The boundaries of the consensus peaks are then calculated by the union of all the peaks that fulfill the reproducibility constrains. To assess the optimal J-index threshold, we tested three different J-index thresholds, namely, the 20%, 40% and 80% (here called J20, J40 and J80, respectively). As last step, optiCLIP identifies the miRNAbinding sites from the identified peaks, using miRBShunter that identifies all potential miRNA::RNA heteroduplexes for both seed and non-seed mode of binding. Heuristic approach was used to adapt the parameters of miRBShunter to datasets with different sizes: where p is the P-value calculated by Homer  (24) , ms is the match score calculated by Homer between the identified motifs and the miRNA sequences, the Fimo score is the threshold for the program FIMO of the MEME suite  (25)  to localize the motifs in the peak.
Specificity in identifying miRNA-binding sites on Ago2 peaks To evaluate the specificity of the identified miRNA-binding sites, we have randomized five times the sequence of the 25 most expressed miRNAs of each dataset by using the Shuffle DNA software (www. bioinformatics. org/ sms2/shuffle dna. html), thus five creating negative control datasets, and used with peaks identified by J20 threshold and the merge method to feed the prediction programs. The number of identified miRNA-binding sites using miRB-Shunter, miRanda and TargetScan by the real sequence of miRNAs (reference dataset) was then compared with each of the negative control datasets (named 'shuffles') Wilcoxon signed-rank test and boxplots were generated by R software.
Precision of miRNA-binding sites The precision of the prediction can be measured by the distance from the miRNA-binding sites and the center of the peaks. The precision has been calculated for miRBShunter, TargetScan and miRanda for all datasets on peaks identified by J20 threshold or merge methods. miRNA-binding sites close to the peak center correspond to high degree of precision  (4, 26). The comparison between the three prediction programs, namely, miRBShunter, miRanda, and Tar-getScan, was performed by using the Kruskal-Wallis test. Boxplots were generated by R software.
Alignment and quantification of the most expressed miRNAs To identify the most expressed miRNAs in the eleven datasets, we selected reads from 18 to 35 nt in length and mapped the reads using Novoalign software with highly stringent parameters: novoalign -c 20 -t 20 -l 5 -s 1 -F STDFQ -o sam -r All -d INDEX GENOME -f input file > output file. Then, we quantified the number of mapped reads for each miRNA, reads count were normalized for each replicate using the quantile normalization. We finally calculated the mean among replicates and ranked them to obtain the 25 most expressed miRNAs for each dataset (Supplemental Table  S1 ).
miRgame: degree of occupancy for Ago2 CLIP-seq peaks The goal of the degree of occupancy is to provide a stratification of the Ago2 peaks, which silencing of the target mR-NAs would be positively correlated. To establish a quantitative model to calculate the degree of occupancy of Ago2 peaks, we made two assumptions: • the identified miRNA-binding sites do not depend on the miRNA expression levels since we identified the binding sites only of the most expressed miRNAs in each dataset (up to the top 25 mostly expressed miRNAs for each dataset); • miRNA-binding sites within the same peak synergistically cooperate. The first assumption implies that we can ignore the concentration levels of each miRNA considered in the analysis. The second assumption implies we can disregard the distance among miRNA-binding sites within the same Ago2 peak.
The marginal contribution score of each coalition for each miRNA-binding site was considered as equal to the relative score calculated by miRanda and the context score for TargetScan, respectively. The use of the logarithmic scale to sum the contribution of each miRNA::RNA heteroduplex within the same binding site would lower down the μ ν score of binding sites targeted by many miRNAs, reasoning that these binding sites would be less regulated by expression changes of single targeting miRNAs, but by the pool of targeting miRNAs, as it has been previously observed  (34). To validate the miRNA degree of occupancy, we used the gene expression profiles whereby single miRNAs were either knocked out or overexpressed. The model miRgame was run on the J20 threshold and merge method peaks and the related binding sites predicted by miRBShunter, miRanda or TargetScan for the considered datasets.
Identification of miRNA-binding sites by J-index method and its sensitivity We collected and analyzed eleven Ago2 CLIP-seq datasets (Table  1 ) generated from human or mouse cells or tissues and analyzed to identify miRNA-binding sites. In Supplemental Figure  S1  is illustrated the computational framework we used for this analysis and named it optiCLIP (optimized pipeline of CLIP-seq data analysis), which is freely downloadable at https://github. com/TrabucchiLab/ optiCLIP (see Materials and Methods for details). We tested three different J-index thresholds, namely, the 20%, 40% and 80% (here called J20, J40 and J80, respectively), by quantifying the number and length of the identified peaks (Figure  1A  and B ). In human datasets, the J20 threshold identified a significantly higher number of peaks in human datasets compared to J40 and J80 thresholds (P < 0.
05, Wilcoxon signed-rank test). In mouse datasets, the number of peaks identified by J20 threshold was higher compared to J40 and J80 thresholds but not significant (P = 0. 062, Wilcoxon signed-rank test) (Figure  1A ). Furthermore, the peaks identified by the J20 threshold had a significantly wider length compared to the two other J-index thresholds (P < 0. 001, Wilcoxon rank sum test) (Figure  1B ). As expected, the number of identified peaks by merge method was higher and the length wider compared to the J20 because less stringent (Figure  1A , B, and Supplemental data 2). Considering the transcriptomic locations of the identified peaks, we observed a significant enrichment of the 3 UTR location in J20 threshold compared to J40 and J80 thresholds and the merge method (P < 0. 01, Wilcoxon signedrank test) (Figure  1C ).
Next, we detected miRNA-binding sites using miRB-Shunter and other four prediction programs namely, Tar-getScan, miRanda, RNAhybrid, and TarPmiR (Supplemental data 2). As shown in Figure  1D , the merge method outperformed the J-index method in number of binding sites, while J20 performed better than the other two J-index thresholds for all prediction programs in human datasets (P < 0. 05, Wilcoxon signed-rank test). Similar results were observed in mouse datasets although the differences among the methods were not significant (P = 0. 062, Wilcoxon signed-rank test). Then, we inspected the transcriptomic location of the identified miRNA-binding sites. As shown in Figure  1E  and Supplemental Figure  S2A, B  and C , miRNA-binding sites were enriched in the 3 UTR in J20 threshold compared to J40 and J80 (Supplemental Table  S2 ), and the merge method for all prediction programs.
Specificity and precision in identifying miRNA-binding sites on Ago2 peaks Concerning the specificity, we performed five sequence randomizations of the top expressed miRNAs in each dataset, identified the miRNA-binding sites in each randomized dataset, and compared them to those found using the original datasets. We used TargetScan, miRanda and miRB-Shunter prediction programs and both J20 threshold and merge methods. Importantly, we found that the number of miRNA-binding sites was significantly much higher in the original datasets compared to the randomized ones (Figure  2A ), indicating that all the three prediction programs are specific. Furthermore, as shown in Figure  2B , all three prediction programs have identified miRNA-binding sites that are close to the center of the related peaks using J20 threshold and merge methods.
Clustering the mode of binding of miRNAs We clustered the data generated by the optiCLIP framework to investigate the miRNA mode of binding in a transcriptome-wide fashion. We plotted the binding sequence position for each miRNA as heatmap of the frequency using the 'hclust' function. Dendrograms were obtained using K-means function. The clustering of miRNA-binding sites was applied to human and mouse datasets using miRBShunter, Tar-getScan, and miRanda prediction programs on peaks identified by J20 threshold and the merge method. As shown in Figure  3A , miRNA-binding sites identified by miRBShunter are focused on the seed sequence just for a subset of miRNAs considering the J20 threshold in human datasets, while in merge method the subset of seed-focused binding sites was smaller (Supplemental Figure  S3A ).
These results indicate that the miRNA mode of binding predicted by miRBShunter show heterogeneity among different miR-NAs. As expected, for TargetScan and miRanda the weight of the miRNA seed in predicting binding sites was much higher respect to the rest of the miRNA sequence for both J20 threshold and the merge method (Supplemental Figures S4A-S7A). In mouse datasets, we observed a variety of mode of binding with no specific subgroup focused on the seed sequence for binding sites identified by miRBShunter in both J20 and merge method (Supplemental Figures  S8A  and S9A ), while for TargetScan and miRanda seed binding was again the enriched mode of binding (Supplemental Figures  S10A-S13A ). Concerning the clustering method, we found similar clustering between hclust and K-means methods for the considered three prediction programs, in both J20 index threshold and merge method for human and mouse datasets (Figure  3 , Supplemental Figures  S3-S13 ).
In Supplemental Figures  S14  and S15 , we plotted the cumulative fractions of the other four miRNAs. Comparing J20 threshold with the merge method, we noticed significant values for miRBShunter considering mmu-let-7a-5p (P = 0. 037, Kolmogorov-Smirnov test) and TargetScan, miRBShunter and RNAhybrid for hsa-miR-124-3p (P = 0. 01, P = 0. 0004 and P = 0. 042, respectively, Kolmogorov-Smirnov test). Comparing either J20 threshold or the merge method with background, mmu-let-7a-5p target repression was significant for miRBShunter (P = 0. 004 and P = 0. 01, respectively, Kolmogorov-Smirnov test), while for TarPmiR and RNAhybrid the merge method outperformed (P = 0. 03 and P = 0. 03, respectively, Kolmogorov-Smirnov test) (Supplemental Figure  S14A ). For hsa-miR1-3p, only the merge method using miRanda program gave significant values (P = 0.
For miRanda these percentages were 5-31% and 6-33% and for TargetScan were 3-27% and 2-32%, re-  Next, we plotted the number of multiple miRNA-binding sites identified by miRBShunter, miRanda and TargetScan per peak as a function of the number of promiscuous miRNA-binding sites for J20 threshold and the merge method (Figure  5C  and Supplemental Figures  S16C, S17C , S18C, S19C and S20C). We found a linear dependency with Spearman correlation score between these two variables for all prediction programs and for both J20 threshold and the merge method. To validate these results, we compared the number of miRNA-binding sites normalized by number of peaks, for each dataset, between the original miRNA sequences and five negative controls (shuffles). The results are shown in Figure  5D  and in Supplemental Figures  S16D, S17D , S18D, S19D and S20D.
Overall, we identified a statistically higher number of miRNA-binding sites per peak in the original datasets compared to each negative control for miRBShunter, TargetScan and miRanda and for both J20 threshold and the merge method (P < 0. 05, Wilcoxon signed-rank test). The same analysis was performed by filtering out the multiple interactions on the promiscuous binding sites and, therefore, considering each promiscuous binding site as one, to avoid the bias toward those miR-NAs that share similar sequences and tend to target the same binding sites in the reference datasets, but likely not in the shuffle controls. As shown in the Supplemental Figure  S21 , we also identified a statistically higher number of miRNA-binding sites per peak in the original datasets in both J20 and merge methods for miRBShunter and mi-Randa.
However, TargetScan did not show any difference between reference and shuffle datasets and mainly found one binding site per peak. Thus, we concluded that only miRBShunter and miRanda specifically and systematically found multiple binding sites. However, TargetScan that applies a much more stringent algorithm mainly found one miRNA-binding site per peak. These data indicate that the identification of multiple miRNA-binding sites was specific and not by chance. Together, these data indicate that on the one hand many Ago2 peaks contain multiple miRNA-binding sites that can potentially interplay to synergistically enhance the miRNAdependent repression. On the other hand, the complexity of the miRNA-dependent silencing of target mRNA is enhanced by the fact that some binding sites can be potentially targeted by many different miRNAs.
miRgame: a quantitative model to calculate the degree of occupancy of Ago2 CLIP-seq peaks To study how the miRNA-binding sites within each Ago2 peak interplay to finely control miRNA-dependent silencing, we derived a model to quantify a miRNA degree of occupancy for each peak. We hypothesize that the degree of occupancy could provide a stratification of Ago2 peaks that would positively correlate to miRNA-dependent mRNA repression. We calculated the degree of occupancy for all peaks containing miRNA-binding sites identified by optiCLIP framework on the 11 Ago2 CLIP-seq datasets of interest and using two additional prediction programs, namely, miRanda and TargetScan, to compare these data. We found a linear distribution of the degree of occupancy for peaks located only in 3 UTR containing miRNA-binding sites when plotted as a function of the number of miRNA-binding sites within each peak (Figure  6A  for J20 threshold and Supplemental Figure  S22A  for the merge method).
Then, we calculated the significance of the downregulation in each level for the mR-NAs containing binding sites for mmu-miR-155-5p, as our best miRNA-mediated repression dataset compared to the other four miRNAs (Figure  4  and Supplemental Figures  S14  and S15 ). The results are presented in Figure  6B  for J20 threshold and in Supplemental Figure  S22B  for the merge method. Concerning the J20 threshold, we found that levels three and four have the strongest target repression for predictions performed by miRBShunter, while for miRanda and TargetScan level four was overall linked with strongest target repression. Similar results were found for the merge method, except for TargetScan that does not show any target repression. We also calculated the significance of the downregulation in each level for the mRNAs containing binding sites for hsa-miR-16-5p, as our second-best miRNA-mediated repression dataset (Figure  4  and Supplemental Figures  S14  and S15 ).
As shown in Supplemental Table  S3 , for this miRNA useful stratification was observed for miRanda merge. For TargetScan J20 significant repression was observed in the level 3. In this latest condition, we further divided the peaks in 2 levels of occupancy degree instead of 4 and found that significant repression was observed in the highest level. In addition, for mmu-miR-155-5p, we compared the repression calculated by miRgame with that by stratifying the dataset with simple counting of miRNA-binding sites per peak (Supplemental Table  S4 ). This analysis indicates that the simple-count method yields similar results of miRgame for miRBShunter and miRanda. However, it fails to stratify e66 Nucleic Acids  Research, 2021, Vol. 49, No. 11  PAGE 12 OF 14 peaks according to the degree of repression for TargetScan.