Gene Expression Analysis of RNA-Seq Data from Human Osteoarthritis Knee

Abstract

Osteoarthritis has always been one of the most researched degenerative disorders pervasively. Moreover, this is been classified under complex disorders with an emphasis on the clinical analysis for the older people. However, a genomic understanding of the hereditary impact of this disease seems lacking. The basic genomic analysis has identified various biomarkers, but lacks prognosis. In this direction, we have retrieved publicly available RNA- seq data from the NCBI Sequence Read Archive (SRA). Utilising the Galaxy browser, the data analysis was performed to identify the differentially expressed genes from SRP IDs

{ERP105501, SRP322624}. The overlaps of these studies revealed 105 upregulated and 127 downregulated genes. These sets were further processed for multi-omic analysis through the databases. We identified specific process involvement through the Gene Ontology followed by hub genes through Protein -protein interactions. In addition to this, the transcription factors were identified from the list through BART search and miRNAs interacting with the identified hub genes through miRNet. This study aims to perform genomic analysis of osteoarthritis by identifying the differentially expressed genes among multiple tissues from human osteoarthritis knee joint and further analyse to find the key genes, biological pathways, transcription factors and miRNAs that can be considered as hub genes of the Osteoarthritis disease in humans.

1. Introduction

2. Materials And Methods

3. Results And Discussion

3.1 DATA CURATION

From the search results, 16 samples [8 osteoarthritis (biological replicates) and 8 non- osteoarthritis] collected from the Posterior Lateral Condyle of knee joint (ERP105501) and 6 samples [3 osteoarthritis and 3 non-osteoarthritis] collected from the Synovial Tissue of knee joint (SRP322624) were selected and imported into galaxy browser.

The transcriptome response of osteoarthritis in cartilage and possible patient subgroups were investigated using large-scale RNA-Seq analysis. The posterior lateral condyle of the knees of 60 patients with osteoarthritis (OA) who underwent complete knee replacement surgery and 10 control non-OA patients who underwent amputation were studied for data collection.

Pigmented villonodular synovitis (PVNS) is a rare disorder that causes significant joint damage and has a high recurrence rate even after surgical resection. It is caused by benign synovial tissue growth. However, a lack of knowledge about the pathophysiology prevents effective treatment. overall layout: In order to examine the PVNS transcriptome, the gene expression patterns of three patients with PVNS and three patients with osteoarthritis (OA) were investigated using integrated RNA sequencing (RNA-seq) and microarray.

Table 4.1.1 Metadata of the selected sample

The samples were further pre-processed for quality improvement in fastp and considered for mapping. Reads with reduced quality features such as very low quality, very short read length were removed. Adapter sequences were automatically detected by fastp for paired end data and were removed. Mismatching base pairs in overlapping region of paired ends were also removed. From the fastQC output data, all samples had a sequence duplication levels close to 70%. But, the number of duplicate reads were substantially higher than what can be inferred from read coverage data. This problem is even worse for RNA-seq data, which are frequently dominated by the transcripts of a small number of genes. For RNA-seq data, the exceptionally high read coverage for the specific highly expressed transcripts can easily result in fastQC read duplication levels of 70% or greater [14].

3.2 MAPPING OF READS TO THE GENOME

The output of HISAT2 were very significant as on an average, above 98% of all reads across all samples were mapped to the genome. The human genome GRCh38/hg38 was used as the reference genome for mapping the reads. The index file for this genome version was available as a built-in in the HISAT2 tool and was used for quick and sensitive mapping.

Table 4.2.1 List of all samples' read mapping levels

3.3 COUNTING READS PER GENES

The featureCounts output summarized on an average of 78% assignment of reads across all samples. Of those Assigned reads, 98% of reads were mapped exactly once to the reference genome. Since only proportions under 70% should be checked for potential contamination, we can move on to the analysis. Less than 20% of the reads from all the samples were mapped to multiple locations on the reference genome. For next-generation sequencing technology short reads, this falls within the expected range.

3.4 DIFFERENTIALLY EXPRESSED FEATURES

