Determining Differential Gene Expression of Diabetic Kidney Disease Infected Podocytes for Genes Associated with Cellular Metabolism and Affected Metabolic Pathways

Sam Lerner

Introduction

The original study, which collected the RNA sequencing data used in this analysis, was conducted to identify a statistically significant, differentially expressed long non-coding RNA (lncRNA) in podocytes from three patients infected with diabetic kidney disease (DKD) and three patients with healthy kidneys and determine the role of the lncRNA in the progression of DKD. The authors found that the differentially expressed lncRNA, ENST00000436340, interacts with PTBP1 which initiates binding to RAB3B, which degrades mRNA. This causes cytoskeleton deformation and GLUT4 translocation which is causal to DKD progression. The authors of the paper conducted total RNA sequencing to identify the ENST00000436340 and found 171 total DEmRNAs. This analysis will use the RNAseq data collected from this study to determine metabolic pathway differences between the two groups.

Load the required libraries for this data analysis R markdown script

setwd(dir = "~/project_data")
library(BiocManager)
library(knitr)
library(DESeq2) 
library(RColorBrewer)
library(pheatmap)
library(ggplot2)
library(tidyverse)
library(biomaRt)
library(clusterProfiler) 
library(org.Hs.eg.db) 
library(DOSE) 
library(ggnewscale) 
library(pathview) 
library(ggupset) 
library(enrichplot)
library(data.table)
library(EnhancedVolcano)
library(GO.db)
library(cowplot)

Part 1. Defining the DESeqDataSet for analysis

First, input the annotation file which was created in excel and contains the sample names, file names, conditions, and replicates and make the annotation file a data frame (1). Then characterize condition and replicate as factors and create the DESeqDataSet (dds). Because HTseq-count was used to generate gene counts, the DESeqDataSetFromHTSeqCount function was used to create the dds (2). To reduce file size and improve the quality of figures generated later in the script, any gene with less than 10 reads is filtered out (3). To counteract the default organization of the two conditions, the relevel function was used to establish healthy kidney samples as the reference (4).

# 1
sampleTable <- read.csv("GSE199838_DKDvsHEALTHY_anno.csv", header = T)
sampleTable <- data.frame(sampleTable)

# 2
sampleTable$Condition <- factor(sampleTable$Condition)
sampleTable$Replicate <- factor(sampleTable$Replicate)
 
dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
                                  directory = ".",
                                  design = ~Condition)

# 3
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# 4
dds$Condition <- relevel(dds$Condition, ref = "Healthy")

Part 2. Differential Gene Expression Analysis

The DESeq function calculates the differential gene expression and different statistics for each gene (baseMean, log2FoldChange, lfcSE, stat, pvalue, and padj) and the genes are reordered by smallest padj value (1). Then, the list of DE genes is merged with the normalized gene counts and the resulting table is saved (2). This saved file is recalled multiple times throughout the script. Finally, the data set is filtered for only differentially expressed genes with log2FC > 1 and padj < 0.05 and the resulting dataset is saved(3).

# 1
dds <- DESeq(dds)
res <- results(dds, contrast=c('Condition', 'DKD', 'Healthy')) 
result <- res[order(res$padj), ]
table(result$padj<0.05)
## 
## FALSE  TRUE 
## 38746   487
# 2
resdata <- merge(as.data.frame(result), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE)
names(resdata)[1] <- "Gene_id"
write.csv(resdata, file="DKDvsHEALTHY_unfiltered_matrix.csv")

# 3
resdata <- as.data.frame(result) %>% dplyr::filter(abs(log2FoldChange) > 1 & padj < 0.05) %>% tibble::rownames_to_column("Gene_id")
write.csv(resdata, "DKDvsHEALTHY_log2fc1_fdr05_DEGlist.csv", row.names = TRUE)

Part 3. Pathway Analysis

