USE CASE

From Raw Data to Therapeutic Insights: End-to-End RNA-Seq Analysis of IL-17A in Osteoarthritis

Written by Axel Martinelli in collaboration with Almaden Genomics.

Osteoarthritis RNA-Seq data analysis

Introduction

In 2021 Mimpen et al. conducted a study on how different members of the Interleukin-17 (IL-17) cytokine family induce Osteoarthritis-like (OA) gene expression changes. In the original manuscript, the authors exposed human chondrocytes (HC) and synovial fibroblasts (SF) from end-stage OA patients to three members of the IL-17 cytokine family (namely IL-17A, IL-17AF and IL17-F) and observed the dysregulation of pathways associated with angiogenesis, immunity and the complement system. They also observed an association with arthritis and musculoskeletal disease gene sets and concluded that IL-17A played a role in OA pathophysiology.

The biological material (HC and SF cells) for the study was obtained from OA patients. The samples were divided into a control group and three treatment groups (IL-17A, IL-17AF and IL-17F) of six samples each for each cell type and then sequenced, for a total of 24 HC and 24 SF RNA-seq samples.

In collaboration with our partners at Almaden Genomics, we’ve chosen this study as the basis for a practical re-analysis—taking the original dataset from raw RNA-seq data to biological insights using a modern, no-code bioinformatics approach. By combining Almaden Genomics’ g.nome® platform with our own Omics Playground, we aim to reproduce and explore the findings in a transparent, interactive, and fully reproducible workflow.

Pre-processing with g.nome®

The first step is retrieving the RNA-seq data. With Almaden Genomics’ g.nome, public datasets can be imported directly using project accession numbers. This makes the setup fast, reproducible, and user-friendly, so you can focus on the biology rather than the technical details.

Once the data is retrieved, g.nome handles preprocessing end-to-end: from quality control and trimming to alignment and gene-level count generation, all within an intuitive, no-code workflow. The result is a clean, analysis-ready dataset, complete with comprehensive quality reports and organized metadata. For more information about the pre-processing step, refer to Almaden’s blog post.

With the data fully prepared, we can now look into the downstream analysis and interactively explore which biological insights emerge using the Omics Playground platform.

Thanks to the seamless integration, this transition happens with a single click, right from within the g.nome platform.

Clustering Analysis

Samples clustered by cell-line (using the tSNE approach), but beyond that it was impossible to distinguish any further sub-structure in the data according to treatment (Fig. 1).

tSNE plot showing the clustering of the samples according to cell type and treatment. Cell type: green, HC; blue, SF. Treatments: circles, control; squares, IL-17A; stars, IL-17AF; triangles, IL-17F.
Figure 1. tSNE plot showing the clustering of the samples according to cell type and treatment. Cell type: green, HC; blue, SF. Treatments: circles, control; squares, IL-17A; stars, IL-17AF; triangles, IL-17F.

When cell samples were separated by cell type and grouped by treatment, a marked distinction in gene expression between treated (particularly IL-17A treated) and control groups could be observed in clustered heatmaps (Fig. 2A and 3A). In both cases, gene cluster S3 showed the largest difference in gene expression profiles between IL-17A treated samples and the control group. Annotation of the S3 gene cluster in both cell types yielded similar results, displaying a strong correlation with inflammatory responses, angiogenesis, various immune responses and ,in the case of HC samples, complement activation (Fig. 2B and 3B). This was in line with the main observations of the published article. 

Clustered heatmap of the SF cell type samples grouped by treatment. Genes are divided into four main clusters (S1-4).
Figure 2. Clustered heatmap of the HC cell type samples grouped by treatment. Genes are divided into four main clusters (S1-4).
Clustered heatmap of the SF cell type samples grouped by treatment. Genes are divided into four main clusters (S1-4).
Figure 3. Clustered heatmap of the SF cell type samples grouped by treatment. Genes are divided into four main clusters (S1-4).

Differential Expression Analysis

This analysis was performed using the intersection of  three different gene expression analysis methods (EdgeR (qlf), DESEQ2 (wald) and limma (trend)) with FDR<0.05 and logFC>|0.5|.

When comparing IL-17A treated vs control group by cell type, 306 genes (196 up- and 110 down-regulated) were differentially expressed in HC samples but only 38 (32 up- and 6 down-regulated) in SF samples (Fig. 4A and 4B). 29 genes were present in both cell type comparisons (Fig. 5). Among these, IL-6, four C-X-C motif chemokine ligands, NFKB inhibitor zeta (NFKBIZ), C-C motif chemokine ligand 2 (CCL2), and IRF1 were upregulated. All these genes are involved in immune responses. Another upregulated gene was MAP3K8, which is known to be overexpressed in rheumatoid arthritis.

