Genome-wide analysis of dysregulated RNA-binding proteins and alternative splicing genes in keloid

Abstract

Introduction: The pathogenesis of keloids remains unclear. Methods: In this study, we analyzed RNA-Seq data (GSE113619) of the local skin tissue of 8 keloid-prone individuals (KPI) and 6 healthy controls (HC) before and 42 days after trauma from the gene expression omnibus (GEO) database. The differential alternative splicing (AS) events associated with trauma healing between KPIs and HCs were identifified, and their functional differences were analyzed by gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) pathways. The co-expression relationship of differentially alternative splicing genes and differentially expressed RNA binding proteins (RBPs) was established subsequently. Results: A total of 674 differential AS events between the KD42 and the KD0 and 378 differential AS events between the HD42 and the HD0 were discovered. Notably, most of the differential genes related to keloids are enriched in actin, microtubule cells, and cortical actin cytoskeletal tissue pathway. We observed a signifificant association between AS genes (EPB41, TPM1, NF2, PARD3) and trauma healing in KPIs and HCs. We also found that the differential expression of healthy controls-specifific trauma healing-related RBPs (TKT, FDPS, SAMHD1) may affect the response of HCs to trauma healing by regulating the AS of downstream trauma healing-related genes such as DCN and DST. In contrast, KPIs also has specifific differential expression of trauma healing related RBPs (S100A9, HspB1, LIMA1, FBL), which may affect the healing response of KPIs to trauma by regulating the AS of downstream trauma healing-related genes such as FN1 and TPM1. Discussion: Our results were innovative in revealing early wound healing-related genes (EPB41, TPM1, NF2, PARD3) in KPI from the perspective of AS regulated by RBPs.

Keywords: RNA-bindig proteins; RNA-sequencing (RNA-seq); alternative splicing; keloid; trauma healing.

FIGURE 1

Differential gene analysis of local skin repair in patients with KPI and HC after 42 days of a wound. (A) Volcano plots to present all DEGs between HD42 (healthy day) and HD0 samples with DESeq. FDR ≤0.05 and FC (fold change) ≥ 2 or ≤0.5. (B) Volcano plots to present all DEGs between KD42 (keloid day) and KD0 samples with DESeq. FDR ≤0.05 and FC (fold change) ≥ 2 or ≤0.5. (C) Principal component analysis (PCA) based on FPKM value of all DEG. The ellipse for each group is the confidence ellipse. (D) Venn diagram showing the up-regulation overlap between KPI and HC. (E) Venn diagram showing the downregulation overlap between KPI and HC. (F) Bar plot showing the most enriched GO biological process results of the Common up-regulation. (G) Bar plot showing the most enriched GO biological process results of the Common down-regulation. (H) The top 5 most enriched GO terms (biological process) were illustrated for specific genes among the different groups. The color scale shows the row-scaled significance (-log10 corrected p-value) of the terms.

FIGURE 2

Analysis of alternative splicing during local skin repair in KPI and HC after 42 days of a wound. (A) Barplot showing the number of regulatory AS detected by SUVA in each group. ‘alt3p’ indicates the model that 5ʹ splice site is shared and 3ʹ splice site is alternative. “alt5p” indicates the model that 3ʹ splice site is shared and 5ʹ splice site is alternative. “olp” indicates a model that both splice sites are different but part of the splice junction are overlapped. “Contain” indicates a model that both splice sites are different but one splice junction is contained in another splice junction. (B) Splice junction constituting RAS events detected by SUVA was annotated to classical AS event types. And the number of each classical AS event type was shown with a barplot. (C) Barplot showing RAS with different pSAR. RAS which pSAR (Reads proportion of SUVA AS event) ≥ 50% were labeled. (D) Principal component analysis (PCA) based on RAS of pSAR ≥50%. The ellipse for each group is the confidence ellipse. (E) The Heatmap showing the splicing ratio of RAS (PSAR ≥50%). (F) Venn diagram showing overlap of RAS id between KPI (day42 and day0) and HC (day42 and day0) groups. (G) Bar plot showing the most enriched GO biological process results for specific RAS (pSAR ≥50%) in the HC group. (H) Bar plot showing the most enriched GO biological process results for specific RAS (pSAR ≥50%) in the KPI group.

FIGURE 3

Analysis of differential AS of wound healing-associated genes in skin tissue from KPI and HC. (A, B) The Heatmap diagram shows the splicing ratio of wound healing and cytoskeleton in the KPI and HC. (C) Visualization of junction reads distribution of erythrocyte membrane protein band 4.1 (EBP41) in AS events clualt3p13395 from different groups. Splice junctions were labeled with SJ reads the number, and the altered exon was marked out with a red box. The splicing events model is shown in the top panel. Boxplot in the bottom panel showing the splicing ratio profile of the splicing event from EBP41. Boxplot showing splicing ratio of clualt3p13395 EPB41 on the right. *p ≤ 0.05, **p ≤ 0.01,***p ≤ 0.001. (D) Visualization of junction reads distribution of encoding tropomyosin (TPM1) in AS events clualt3p54512 from different groups. Splice junctions were labeled with SJ reads the number, and the altered exon was marked out with a red box. The splicing events model is shown in the top panel. Boxplot in the bottom panel showing the splicing ratio profile of the splicing event from EBP41. Boxplot showing splicing ratio of clualt3p54512 TPM1 on the right. *p ≤ 0.05, **p ≤ 0.01,***p ≤ 0.001.

FIGURE 4

Differentially expressed RBPs potentially affect AS of wound healing-associated genes in KPI and HC. (A) Venn diagram showing the overlap of DEG and RBP. (B) The Heatmap showing the expression profile of HC-specific DERBP. (C) The Heatmap showing the expression profile of KPI-specific DERBP.(D) Bar plot showing the most enriched GO biological process results of RAS co-expressed by specific DERBP in the HC group. Cutoffs of p ≤ 0.01 and Pearson coefficient ≥0.6 or ≤ −0.6 were applied to identify the co-expression pairs. (E) Bar plot showing the most enriched GO biological process results of RAS co-expressed by specific DERBP in the KPI group. Cutoffs of p ≤ 0.01 and Pearson coefficient ≥0.6 or ≤ −0.6 were applied to identify the co-expression pairs. (F) Network diagram showing the wound healing pathway of RAS co-expressed by specific DERBP in HC. (G) Network diagram showing the wound healing pathway of RAS co-expressed by specific DERBP in KPI. The red marked DERBP stands for Pearson coefficient ≥0.6.