To set up the unfiltered matrix for pathway analysis, import ensembl id and entrez id, and match the two imports (1). Then merge the annotated gene list with resdata by the gene IDs(2). To prepare the background gene set, specify for log2FoldChange, name the gene list, remove any values stored as NA, and sort the list in decreasing order of log2FoldChange (3). To prepare the significant gene set, select for only genes with padj < 0.05, specify for log2FoldChange, name the gene list, remove any values stored as NA, and filter to allow for only genes with an abolute log2FoldChange > 1 (4). Call for the human database from BioMart and conduct enrichGO analysis, output the analysis results, and save the data frame (5). Create a bar plot representing the most differentially expressed pathways from the enrichGO analysis and save both plots (6). Create a cnet plot to show the genes that are associated with the GO Terms (7).

# 1
pathway_resdata <- read.csv("DKDvsHEALTHY_unfiltered_matrix.csv", header = T)
pathway_resdata <- pathway_resdata[, c(2:14)]
ensembl <- useEnsembl(biomart = "genes",
                      dataset = "hsapiens_gene_ensembl",
                      mirror = "useast")
biomart.output <- getBM(mart = ensembl, values = pathway_resdata$Gene_id, 
                      filters = "ensembl_gene_id_version",
                      attributes = c("ensembl_gene_id_version", 
                                     "entrezgene_id"))

# 2
names(biomart.output)[1] <- "Gene_id"
resdata_gene <- dplyr::left_join(pathway_resdata, biomart.output,
                        by=c('Gene_id'='Gene_id'))

# 3
original_gene_list <- resdata_gene$log2FoldChange
names(original_gene_list) <- resdata_gene$entrezgene_id
gene_list<-na.omit(original_gene_list)
gene_list = sort(gene_list, decreasing = TRUE)

# 4
sig_genes_df = subset(resdata_gene, padj < 0.05)
genes <- sig_genes_df$log2FoldChange
names(genes) <- sig_genes_df$entrezgene_id
genes <- na.omit(genes)
genes <- names(genes)[abs(genes) > 1]

