{'fraction_non_alphanumeric': 0.035247359062859536, 'fraction_numerical': 0.022016525468047277, 'mean_word_length': 5.28425895497864, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transcriptome signature of cellular senescence.

Cellular senescence, an integral component of aging and cancer, arises in response to diverse triggers, including telomere attrition, macromolecular damage and signaling from activated oncogenes. At present, senescent cells are identified by the combined presence of multiple traits, such as senescence-associated protein expression and secretion, DNA damage and β-galactosidase activity; unfortunately, these traits are neither exclusively nor universally present in senescent cells. To identify robust shared markers of senescence, we have performed RNA-sequencing analysis across eight diverse models of senescence triggered in human diploid fibroblasts (WI-38, IMR-90) and endothelial cells (HUVEC, HAEC) by replicative exhaustion, exposure to ionizing radiation or doxorubicin, and expression of the oncogene HRASG12V. The intersection of the altered transcriptomes revealed 50 RNAs consistently elevated and 18 RNAs consistently reduced across all senescence models, including many protein-coding mRNAs and some non-coding RNAs. We propose that these shared transcriptome profiles will enable the identification of senescent cells in vivo, the investigation of their roles in aging and malignancy and the development of strategies to target senescent cells therapeutically.

Progress to identify senescent cells in order to exploit them therapeutically has been hampered by a lack of robust and universal measurable traits. Thus far, senescence has been studied in a range of cell types induced by diverse triggers such as replicative exhaustion, DNA damage, oxidation and other stress conditions like signaling through oncoproteins. Due to this heterogeneity, finding broad biomarkers of senescence has been challenging (16,48) and senescent cells are currently found through the combined detection of multiple biochemical markers such as p16, p53, p21 and SA-␤Gal, despite the fact that they are not exclusively nor consistently induced in senescence (17)(18)(19)(20).
Proliferating (P) WI-38 and IMR-90 fibroblasts were used at population doubling level (PDL)15-PDL25, and for replicative senescence (S), cells were cultured until replicative exhaustion (PDL50-PDL59). For ionizing radiation (IR)-induced senescence, proliferating (PDL25) WI-38 and IMR-90 cells were exposed to 10 Gy, and HAECs and HUVECs to 4 Gy; cells were harvested 10 days after IR. For Doxorubicin (Dox)-induced senescence, WI-38 cells (PDL25) were treated with 2 g/ml Dox (Sigma) for 24 h and harvested 7 days later. For oncogene-induced senescence (OIS), WI-38 cells (PDL25) were transduced at 10 MOI (multiplicity of infection) using lentivirus expressing HRAS G12V or an empty vector (EV) and treated with puromycin (1 g/ml) for 5 days. Three additional senescence comparison groups (WI-38 fibroblasts that were either proliferating or rendered senescent by replicative exhaustion, Dox treatment or exposure to IR) were generated in-house using the same workflow, library preparation, genome alignment and normalization methods, P-value adjustment and threshold cutoff as the above comparison groups. Senescence-associated ␤-galactosidase (SA-␤Gal) activity was assessed using a Senescence ␤-galactosidase staining kit (Cell Signaling).

RNA isolation and sequencing
Total RNA was extracted using TriPure (Roche), and the original six senescence models (HUVEC IR+/IR-, HAEC IR+/IR-, IMR-90 IR+/IR-, IMR-90 S/P, WI-38 Dox+/Dox-and WI-38 HRAS G12V /EV) were sequenced by Genewiz. Paired-end cDNA libraries were created after rRNA depletion using Illumina Ribo-Zero and sequenced using an Illumina Hiseq 4000 with a 2× 150-bp configuration single index per lane. Library quality was guaranteed with 80% of bases having a Q-score >30. The remaining three models of senescence, WI-38 S/P, WI-38 IR+/IR-and WI-38 Dox+/Dox-, were prepared by RNA extraction and a paired-end cDNA library created after rRNA depletion using Illumina Ribo-Zero; they were sequenced in-house using the Illumina Hiseq 2500 with a 2 × 150-bp configuration, single index per lane. The original six groups were sequenced at a depth of 150 million paired-end reads, and the latter three groups were sequenced at a depth of 220 million paired-end reads, with minimum of 75 million unidirectional reads for the former and 110 million unidirectional reads for the latter, respectively. RNA-seq analyses of all senescence models are in Supplementary Tables S1-S9 (GSE130727).

Bioinformatic Analysis
Aligned gene counts were analyzed using the edgeR package (Bioconductor release version 3.8) in R (version 3.5.1). Trimmed Mean of M values (TMM) was used as a scaling normalization method for the differential expression analysis of the RNA-seq data. Common dispersion and tagwise dispersion were estimated using the quantile-adjusted conditional maximum likelihood (qCML) method. A pairwise in-group plot of the counts per million (CPM) using mapped transcripts after normalization served as quality control for sample bias and abundance bias (Sup 2A-2K). Benjamini-Hochberg adjusted P-value (false discovery rate or FDR) < 0.15 was used as the significance threshold for a given individual group comparison after considering lower thresholds that were too permissive and higher thresholds that yielded too few shared transcripts. The Venn diagram of the multiple comparison overlaps and subsets was generated using the Vennerable package in R. Ingenuity Pathway Analysis (winter 2018 version) was conducted using the maximum fold change for those genes found to be consistently elevated or reduced across all senescence models.

RT-qPCR analysis
Total RNA was used for reverse transcription (RT) using random hexamers and reverse transcriptase (Invitrogen); real-time quantitative (q)PCR analysis was then performed with SYBR Green master mix (Kapa Biosystems) and genespecific primers using Applied Biosystems 7300 instrument. Forward and reverse primers used to validate the RNA-seq results are listed in Supplementary Table S10.
Senescence was monitored by assessing senescenceassociated ␤-galactosidase (SA-␤Gal) activity; senescent cells displayed the characteristic flattened morphology and enhanced SA-␤Gal activity in all models, except for HU-VECs, where the staining was negative after IR under the conditions tested ( Figure 1C). Senescence was also characterized by western blot analysis of one or several protein markers [p21 (CDKN1A), p16 (CDKN2A), p53 (TP53)] in both WI-38 and IMR-90 cells following replicative exhaustion, and in OIS-and Dox-induced senescence in WI-38 cells ( Figure 1A and B). RNA reverse transcription (RT) followed by real-time quantitative (q)PCR analysis revealed increased p21 mRNA levels in IMR-90 (S) and in WI-38 OIS relative to the respective controls ( Figure 1A and B). Additional markers in these paradigms of cell senescence have been shown elsewhere (21), including higher p21 mRNA levels in IR-and Dox-induced senescence in WI-38 fibroblasts, and IR-induced senescence of IMR-90 fibroblasts, HAECs and HUVECs. These well-established models of cellular senescence were used for further analysis.

RNA expression is most strongly influenced by the cell origin
In order to identify shared features of the transcriptomes of senescent cells, we designed six comparison groups using four cell lines for RNA sequencing, as shown in the workflow (Figure 2A). In this analysis, we compared replicative (WI-38, IMR-90), IR-induced (WI-38, IMR-90, HAEC, HUVEC), Dox-induced (WI-38) and oncogeneinduced (WI-38) senescence. After whole-cell RNA extraction from proliferating and senescent cells, libraries were prepared. RNA-seq performed using an Illumina should be HiSeq 4000 instrument at a depth of 75 million unidirectional reads for a total of 150 million paired-end reads per sample; after completion, bioinformatic analyses were performed (Supplementary Tables S1-S9 and Figures S1 and  S2).
Pearson correlation coefficient analyses to compare the individual sequenced samples against one another showed that the bulk differences in RNA expression were depen-dent on the cell type of origin, in this case endothelial cells (HUVECs, HAECs) or fibroblasts (WI-38, IMR-90), rather than on the method used to induce senescence ( Figure 2B). This finding is consistent with the diversity of senescent programs depending on the tissue of origin and supports the view that cellular senescence is not a single phenotypic state (16).

Shared differentially expressed RNAs
RNA-seq analysis was followed by identification of differentially expressed transcripts meeting the threshold of FDR P < 0.15 in each comparison and further comparing all six groups to assess commonly downregulated RNAs. This FDR threshold was chosen after first considering more stringent thresholds that rendered either no transcripts or too few transcripts for subsequent analysis, as well as less stringent thresholds that yielded longer lists of transcripts that failed subsequent validation. A Venn diagram represents the 112 RNA transcripts significantly downregulated that overlap among the comparisons ( Figure 3A). We sought to reduce possible technical bias and increase the number of cell models tested in order to have the most stringent possible list of RNAs. To this end, we included three additional senescence comparison groups generated using WI-38 cells that were either proliferating or rendered senescent by replicative exhaustion, Dox treatment or exposure to IR ('Materials and Methods' section). After merging all the datasets, the list of shared downregulated RNAs was reduced to 18 transcripts ( Figure 3B); most of them were protein-coding transcripts, including MCUB, PTMA, HIST1H1D and HIST1H1A mRNAs. Only one non-coding RNA was downregulated in all comparisons, lincRNA AC074135.1.
Using the same strategy, we uncovered 251 transcripts that increased during senescence ( Figure 4A), and further analysis with the external three datasets (as above) revealed 50 additional transcripts common across all senescence models ( Figure 4B). Among these, mitochondrial DNAencoded tRNAs (tRNA-TA, -TN, -TC, -TY) (22), two were lncRNAs and three were antisense transcripts ( Figure 4B). However, the majority were protein-coding transcripts, including SRPX, SRPX 2 and STAT1 mRNAs. Together, the data in Figures 3 and 4 identified shared upregulated (50) or downregulated (18) RNAs in the indicated senescence models (16). A heat map representation of these 68 RNAs (50+18) is shown in Supplementary Figure S3. Ingenuity pathway analysis (IPA) revealed several metabolic and inflammatory pathways such as granzyme A signaling, cysteine biosynthesis and interferon signaling (Supplementary Figure S4A). Network analysis found that several of the encoded proteins were part of shared regulatory networks in the nucleus and cytoplasm (Supplementary Figure S4B).

Validation of shared transcripts
The changes in shared downregulated and upregulated transcripts (Figures 3 and 4) were validated by using RT-qPCR analysis employing gene-specific primer pairs (Supplementary Table S10). Among the many RNAs shared across all senescence models, we validated 16 reduced and    12 increased transcripts ( Figure 5A and B) in replicative, IR-induced, Dox-induced and oncogene-induced senescence in WI-38 cells, as well as in replicative and IR-induced senescence in IMR-90 cells and IR-induced senescence in HUVEC and HAEC; additional validation in biological replicates is shown in Supplementary Figure S5. We included p16 and p21 mRNAs as canonical controls during the validation of individual transcripts. Interestingly, however, while p16 mRNA was consistently elevated in most senescent populations, this increase was marginal in Dox-induced and IR-induced WI-38 cell senescence. On the other hand, p21 mRNA was elevated in WI-38 senescence, but was minimally or not at all changed in senescent IMR-90 fibroblasts or in HUVECs or HAECs. By contrast, the individually tested RNAs showing differential expression across all senescence models were successfully confirmed. For instance, FAM129A, HIST1H1E, HIST1H1D, HIST1H1A and ITPRIPL1 mRNAs all showed reduced abundance ( Figure 5A); interestingly, the lncRNA PURPL was not only consistently upregulated in every senescence model, but also vastly surpassed the increases seen for positive controls p16 and p21 mRNA ( Figure 5B). Together, these findings suggest that these 68 transcripts may serve as superior RNA biomarkers of senescence rather than the classical p16 and p21 mRNAs.

Shared transcripts contribute to defining a senescence signature
Next, we sought to determine if the 68 differentially expressed RNAs might help to discriminate between senescent and non-senescent cells. Unbiased principal component analysis (PCA) using the normalized log CPM for all the mapped genes in the sequenced samples revealed a clustering pattern similar to that established by the Pearson Coefficient heat map ( Figure 2B), whereby the majority of the genetic variance lays upon the cell type of origin ( Figure  6A). Importantly, however, when the PCA was recalculated including only those 68 transcripts, we observed a clear separation based on senescent versus non-senescent cell state, in addition to cell type ( Figure 6B).
We were interested in finding the minimum number of transcripts, among the 68, needed to achieve a discriminatory pattern between senescent and non-senescent cells using the RNA-seq log CPM values. The computational intensity and practical limitations of permutating all the unique combinations of 68 transcripts (∼10 20 possible combinations) prompted us to use statistical modeling to reach this goal. Supervised logarithmic (Binomial) regression with LASSO (least absolute shrinkage and selection operator) regularization was employed, assigning non-senescent samples a truth value of 0 and senescent samples a truth value of 1. LASSO is a powerful machine learning and statistical method that is applicable for performing dimensional reduction of data for downstream analysis and may also help in feature selection and in revealing hidden biological information. Using leave-one-out (LOO) crossvalidation to evaluate performance, we were able to decrease the number of transcripts to 5 at an optimal value of λ (penalty coefficient), while still observing a low misclassification error rate (Supplementary Figure S6). The five tran-scripts are SLCO2B1, CLSTN2 and PTCHD4 mRNAs, as well as LINC02154 and PURPL lncRNAs. We then recalculated the PCA including CPM values of only those five RNA transcripts and found that the separation axis between senescent and non-senescent samples was preserved ( Figure 6C), while there was a decrease in the separation based on cell of origin. It is also important to note that while these five transcripts were computationally selected as the best variables for this statistical modeling approach, their specific biological relevance in senescence remains to be investigated.

DISCUSSION
We report global profiles of RNAs expressed in eight experimental models of cellular senescence commonly studied using WI-38 and IMR-90 fibroblasts, HUVECs and HAECs (Figure 1). We identified 50 elevated and 18 reduced transcripts shared among all comparison pairs (Figures 3-5). Computational analyses including Pearson correlation coefficient and PCA indicated that the senescence profiles are more dependent on the cell of origin than on the mode of induction (Figures 2 and 6). For example, IR-induced senescence in WI-38 and IMR-90 fibroblasts is similar and closer to senescence triggered by other stimuli in fibroblasts than it is to IR-induced endothelial senescence (Figure 6). The dissimilarity between senescent fibroblasts and endothelial cells could be due to differences in the dynamic progression of these cell types toward senescence, a dimension of senescence that has not been studied systematically to-date and warrants investigation. The transcripts listed in CellAge (23) were largely (>90%) shared with the aggregate of all differentially expressed mRNAs in at least one cell type and a given senescence trigger (Supplementary Tables  S1-S9).
Interestingly, histone-encoding HIST1H1D, HIST1H1A, HIST1H1E and HIST2H2AB mRNAs were less abundant in senescent cells, in keeping with earlier evidence that the general expression of histones (38)(39)(40), and in particular those encoded by these mRNAs (41)(42)(43), is reduced in senescent cells. These data are in agreement with an epigenetic role for histones in the induction and maintenance of senescence. Indeed, H2AFJ is required for the production of SASP factors by fibroblasts and old human epidermis and was previously found to accumulate in senescent cells (44).
The lncRNA PURPL (p53-upregulated regulator of p53 levels) was one of the transcripts most highly and consistently elevated among senescent cells (Figure 4 and 5). PURPL production was recently shown to be transcrip- tionally controlled by the transcription factor p53, which is generally elevated in senescent cells. However, PURPL in turn helps to maintain low basal p53 levels in colon cancer cells (HCT116) by interacting with MYBBP1A and thereby preventing the assembly of the stable p53-MYBBP1A complex (45). In addition, given that PURPL promotes tumorigenicity by acting as a pro-survival factor in these cells (45), it is possible that PURPL may also contribute to the pro-survival phenotype of senescent cells. Our systematic analysis further reveals that PURPL may be a more robust marker transcript than p16 mRNA or p21 mRNA, as the latter two mRNAs did not reach significance cutoffs in all senescence models, as determined by RNA-seq and RT-qPCR analyses.
In closing, our study highlights the value of computational and statistical modeling methods, not simply as practical analytical tools, but also as a means to identify latent structures hidden beyond high-dimensional biological data such as transcriptomes. The subsets of transcripts we have identified, coding and non-coding, display shared expression patterns across a range of senescent cell models. Given their robust and consistent expression patterns, we propose that they may be involved in the establishment and/or maintenance of senescence, although it remains to be studied whether they regulate senescence and by what mechanisms. Finally, while cultured models of cellular senescence recapitulate the phenotypic traits of senescent cells in tissues and organs, future analyses will focus on determining whether these RNAs are also differentially expressed in physiologic and pathologic settings in vivo. Their potential utility to target senescent cells in 'senolytic' interventions will also be the goals of upcoming studies. With rising appreciation that senescent cells modulate age-associated processes (16,46), and that key pathways active in senescent cells (e.g. immune and lysosomal function) are also active in aging cells (47), the senescence-associated transcripts identified here constitute informative molecular parameters in studies of aging health and disease.