Volcano plots for IL-17A treated vs control pairwise comparisons by cell type. A: HC cell type; B: SF cell type.
Figure 4. Volcano plots for IL-17A treated vs control pairwise comparisons by cell type. A: HC cell type; B: SF cell type.
Venn diagram of the intersection of the dysregulated genes in the IL-17A vs control comparisons by cell type. Genes are separated into up- (top numbers) and down-regulated genes (bottom numbers). A: SF cell type; B: HC cell type.
Figure 5. Venn diagram of the intersection of the dysregulated genes in the IL-17A vs control comparisons by cell type. Genes are separated into up- (top numbers) and down-regulated genes (bottom numbers). A: SF cell type; B: HC cell type.

Gene Set Enrichment Analysis

The expression profiles of the IL-17A vs control groups for the HC and SF cell types were compared against various public datasets using the fgsea algorithm with an FDR<0.05 and logFC>0.2. When tested against the “DISEASE” database in the platform, both comparisons yielded a positive correlation against osteoarthritis, arthritis, rheumatoid arthritis and inflammatory and autoimmune diseases such as scleroderma, Crohn’s disease and discoid lupus erythematosus. Results are shown for the HC cell type profile only (Fig. 6).

Figure 6. Top 10 most correlated disease gene expression profiles with IL-17A treatment in HC cells.
Figure 6. Top 10 most correlated disease gene expression profiles with IL-17A treatment in HC cells.

A GO enrichment analysis showed the significant upregulation (q<0.01) of various gene sets associated with immune responses and inflammation in both HC and SF IL-17A treated samples and to a lesser extent in IL-17AF treated HC samples (Fig. 7). 

Activation matrix generated with Omics Playground
Figure 7. Activation matrix showing the most significant up- (red) and down-regulated (blue) GO terms across all pairwise comparisons. 

Next, I used the drug connectivity panel in the platform to compare the pairwise comparisons against the L1000 drug activity database. I focused on drugs generating significantly (q<0.05) inversely correlated gene expression profiles that could act as potential inhibitors of the IL-17 treatment. The most prominent mode of action identified among potential inhibitors was the inhibition of p38 MAPK (p38 mitogen-activated protein kinases), followed by CDK (cyclin dependent kinases) and MEK (mitogen-activated protein kinase kinases)  inhibition (Fig. 8A). MEKs play a central role in mediating OA pathology and have long been proposed as a therapeutic target. The p38 MAPK pathway has been strongly implicated in the pathology of rheumatoid arthritis and its inhibition suppressed the expression of proinflammatory cytokines in cartilage tissues from OA patients. Finally, CDK9 inhibition has been shown to have beneficial effects both in vitro and in mouse models of OA (1 and 2). 

When looking at individual drugs, I observed a significant negative correlation of all IL-17 treatments with several compounds, of which five stood out (Fig. 8B). The first is fostamatinib, an SYK inhibitor which was trialed on OA patients (3) but appears to have limited efficacy ex vivo (4). Compound AZ-628, a receptor-interacting protein kinase-3 (RIP3) inhibitor, also showed an antagonist profile across the board and has been recently shown to abolish OA pathogenesis (5). Alvocidib was shown to inhibit the development of post-traumatic OA in mice (6), while triptolide has been extensively studied as a treatment for OA (7). Finally, PP-2, an Src kinase inhibitor, has been shown to regulate chondrocyte differentiation, indicating the inhibition of Src kinase as a potential therapeutic approach for skeletal diseases such as OA.

Drug connectivity analysis. A: Most significant drug modes of action (MOA) observed in the drug connectivity map analysis. Negative normalised enrichment scores (NES) indicate MOA that counteract the IL-17 induced profile and can thus have therapeutic applications. The most negatively correlated MOA are highlighted. B: Activation matrix of the most negatively correlated drug profiles across all pairwise comparisons. Three compounds tested in OA patients or OA models are highlighted.
Figure 8. Drug connectivity analysis. A: Most significant drug modes of action (MOA) observed in the drug connectivity map analysis. Negative normalised enrichment scores (NES) indicate MOA that counteract the IL-17 induced profile and can thus have therapeutic applications. The most negatively correlated MOA are highlighted. B: Activation matrix of the most negatively correlated drug profiles across all pairwise comparisons. Three compounds tested in OA patients or OA models are highlighted. 

Biomarkers Discovery

Using the “Find Biomarkers” panel in tha platform, I produced a decision tree to help find genetic markers for different treatment groups (Fig. 9A). While this is only one of the possible decision trees, it indicates how the activation of components of the immune system is a hallmark of OA-like changes. Chemokine CXCL1 is able to distinguish treated and untreated samples and the difference in expression across all treatment groups compared to the control group is pronounced (Fig. 9B). The same is true for the NFKB Inhibitor Zeta gene (NFKBIZ) and, to a lesser extent, CXCL3 and IL-6. In order to distinguish between the different IL-17 treatments, expression levels of CXCL1, CXCL3, IL-6 and NFKBIZ appear to be indicative and point to the greater impact that IL-17A has on gene expression changes than the alternative isoforms. 