# 5
keytypes(org.Hs.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MAP"         
## [16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"        
## [21] "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [26] "UNIPROT"
go_enrich <- enrichGO(gene = genes,
                      universe = names(gene_list),
                      OrgDb = org.Hs.eg.db,
                      keyType = 'ENTREZID',
                      readable = T,
                      ont = "BP",
                      pvalueCutoff = 0.05,
                      qvalueCutoff = 0.05)
cluster_summary <- data.frame(go_enrich)
write.csv(cluster_summary, "clusterProfiler_DKDvsHEALTHY.csv")

# 6
cluster_barplot <- barplot(go_enrich, showCategory = 10,
                   title = "GO Biological Pathways",
                   font.size = 8)
ggsave("DKDvsHealthy_all-ClusterBarplot.png")
cluster_barplot


Figure 1. A cluster barplot displaying the ten most representative GO Terms, selected by the highest gene ratio (number of associated genes/ total number of significant genes). The p-adjusted values for each GO Term is shown by the color with red being padj = 0.0025 and blue being padj = 0.0125.

# 7
cnetplot(go_enrich, showCategory = 4,
         categorySize="pvalue", foldChange=gene_list, 
         circular = TRUE, colorEdge = TRUE,
         cex_label_gene = 0.4, vertex.label.font=6,
         cex_label_category = 0.5)

ggsave("DKDvsHealthy_all-ClusterBarplot.png")


Figure 2. A circular cnet plot that shows the relationship between the differentially expressed genes and their GO terms. The different colored lines correspond to the GO term and the color of the gene points represents the fold change of the gene expression.

Part 4. Examination of Metabolic Pathways

The first line of this code chunk establishes the organism (homo sapien). The gseKEGG analysis uses the KEGG database which contains genes associated with specific pathways and calculates the pathway association of the differentially expressed genes (1). The generated kk2 list is used to determine which pathways are most significantly different between DKD and normal cells and the pathways related to metabolism are selected (2). A GSEA plot is made with the selected metabolic pathways (3).

# 1
kegg_organism = "hsa"
kk2 <- gseKEGG(geneList     = gene_list,
               organism     = kegg_organism,
               nPermSimple  = 10000,
               minGSSize    = 3,
               maxGSSize    = 800,
               pvalueCutoff = 0.05,
               pAdjustMethod = "none",
               keyType       = "kegg") #I’m not sure if it takes “ENTREZID” or requires geneid
write.csv(kk2, "KEGG_DKDvsHEALTHY.csv")

# 2
kk2_met_pathways <- kk2 %>%
  filter(Description == 'Nucleotide metabolism' | Description == 'Oxidative phosphorylation' | Description == 'Purine metabolism' | Description == 'Pyrimidine metabolism')

# 3
finalgsea <- gseaplot2(kk2_met_pathways, geneSetID = 1:4)
finalgsea


Figure 3. A Gene Set Enrichment Analysis Plot of four different metabolic pathways. The Running Enrichment Score (y-axis) tracks the increases and decreases of DKD differential expression for each gene associated with the labeled pathways ranked in order of most significantly different to least (x-axis). The actual Enrichment Score for each pathway is the peak on the running Enrichment Score. The vertical lines below the Running Enrichment Scores represent each gene in the corresponding pathway. The bottom of the plot, Ranked List Metric, shows how the genes are distributed between the DKD and Healthy phenotypes with DKD being on the left side and Heathy on the right side.

Part 5. Creation of Gene List of Differentially Expressed Metabolism Associated Genes

Import the data, use “getBM” to label each gene with its name, add the names to the resdata matrix, and organize the resdata matrix so the gene name is the first column (1). Use the “useEnsembl” function to select the gene dataset (2). Use the “GOBOFFSPRING” function from the GO.db library to create a list of GO ids in the cellular metabolism GO family and use this list to create a list of the genes associated with these GO terms (3). Filter the “resdata_named” dataset by the metabolism GO gene list so only the genes associated with cellular metabolism are in the “resdata_met” dataset (4).

# 1
resdata_unfiltered <- read.csv("DKDvsHEALTHY_unfiltered_matrix.csv", header = T)
biomart.output <- getBM(mart = ensembl, values = resdata_unfiltered$Gene_id, 
                      filters = "ensembl_gene_id_version",
                      attributes = c("external_gene_name",
                                     "ensembl_gene_id_version"))
resdata_named <- dplyr::left_join(resdata_unfiltered, biomart.output,
                        by=c('Gene_id'='ensembl_gene_id_version'))
resdata_named <- resdata_named[, c(15, 2:14)]

# 2
ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")

# 3
metabolism_and_offspring <- c(GOBPOFFSPRING[["GO:0044237"]], "GO:0044237")
ensembl_m_and_off <- getBM(attributes=c('ensembl_gene_id_version'),
                   filters = 'go', values = metabolism_and_offspring, mart = ensembl)

# 4
resdata_named<- data.table(resdata_named)
met_genes <- data.table(ensembl_m_and_off)
resdata_met <- subset(resdata_named, Gene_id %in% met_genes$ensembl_gene_id_version)


Part 6. Generation of Volcano Plots for All Genes and Metabolism Associated Genes

This section generates a volcano plot for all the genes and another for only the metabolism associated genes, so consider both datasets when reviewing the code. Create a threshold variable defined by having a statistically significant p-adjusted value (< 0.05) and count the number of statistically significant (p-adjusted value < 0.05) differentially expressed genes (log2FC > 1) in both data sets (1). Create volcano plots with an FC cutoff of 1 and a p-adjusted cutoff of 0.05 (2).

# 1
threshold_padj_all <- resdata_named$padj < 0.05
threshold_padj_met <- resdata_met$padj < 0.05
resdata_named$threshold <- threshold_padj_all
resdata_met$threshold <- threshold_padj_met

resdata_all_sig_DEG <- as.data.frame(resdata_named) %>% dplyr::filter(abs(log2FoldChange) > 1 & padj < 0.05)
resdata_met_sig_DEG <- as.data.frame(resdata_met) %>% dplyr::filter(abs(log2FoldChange) > 1 & padj < 0.05) 

resdata_all_up <- as.data.frame(resdata_all_sig_DEG) %>% dplyr::filter(log2FoldChange > 1)
resdata_all_down <- as.data.frame(resdata_all_sig_DEG) %>% dplyr::filter(log2FoldChange < -1)
resdata_met_up <- as.data.frame(resdata_met_sig_DEG) %>% dplyr::filter(log2FoldChange > 1)
resdata_met_down <- as.data.frame(resdata_met_sig_DEG) %>% dplyr::filter(log2FoldChange < -1) 

# 2
EnhancedVolcano(resdata_named,
                lab = NA,
                x = 'log2FoldChange',
                y = 'padj', 
                title = "DKD versus Healthy (all genes)",
                subtitle = "Differential expression", 
                FCcutoff = 1,
                pCutoff = 0.05,
                axisLabSize = 12, 
                col=c('black', 'grey', 'black', 'red3'),
                legendLabels=c('Not sig.','Log2FC','padj', 'padj & Log2FC'),
                legendPosition = 'right',
                legendLabSize = 16,
                legendIconSize = 5.0,
                gridlines.major = FALSE,
                gridlines.minor = FALSE,
                ylim = c(0,10))

ggsave("DKDvsHealthy_all-Enhancedvolano.png")

Figure 4. A volcano plot showing the statistically significant (p-adjusted value < 0.05) differentially expressed (Log2FC > 1) of all the genes in the dataset. The points in red in the upper right portion of the graph are significant, upregulated genes in DKD cells and the points in red in the upper left are significant, downregulated genes in DKD cells. There are 487 total significant differentially expressed genes, 304 downregulated and 183 upregulated. The gene RNF5 had a -log2P of >30 and a log2FC of ~-5, but the y limits cut this out to make the graph easier to read.

EnhancedVolcano(resdata_met,
                lab = as.character(resdata_met$external_gene_name),
                x = 'log2FoldChange',
                y = 'padj', 
                selectLab = c('CCL19', 'AGTR2', 'NFKBIL1', 'HBB', 'NEUROG1'),
                title = "DKD versus Healthy (Metabolism Associated Genes)",
                subtitle = "Differential expression", 
                FCcutoff = 1,
                pCutoff = 0.05,
                labSize = 4, 
                axisLabSize = 12, 
                col=c('black', 'grey', 'black', 'red3'),
                labCol = 'black',
                labFace = 'bold',
                boxedLabels = TRUE,
                legendLabels=c('Not sig.','Log2FC','padj', 'padj & Log2FC'),
                legendPosition = 'right',
                legendLabSize = 16,
                legendIconSize = 5.0,
                drawConnectors = TRUE,
                widthConnectors = 1.0,
                colConnectors = 'black', 
                gridlines.major = FALSE,
                gridlines.minor = FALSE,
                ylim = c(0,10))

ggsave("DKDvsHealthy_met-Enhancedvolano.png")

Figure 5. A volcano plot of genes associated with metabolic pathways. This uses the same layout as Figure 4, except this only includes the genes in the cellular metabolism GO family. The labeled genes are the genes most significantly and differentially upregulated or downregulated in DKD cells. There are 87 total significant differentially expressed metabolism associated genes, 56 downregulated and 31 upregulated. Like in Figure 4, the gene RNF5 had an extremely low p-adjusted value and was cutoff by the y-axis limits to make the graph easier to read.


Part 7. Examination of the Most Significant Differentially Expressed Metabolism Associated Genes

The “expression_plot” data set is created by organizing the sample columns all in the same column under the name “name” and the data set is merged via the “left_join” function with “sampleTable” by “name” (1). Then for each of the labeled genes in Figure 5, as well as RNF5, an expression data set is made by filtering the “expression_plot” data set by name of the gene (2). For the six genes of interest, the “Condition” column is converted into a factor type and the levels are defined for organization of the graphs. Then individual box plots for each gene is made with the normalized gene counts on the y-axis and grouped by sample condition (3). The gene count box plots are displayed together in grid format (4). Finally, a table is made with brief descriptions of the six genes (5).

# 1
expression_plot <- pivot_longer(resdata_met,
                                cols = 9:14,
                                values_to = "normalized_counts")
colnames(sampleTable)[1] = "name"
expression_plot <- left_join(x = expression_plot, 
                             y = sampleTable, 
                             by = "name")

# 2
CCL19_exp <- expression_plot %>%
  filter(external_gene_name == "CCL19")
AGTR2_exp <- expression_plot %>%
  filter(external_gene_name == "AGTR2")
NFKBIL1_exp <- expression_plot %>%
  filter(external_gene_name == "NFKBIL1")
RNF5_exp <- expression_plot %>%
  filter(external_gene_name == "RNF5")
HBB_exp <- expression_plot %>%
  filter(external_gene_name == "HBB")
NEUROG1_exp <- expression_plot %>%
  filter(external_gene_name == "NEUROG1")

# 3
CCL19_exp$Condition <- factor(CCL19_exp$Condition, levels = c("DKD", "Healthy"))
CCL19_plot <- ggplot(CCL19_exp) +
  geom_boxplot(aes(x=Condition, 
                   y=normalized_counts, 
                   fill=Condition)) + 
  theme_bw() + scale_fill_manual(values = c("deepskyblue1", "firebrick3")) +
        ggtitle("CCL19") +
        xlab("Condition") + 
        ylab("Normalized Counts") +
        theme(legend.position = "none",
              panel.grid = element_blank(), 
              plot.title = element_text(size = rel(1.25), hjust = 0.5), 
              axis.title = element_text(size = rel(0.75)),
              axis.text = element_text(size = rel(0.75))) 

AGTR2_exp$Condition <- factor(AGTR2_exp$Condition, levels = c("DKD", "Healthy"))
AGTR2_plot <- ggplot(AGTR2_exp) +
  geom_boxplot(aes(x=Condition, 
                   y=normalized_counts, 
                   fill=Condition)) + 
  theme_bw() + scale_fill_manual(values = c("deepskyblue1", "firebrick3")) +
        ggtitle("AGTR2") +
        xlab("Condition") + 
        ylab("") +
        theme(legend.position = "none",
              panel.grid = element_blank(), 
              plot.title = element_text(size = rel(1.25), hjust = 0.5), 
              axis.title = element_text(size = rel(0.75)),
              axis.text = element_text(size = rel(0.75)))

NFKBIL1_exp$Condition <- factor(NFKBIL1_exp$Condition, levels = c("DKD", "Healthy"))
NKFBIL1_plot <- ggplot(NFKBIL1_exp) +
  geom_boxplot(aes(x=Condition, 
                   y=normalized_counts, 
                   fill=Condition)) + 
  theme_bw() + scale_fill_manual(values = c("deepskyblue1", "firebrick3")) +
        ggtitle("NFKBIL1") +
        xlab("Condition") + 
        ylab("") +
        theme(legend.position = "none",
              panel.grid = element_blank(), 
              plot.title = element_text(size = rel(1.25), hjust = 0.5), 
              axis.title = element_text(size = rel(0.75)),
              axis.text = element_text(size = rel(0.75)))

RNF5_exp$Condition <- factor(RNF5_exp$Condition, levels = c("DKD", "Healthy"))
RNF5_plot <- ggplot(RNF5_exp) +
  geom_boxplot(aes(x=Condition, 
                   y=normalized_counts, 
                   fill=Condition)) + 
  theme_bw() + scale_fill_manual(values = c("deepskyblue1", "firebrick3")) +
        ggtitle("RNF5") +
        xlab("Condition") + 
        ylab("Normalized Counts") +
        theme(legend.position = "none",
              panel.grid = element_blank(), 
              plot.title = element_text(size = rel(1.25), hjust = 0.5), 
              axis.title = element_text(size = rel(0.75)),
              axis.text = element_text(size = rel(0.75))) 

HBB_exp$Condition <- factor(HBB_exp$Condition, levels = c("DKD", "Healthy"))
HBB_plot <- ggplot(HBB_exp) +
  geom_boxplot(aes(x=Condition, 
                   y=normalized_counts, 
                   fill=Condition)) + 
  theme_bw() + scale_fill_manual(values = c("deepskyblue1", "firebrick3")) +
        ggtitle("HBB") +
        xlab("Condition") + 
        ylab("") +
        theme(legend.position = "none",
              panel.grid = element_blank(), 
              plot.title = element_text(size = rel(1.25), hjust = 0.5), 
              axis.title = element_text(size = rel(0.75)),
              axis.text = element_text(size = rel(0.75))) 

NEUROG1_exp$Condition <- factor(NEUROG1_exp$Condition, levels = c("DKD", "Healthy"))
NEUROG1_plot <- ggplot(NEUROG1_exp) +
  geom_boxplot(aes(x=Condition, 
                   y=normalized_counts, 
                   fill=Condition)) + 
  theme_bw() + scale_fill_manual(values = c("deepskyblue1", "firebrick3")) +
        ggtitle("NEUROG1") +
        xlab("Condition") + 
        ylab("") +
        theme(legend.position = "none",
              panel.grid = element_blank(), 
              plot.title = element_text(size = rel(1.25), hjust = 0.5), 
              axis.title = element_text(size = rel(0.75)),
              axis.text = element_text(size = rel(0.75))) 
# 4
boxplot_grid <- plot_grid(CCL19_plot,
                          AGTR2_plot,
                          NKFBIL1_plot,
                          RNF5_plot,
                          HBB_plot,
                          NEUROG1_plot,
                          ncol = 3,
                          rel_widths = c(1, 1, 1),
                          labels = NA)
ggsave("DKDvsHealthy_met-boxplots.png")
boxplot_grid


Figure 6. A grid of box plots displaying the normalized gene counts of six different significant differentially expressed genes associated with metabolic pathways. The first row of plots show upregulated genes in DKD cells and the bottom row of genes shows downregulated genes in DKD cells.

# 5
biomart.met_genes <- getBM(mart = ensembl, values = c("CCL19", "AGTR2", "NFKBIL1", "RNF5", "HBB", "NEUROG1"), 
                      filters = "external_gene_name",
                      attributes = c("external_gene_name",
                                     "description"))
biomart.met_genes
##   external_gene_name
## 1            NFKBIL1
## 2               RNF5
## 3                HBB
## 4              AGTR2
## 5              CCL19
## 6            NEUROG1
##                                                         description
## 1          NFKB inhibitor like 1 [Source:HGNC Symbol;Acc:HGNC:7800]
## 2         ring finger protein 5 [Source:HGNC Symbol;Acc:HGNC:10068]
## 3        hemoglobin subunit beta [Source:HGNC Symbol;Acc:HGNC:4827]
## 4  angiotensin II receptor type 2 [Source:HGNC Symbol;Acc:HGNC:338]
## 5 C-C motif chemokine ligand 19 [Source:HGNC Symbol;Acc:HGNC:10617]
## 6                   neurogenin 1 [Source:HGNC Symbol;Acc:HGNC:7764]

Table 1. Brief descriptions of the six metabolism associated genes of interest.

Session Info

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] cowplot_1.1.1               GO.db_3.16.0               
##  [3] EnhancedVolcano_1.16.0      ggrepel_0.9.3              
##  [5] data.table_1.14.8           enrichplot_1.18.4          
##  [7] ggupset_0.3.0               pathview_1.38.0            
##  [9] ggnewscale_0.4.8            DOSE_3.24.2                
## [11] org.Hs.eg.db_3.16.0         AnnotationDbi_1.60.2       
## [13] clusterProfiler_4.6.2       biomaRt_2.54.1             
## [15] lubridate_1.9.2             forcats_1.0.0              
## [17] stringr_1.5.0               dplyr_1.1.2                
## [19] purrr_1.0.1                 readr_2.1.4                
## [21] tidyr_1.3.0                 tibble_3.2.1               
## [23] tidyverse_2.0.0             ggplot2_3.4.2              
## [25] pheatmap_1.0.12             RColorBrewer_1.1-3         
## [27] DESeq2_1.38.3               SummarizedExperiment_1.28.0
## [29] Biobase_2.58.0              MatrixGenerics_1.10.0      
## [31] matrixStats_0.63.0          GenomicRanges_1.50.2       
## [33] GenomeInfoDb_1.34.9         IRanges_2.32.0             
## [35] S4Vectors_0.36.2            BiocGenerics_0.44.0        
## [37] knitr_1.42                  BiocManager_1.30.20        
## 
## loaded via a namespace (and not attached):
##   [1] shadowtext_0.1.2       fastmatch_1.1-3        systemfonts_1.0.4     
##   [4] BiocFileCache_2.6.1    plyr_1.8.8             igraph_1.4.2          
##   [7] lazyeval_0.2.2         splines_4.2.2          BiocParallel_1.32.6   
##  [10] digest_0.6.31          yulab.utils_0.0.6      htmltools_0.5.5       
##  [13] GOSemSim_2.24.0        viridis_0.6.2          fansi_1.0.4           
##  [16] magrittr_2.0.3         memoise_2.0.1          tzdb_0.3.0            
##  [19] Biostrings_2.66.0      annotate_1.76.0        graphlayouts_0.8.4    
##  [22] timechange_0.2.0       prettyunits_1.1.1      colorspace_2.1-0      
##  [25] blob_1.2.4             rappdirs_0.3.3         textshaping_0.3.6     
##  [28] xfun_0.39              crayon_1.5.2           RCurl_1.98-1.12       
##  [31] jsonlite_1.8.4         graph_1.76.0           scatterpie_0.1.9      
##  [34] ape_5.7-1              glue_1.6.2             polyclip_1.10-4       
##  [37] gtable_0.3.3           zlibbioc_1.44.0        XVector_0.38.0        
##  [40] DelayedArray_0.24.0    Rgraphviz_2.42.0       scales_1.2.1          
##  [43] DBI_1.1.3              Rcpp_1.0.10            viridisLite_0.4.1     
##  [46] xtable_1.8-4           progress_1.2.2         gridGraphics_0.5-1    
##  [49] tidytree_0.4.2         bit_4.0.5              httr_1.4.5            
##  [52] fgsea_1.24.0           pkgconfig_2.0.3        XML_3.99-0.14         
##  [55] farver_2.1.1           sass_0.4.5             dbplyr_2.3.2          
##  [58] locfit_1.5-9.7         utf8_1.2.3             labeling_0.4.2        
##  [61] ggplotify_0.1.0        tidyselect_1.2.0       rlang_1.1.0           
##  [64] reshape2_1.4.4         munsell_0.5.0          tools_4.2.2           
##  [67] cachem_1.0.7           downloader_0.4         cli_3.6.1             
##  [70] generics_0.1.3         RSQLite_2.3.1          gson_0.1.0            
##  [73] evaluate_0.20          fastmap_1.1.1          ragg_1.2.5            
##  [76] yaml_2.3.7             ggtree_3.6.2           bit64_4.0.5           
##  [79] tidygraph_1.2.3        KEGGREST_1.38.0        ggraph_2.1.0          
##  [82] nlme_3.1-160           KEGGgraph_1.58.3       aplot_0.1.10          
##  [85] xml2_1.3.3             compiler_4.2.2         rstudioapi_0.14       
##  [88] filelock_1.0.2         curl_5.0.0             png_0.1-8             
##  [91] treeio_1.22.0          tweenr_2.0.2           geneplotter_1.76.0    
##  [94] bslib_0.4.2            stringi_1.7.12         highr_0.10            
##  [97] lattice_0.20-45        Matrix_1.5-1           vctrs_0.6.2           
## [100] pillar_1.9.0           lifecycle_1.0.3        jquerylib_0.1.4       
## [103] bitops_1.0-7           patchwork_1.1.2        qvalue_2.30.0         
## [106] R6_2.5.1               gridExtra_2.3          codetools_0.2-18      
## [109] MASS_7.3-58.1          withr_2.5.0            GenomeInfoDbData_1.2.9
## [112] parallel_4.2.2         hms_1.1.3              grid_4.2.2            
## [115] ggfun_0.0.9            HDO.db_0.99.1          rmarkdown_2.21        
## [118] ggforce_0.4.1

References

Blighe K, Rana S, Lewis M (2022). EnhancedVolcano: Publication-ready volcano plots with enhanced colouring and labeling. R package version 1.16.0, https://github.com/kevinblighe/EnhancedVolcano.

Hu, J., Wang, Q., Fan, X. et al. Long noncoding RNA ENST00000436340 promotes podocyte injury in diabetic kidney disease by facilitating the association of PTBP1 with RAB3B. Cell Death Dis 14, 130 (2023). https://doi.org/10.1038/s41419-023-05658-7

Love, M.I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8