The DeSeq2 tool was used to perform estimation on the differential expression of genes across all samples. As a result, we arrived at 61,498 unique lines both from Synovial tissue samples and PLC samples that can identify unique genes and also the expression levels of the gene by various parameters.

Internal normalisation is carried out by DESeq2 in which the geometric mean for each gene across all samples is determined. Then, this mean is divided by the gene counts in each sample. The size factor for a sample is the median of these ratios in that sample. Estimates of dispersion have a direct relationship with variance and and inverse relationship with mean. The dispersion is higher for small mean counts and smaller for large mean counts based on this connection. The variance of the genes will be the only factor affecting the dispersion estimates for genes with the same mean. As a result, for a given mean value, the dispersion estimates represent the variance in gene expression. Figure 4.5.1best describe the data fit to the model and thus the data we consider is proven fit to perform differential expression analysis.

Not all the lines out of DeSeq2 shall be statistically sound. It is important to screen the list to potential entries only. One of the columns in the output is the Fold change in log2. The input sequence of factor levels is crucial because the log2 fold changes are based on the principal factor level 1 vs factor level 2 comparison. In this instance, DESeq2 calculates fold changes of "Osteoarthritis" samples against "normal" from the first factor "Osteoarthritis," where the values represent up- or downregulation of genes in Osteoarthritis samples. p-value is the measure of the statistical significance of this change, therefore we filter the entry lines for significant p-value<0.05. Out of 61,498 lines, only 30,375 (the rest being scaffold genes in the chr.gtf file which was used because it doesn’t include unplaced or unlocalized contigs, but only assembled chromosomes) were annotated as genes belonging to Homo Sapiens. Further applying the foresaid filters, 3,305 lines (5.37%) were kept in Synovial data and 4,941 lines (8.03%) were kept in PLC data. Now, we will only choose genes with a fold change (FC) of greater than 2 or less than 0.5. We filter for "(abs(log 2 FC) > 1" (which denotes FC > 2 or FC 0.5) because the DESeq2 output file contains "(log 2 FC)," rather than FC itself. From 3,305 lines, this filter brought down the lines to 2394 (72.44%) in Synovial data and 2,239 (45.31%) in PLC data.

All the 2394 and 2239 gene IDs from the final list were annotated with reference to the human genome GRCh38/hg38 and the gene name was appended to each line. Now we have a list of unique genes differentially expressed in the sample classes. Out of 2394 genes in OA Synovial tissues, 1159 genes were estimated to be upregulated and 940 were downregulated. Likewise, in OA PLC tissues, 662 were upregulated and 1087 were downregulated.

Venn diagrams [15] were drawn to find the overlapping genes between the two data classes and arrived at 105 overlapping upregulated genes and 127 overlapping downregulated genes as seen in Figure 4.5.3. The list of these genes was extracted for further analysis based on upregulated and downregulated genes.

3.5 GENE ONTOLOGY

Gene ontology enrichment analysis for biological process listed that for the upregulated gene set, 20 biological processes were implicated while for the downregulated gene set, as large as 269 GO terms were implicated. Given the percentage of genes in the entire genome that are assigned to a given GO dataset, the P-value is the likelihood of the query genes from the reference dataset. In other words, the background distribution of annotation is compared to the GO keywords shared by the genes in the query list. The p-value was set to 95% confidence and was transformed -log10() form and the corresponding GO terms with the most significant adjusted p-values were screened. From the sorted list of GO terms with high significance, GO terms that come under specific growth and development processes were selected as biomarkers (figure 4.6.1 and figure 4.6.2) since osteoarthritis is a continuous degenerative disease.

3.6 PREDICTION OF TRANSCRIPTION FACTORS

BART displayed 918 Transcription factors regulating the expression of both the gene sets of our study. 4 proteins in the downregulated gene set and 5 in the upregulated gene set act as transcription factor in our OA phenotype and their significance is validated with their corresponding z-scores.

Table 4.7.1 Transcription Factors among our gene sets regulating expression of key genes as Biomarkers

3.7 PROTEIN-PROTEIN INTERACTION ANALYSIS

22 interactions among the upregulated gene set and 105 interactions among the downregulated gene set were found from the PPI network analysis using STRING. Out of these FOS(4), GRIA2(4), ERBB4(4) were the highly interacted proteins and thus considered as hub genes, similarly in the downregulated gene set MELK(11), ASF1B(10), FOXM1(9) considered as hub genes.