Biomarkers discovery. A: One of the possible decision trees for identifying biomarkers distinguishing the various sample groups in the dataset. B: Expression boxplots of the potential biomarker genes identified by the platform. Omics Playground by BigOmics
Figure 9. Biomarkers discovery. A: One of the possible decision trees for identifying biomarkers distinguishing the various sample groups in the dataset. B: Expression boxplots of the potential biomarker genes identified by the platform.

Main Findings

Consistent with the main finding from the original article, the re-analysis of the dataset by Mimpen and colleagues identified OA-like signatures in the transcriptomes of chondrocytes exposed to IL-17A.

The signatures were characterised by the upregulation of genes associated with immune responses and inflammation and matched the public expression profiles of various arthritic and autoimmune conditions. The upregulation of MAP3K8 was also observed. This was correlated to the results from the drug connectivity analysis, which suggested MEK inhibitors as potential therapies for the OA-like profiles. Indeed, MEK inhibition as well as suppression of p38 MAPK  have been tested as treatments for OA, which point to the reliability of IL-17A treatment in mimicking OA transcriptional changes.

Finally, several chemokines and immune-related genes were highlighted as potential biomarkers of IL-17 exposure and also confirmed the greater transcriptional responses induced following IL-17A treatment compared to the IL-17AF and IL-17F isoforms.

Conclusion

With a fully integrated workflow between g.nome and Omics Playground, you can smoothly transition from raw RNA-seq data to high-impact biological insights: all within a streamlined, connected environment. From data preprocessing and QC in g.nome to downstream analysis in Omics Playground, everything flows in just a click.

Discover how g.nome and Omics Playground work together to accelerate your RNA-Seq data analysis and keep your focus where it matters.

Unlock the full potential of your RNA-seq data with g.nome and Omics Playground!

Request a Demo

References

(1) Xue, S., Zhu, L., Wang, C., Jiang, Y., Lu, H., Liu, Y., Shao, Q., Xue, B., Sang, W., & Ma, J. (2019). CDK9 attenuation exerts protective effects on catabolism and hypertrophy in chondrocytes and ameliorates osteoarthritis development. Biochemical and Biophysical Research Communications, 517(1), 132–139. https://doi.org/10.1016/j.bbrc.2019.07.032

(2) Haudenschild, D. R., Carlson, A. K., Zignego, D. L., Yik, J. H., Hilmer, J. K., & June, R. (2019). Inhibition of early response genes prevents changes in global joint metabolomic profiles in mouse post-traumatic osteoarthritis. Osteoarthritis and Cartilage, 27(3), 504–512. https://doi.org/10.1016/j.joca.2018.11.006

(3) Tanaka, Y., Millson, D., & Nakayamada, S. (2020). Safety and efficacy of fostamatinib in rheumatoid arthritis patients with an inadequate response to methotrexate in phase II OSKIRA-ASIA-1 and OSKIRA-ASIA-1X study. Rheumatology, 60(6), 2884–2895. https://doi.org/10.1093/rheumatology/keaa732

(4) Kjelgaard-Petersen, C. F., Christensen, T. R., Hagglund, P., Karsdal, M. A., Thudium, C. S., & Bay-Jensen, A. (2017). AB0062 Ex vivo back-translation of fostamatinib’s effect on joint ecm turnover shows significant effect on bone but no effect on the synovium. https://doi.org/10.1136/annrheumdis-2017-eular.6788

(5) Yang, S., Noh, H. H., Lee, H., Park, H., Ha, Y., Park, S. W., Lee, H., Kim, S. H., Kang, H. W., Eyun, S., Yang, S., & Kim, Y. C. (2020b). TRIM24-RIP3 axis perturbation accelerates osteoarthritis pathogenesis. Annals of the Rheumatic Diseases, 79(12), 1635–1643. https://doi.org/10.1136/annrheumdis-2020-217904

(6) Haudenschild, D. R., Carlson, A. K., Zignego, D. L., Yik, J. H., Hilmer, J. K., & June, R. (2019b). Inhibition of early response genes prevents changes in global joint metabolomic profiles in mouse post-traumatic osteoarthritis. Osteoarthritis and Cartilage, 27(3), 504–512. https://doi.org/10.1016/j.joca.2018.11.006

(7) Lee, H., Huang, X., Son, Y., & Yang, S. (2021). Therapeutic Single Compounds for Osteoarthritis Treatment. Pharmaceuticals, 14(2), 131. https://doi.org/10.3390/ph14020131