4. Conclusion

Declarations

Conflict of Interest: None

Author’s Contribution:  KK designed the work . NR and EE carried out the work.

Acknowledgement: We sincerely acknowledge the bioinformatics facility provided by Bioinformatics department, SASTRA University.

References

1.            Fernandez-Moreno, M., Rego, I., Carreira-Garcia, V., & Blanco, F., 2008. Genetics in Osteoarthritis. Current Genomics, 9(8), 542-547

2.            Valdes, A., & Spector, T., 2008. The Contribution of Genes to Osteoarthritis. Rheumatic Disease Clinics of North America, 34(3), 581-603

3.            Williams, C. J., & Jimenez, S. A., 1993. Heredity, genes and osteoarthritis. Rheumatic diseases clinics of North America, 19(3), 523–543

4.            Paras Hospitals., 2022. Retrieved 8 July 2022, from https://www.parashospitals.com/blogs/osteoarthritis-knee

5.            Afgan, E., Nekrutenko, A., Grüning, B., Blankenberg, D., Goecks, J., & Schatz, M. et al., 2022. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2022 update. Nucleic Acids Research, 50(W1), W345–W351

6.            Leinonen, R., Sugawara, H., & Shumway, M., 2010. The Sequence Read Archive. Nucleic Acids Research, 39(Database), D19-D21

7.            Chen, S., Zhou, Y., Chen, Y., & Gu, J., 2018. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890

8.            Kim, D., Langmead, B., & Salzberg, S. L., 2015. HISAT: a fast spliced aligner with low memory requirements. Nature Methods, 12(4), 357–360

9.            Liao, Y., Smyth, G. K., & Shi, W., 2013. featureCounts: an efficient general-purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7), 923–930

10.         Love, M. I., Huber, W., & Anders, S., 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12)

11.         Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, Thomas PD. PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res. Jan 2019, 47(D1), D419–D426

12.         Ma, W., Wang, Z., Zhang, Y., Magee, N., Feng, Y., & Shi, R. et al., 2021. BARTweb: a web server for transcriptional regulator association analysis. NAR Genomics and Bioinformatics, 3(2)

13.         Wang, Z., Civelek, M., Miller, C., Sheffield, N., Guertin, M., & Zang, C., 2018. BART: a   transcription   factor   prediction   tool   with   query   gene   sets    or    epigenomic profiles. Bioinformatics, 34(16), 2867–2869

14.         Why does FASTQC show unexpectedly high sequence duplication levels (PCR-duplicates)? DNA                            Technologies      Core.,    2022.     Retrieved             10           July               2022,     from https://dnatech.genomecenter.ucdavis.edu/faqs/why-does-fastqc-show-unexpectedly-high- sequence-duplication-levels-pcr-duplicates/

15.         Heberle, H., Meirelles, G., da Silva, F., Telles, G., & Minghim, R., 2015. InteractiVenn: a web-based tool for the analysis of sets through Venn diagrams. BMC Bioinformatics, 16(1)

16.         Chang, L., Yao, H., Yao, Z., Ho, K., Ong, M., & Dai, B. et al., 2021. Comprehensive Analysis of Key Genes, Signaling Pathways and miRNAs in Human Knee Osteoarthritis: Based on Bioinformatics. Frontiers In Pharmacology, 12

17.         Soul, J., Dunn, S., Anand, S., Serracino-Inglott, F., Schwartz, J., Boot-Handford, R., & Hardingham, T., 2017. Stratification of knee osteoarthritis: two major patient subgroups identified by genome-wide expression analysis of articular cartilage. Annals Of the Rheumatic Diseases, 77(3), 423-423

18.         Zhao, Y., Lv, J., Zhang, H., Xie, J., Dai, H., & Zhang, X. (2021). Gene Expression Profiles Analyzed Using Integrating RNA Sequencing, and Microarray Reveals Increased Inflammatory Response, Proliferation, and Osteoclastogenesis in Pigmented Villonodular Synovitis. Frontiers In Immunology, 12