#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# install.packages("tidyverse")
# install.packages("gplots")
#BiocManager::install("edgeR")
#BiocManager::install("tximport")
#BiocManager::install("tximportData")
#BiocManager::install("readr")
#BiocManager::install("GenomicFeatures")
#BiocManager::install("txdbmaker")
#BiocManager::install("org.Mm.eg.db")
library(tximport)
library(readr)
library(GenomicFeatures)
library(dplyr)
library(org.Mm.eg.db)
library(biomaRt)
library(org.Hs.eg.db)
library(edgeR)
library(tidyverse)
library(gplots)
library(ggridges)
library(ggplot2)
library(enrichplot)
library(clusterProfiler)
Make .RDS tximport file
# Example: List of file paths to quant.sf files
files <- c("SRR29804826/quant.sf", "SRR29804827/quant.sf", "SRR29804828/quant.sf",
"SRR29804829/quant.sf", "SRR29804830/quant.sf", "SRR29804831/quant.sf")
names(files) <- c("SRR29804826", "SRR29804827", "SRR29804828",
"SRR29804829", "SRR29804830", "SRR29804831")
list(files)
## [[1]]
## SRR29804826 SRR29804827 SRR29804828
## "SRR29804826/quant.sf" "SRR29804827/quant.sf" "SRR29804828/quant.sf"
## SRR29804829 SRR29804830 SRR29804831
## "SRR29804829/quant.sf" "SRR29804830/quant.sf" "SRR29804831/quant.sf"
# Navigate to tx2gene file GRCm39.104
# setwd("C:/Users/micha/OneDrive/Desktop/GTF_File")
tx2gene_mouse_GRCm39.104 <- read.csv("tx2gene_mouse_GRCm39.104.csv")
head(tx2gene_mouse_GRCm39.104)
## seqnames start end width strand TXNAME GENEID
## 1 1 3143476 3144545 1070 + ENSMUST00000193812 ENSMUSG00000102693
## 2 1 3172239 3172348 110 + ENSMUST00000082908 ENSMUSG00000064842
## 3 1 3322980 3323459 480 + ENSMUST00000192857 ENSMUSG00000102851
## 4 1 3536810 3583776 46967 + ENSMUST00000161581 ENSMUSG00000089699
## 5 1 3602018 3602943 926 + ENSMUST00000192183 ENSMUSG00000103147
## 6 1 3750378 3752011 1634 + ENSMUST00000193244 ENSMUSG00000102348
# This is important because tximport expects the input for the tx2gene argument to be a two-column data # # frame where the first column contains transcript IDs and the second column contains gene IDs.
tx2gene_mouse_GRCm39.104 <- tx2gene_mouse_GRCm39.104[, c("TXNAME", "GENEID")]
head(tx2gene_mouse_GRCm39.104)
## TXNAME GENEID
## 1 ENSMUST00000193812 ENSMUSG00000102693
## 2 ENSMUST00000082908 ENSMUSG00000064842
## 3 ENSMUST00000192857 ENSMUSG00000102851
## 4 ENSMUST00000161581 ENSMUSG00000089699
## 5 ENSMUST00000192183 ENSMUSG00000103147
## 6 ENSMUST00000193244 ENSMUSG00000102348
ID Conversion
# Sanity check
head(rosetta_txi_output_remove_dup, 5)
## ENSEMBL ENTREZID SYMBOL GENENAME UNIPROT
## 1 ENSMUSG00000000001 14679 Gnai3 G protein subunit alpha i3 A2AE36
## 6 ENSMUSG00000000003 54192 Pbsn probasin F8SH42
## 10 ENSMUSG00000000028 12544 Cdc45 cell division cycle 45 F8WJ72
## 17 ENSMUSG00000000037 107815 Scml2 Scm polycomb group protein like 2 B1AVB5
## 22 ENSMUSG00000000049 11818 Apoh apolipoprotein H Q01339
# Sanity check
head(matched.rosetta_txi_output_remove_dup,5)
## ENSEMBL ENTREZID SYMBOL GENENAME UNIPROT
## 1 ENSMUSG00000000001 14679 Gnai3 G protein subunit alpha i3 A2AE36
## 6 ENSMUSG00000000003 54192 Pbsn probasin F8SH42
## 10 ENSMUSG00000000028 12544 Cdc45 cell division cycle 45 F8WJ72
## 17 ENSMUSG00000000037 107815 Scml2 Scm polycomb group protein like 2 B1AVB5
## 22 ENSMUSG00000000049 11818 Apoh apolipoprotein H Q01339
Exploratory Analysis
Performing EDA on tximport results helps to ensure data quality and
prepares the dataset for more reliable differential expression analysis.
This step is integral to identifying and addressing potential issues
that could affect downstream results.
# Load raw txi results
txi <- readRDS("tximport_results.RDS")
txi_counts <- txi$counts
# Load in sample metadata
metadata <- read.csv(file = "metadata.csv",
header = TRUE,
row.names = 1)
head(txi_counts, 5)
## SRR29804826 SRR29804827 SRR29804828 SRR29804829 SRR29804830
## ENSMUSG00000000001 860.058 752.006 891.041 451.121 552.919
## ENSMUSG00000000003 0.000 0.000 0.000 0.000 0.000
## ENSMUSG00000000028 48.001 51.000 36.002 57.000 88.001
## ENSMUSG00000000037 82.992 72.957 106.194 39.012 63.003
## ENSMUSG00000000049 1.000 0.000 1.000 0.000 1.000
## SRR29804831
## ENSMUSG00000000001 489.607
## ENSMUSG00000000003 0.000
## ENSMUSG00000000028 58.000
## ENSMUSG00000000037 58.001
## ENSMUSG00000000049 1.000
head(metadata, 5)
## Assay.Type AvgSpotLen Bases BioProject BioSample
## SRR29804826 RNA-Seq 300 6091030500 PRJNA1135118 SAMN42466452
## SRR29804827 RNA-Seq 300 6068516100 PRJNA1135118 SAMN42466453
## SRR29804828 RNA-Seq 300 6634182600 PRJNA1135118 SAMN42466454
## SRR29804829 RNA-Seq 300 6955524600 PRJNA1135118 SAMN42466455
## SRR29804830 RNA-Seq 300 7595194500 PRJNA1135118 SAMN42466456
## Bytes cell_type Center.Name Collection_Date Consent
## SRR29804826 1954304553 cardiomyocytes DAPING HOSPITAL missing public
## SRR29804827 1958679557 cardiomyocytes DAPING HOSPITAL missing public
## SRR29804828 2157785695 cardiomyocytes DAPING HOSPITAL missing public
## SRR29804829 2277133369 cardiomyocytes DAPING HOSPITAL missing public
## SRR29804830 2487797898 cardiomyocytes DAPING HOSPITAL missing public
## DATASTORE.filetype DATASTORE.provider
## SRR29804826 fastq,run.zq,sra gs,ncbi,s3
## SRR29804827 fastq,run.zq,sra gs,ncbi,s3
## SRR29804828 fastq,run.zq,sra gs,ncbi,s3
## SRR29804829 fastq,run.zq,sra gs,ncbi,s3
## SRR29804830 fastq,run.zq,sra gs,ncbi,s3
## DATASTORE.region Experiment genotype
## SRR29804826 gs.us-east1,ncbi.public,s3.us-east-1 SRX25304157 WT
## SRR29804827 gs.us-east1,ncbi.public,s3.us-east-1 SRX25304156 WT
## SRR29804828 gs.us-east1,ncbi.public,s3.us-east-1 SRX25304155 WT
## SRR29804829 gs.us-east1,ncbi.public,s3.us-east-1 SRX25304154 WT
## SRR29804830 gs.us-east1,ncbi.public,s3.us-east-1 SRX25304153 WT
## geo_loc_name_country geo_loc_name_country_continent geo_loc_name
## SRR29804826 uncalculated uncalculated missing
## SRR29804827 uncalculated uncalculated missing
## SRR29804828 uncalculated uncalculated missing
## SRR29804829 uncalculated uncalculated missing
## SRR29804830 uncalculated uncalculated missing
## Instrument Library.Name LibraryLayout LibrarySelection
## SRR29804826 Illumina NovaSeq 6000 GSM8393301 PAIRED cDNA
## SRR29804827 Illumina NovaSeq 6000 GSM8393300 PAIRED cDNA
## SRR29804828 Illumina NovaSeq 6000 GSM8393299 PAIRED cDNA
## SRR29804829 Illumina NovaSeq 6000 GSM8393298 PAIRED cDNA
## SRR29804830 Illumina NovaSeq 6000 GSM8393297 PAIRED cDNA
## LibrarySource Organism Platform ReleaseDate
## SRR29804826 TRANSCRIPTOMIC Mus musculus ILLUMINA 2024-07-13T00:00:00Z
## SRR29804827 TRANSCRIPTOMIC Mus musculus ILLUMINA 2024-07-13T00:00:00Z
## SRR29804828 TRANSCRIPTOMIC Mus musculus ILLUMINA 2024-07-13T00:00:00Z
## SRR29804829 TRANSCRIPTOMIC Mus musculus ILLUMINA 2024-07-13T00:00:00Z
## SRR29804830 TRANSCRIPTOMIC Mus musculus ILLUMINA 2024-07-13T00:00:00Z
## create_date version Sample.Name source_name SRA.Study
## SRR29804826 2024-07-12T13:39:00Z 1 GSM8393301 cardiomyocytes SRP519546
## SRR29804827 2024-07-12T13:40:00Z 1 GSM8393300 cardiomyocytes SRP519546
## SRR29804828 2024-07-12T13:40:00Z 1 GSM8393299 cardiomyocytes SRP519546
## SRR29804829 2024-07-12T13:40:00Z 1 GSM8393298 cardiomyocytes SRP519546
## SRR29804830 2024-07-12T13:44:00Z 1 GSM8393297 cardiomyocytes SRP519546
## treatment
## SRR29804826 alpha-KG injection
## SRR29804827 alpha-KG injection
## SRR29804828 alpha-KG injection
## SRR29804829 PBS injection
## SRR29804830 PBS injection
# CHANGE AS NEEDED
# In order to make sure each sample is annotated with the correct information,
# we need to define our groups/treatments and set these as factors.
# Donor <- factor(metadata$Donor)
# Arsenic <- factor(metadata$Arsenic)
# To have untreated as reference you need to put levels argument
# The FIRST level listed is treated as the reference group (order matters)
treatment <- factor(metadata$treatment, levels = c("PBS injection", "alpha-KG injection"))
groups <- data.frame(Sample = rownames(metadata), treatment)
# Check it out to make sure everything is in the same order as your samples
cbind(groups, colnames(txi_counts)) # Looks good!
## Sample treatment colnames(txi_counts)
## 1 SRR29804826 alpha-KG injection SRR29804826
## 2 SRR29804827 alpha-KG injection SRR29804827
## 3 SRR29804828 alpha-KG injection SRR29804828
## 4 SRR29804829 PBS injection SRR29804829
## 5 SRR29804830 PBS injection SRR29804830
## 6 SRR29804831 PBS injection SRR29804831
# Let's rename the counts columns to something more informative
groups$ColNames <- paste(groups$treatment,
sep = "_")
colnames(txi_counts) <- groups$ColNames
colnames(txi_counts)
## [1] "alpha-KG injection" "alpha-KG injection" "alpha-KG injection"
## [4] "PBS injection" "PBS injection" "PBS injection"
# Let's take a look at our data!
head(txi_counts)
## alpha-KG injection alpha-KG injection alpha-KG injection
## ENSMUSG00000000001 860.058 752.006 891.041
## ENSMUSG00000000003 0.000 0.000 0.000
## ENSMUSG00000000028 48.001 51.000 36.002
## ENSMUSG00000000037 82.992 72.957 106.194
## ENSMUSG00000000049 1.000 0.000 1.000
## ENSMUSG00000000056 1460.882 1033.427 1168.653
## PBS injection PBS injection PBS injection
## ENSMUSG00000000001 451.121 552.919 489.607
## ENSMUSG00000000003 0.000 0.000 0.000
## ENSMUSG00000000028 57.000 88.001 58.000
## ENSMUSG00000000037 39.012 63.003 58.001
## ENSMUSG00000000049 0.000 1.000 1.000
## ENSMUSG00000000056 868.843 1077.856 855.352
# How many genes and samples do we have?
dim(txi_counts) # of genes and # of samples
## [1] 35380 6
Looking at technical variation of raw data
# Let's start looking at technical variation within our data
# Total counts per sample (library size)
# Adjust the margins to make room for x-axis labels
par(mar = c(8, 4, 4, 2)) # Increase the bottom margin to fit labels
# Create the bar plot
barplot(colSums(txi_counts),
las = 2, # Rotate x-axis labels to be vertical
col = "lightsteelblue2",
cex.names = 0.7) # Adjust the size of the x-axis labels
# Are there any systematic differences in total counts by treatment?
# Treatment
boxplot(colSums(txi_counts) ~ treatment,
col=c("grey","lightsteelblue2"),
ylab="Counts (Sum)",
xlab="Treatment")
# Perform statistical analysis
t.test(colSums(txi_counts) ~ treatment)
##
## Welch Two Sample t-test
##
## data: colSums(txi_counts) by treatment
## t = 4.662, df = 3.7441, p-value = 0.01121
## alternative hypothesis: true difference in means between group PBS injection and group alpha-KG injection is not equal to 0
## 95 percent confidence interval:
## 1148629 4771366
## sample estimates:
## mean in group PBS injection mean in group alpha-KG injection
## 20539817 17579819
Normality of data
# Many stats tests assume a normal distribution of data. Do our raw counts look
# normally distributed?
# Count distribution
# 0 values are an issue, need to exclude low abundance genes
hist(log2(rowSums(txi_counts)), col="grey",
main="Expressed Genes",
xlab="Counts (raw)")
# Let's see if we can fix up that distribution.
# Excluding low abundance genes
MinVals <- apply(txi_counts, 1, min) # Find the minimum value in each row
# How many genes have 0 counts?
sum(MinVals == 0) # of genes. Holy moley, that's a lot that we don't want.
## [1] 20315
Exp <- txi_counts[MinVals > 0, ] # Remove rows (transcripts) with zero counts
# Now that zeros are removed
# The most common way to begin to process your data is through log2 transformation.
# This will give our counts data the normal distribution that our subsequent
# gene expression analysis will expect
# Log2 transformation
ExpLog2 <- log2(Exp)
hist(rowMeans(ExpLog2), xlab = "Counts (log2)", main = "Expressed Genes",
sub = "Removed Zeros",
col="lightsteelblue2")
# That's looking a lot better!
# Now let's take a look at these counts across our samples
# Plot "raw" (our log2 transformed) sample counts
boxplot(ExpLog2, ylab = "log2 counts", main = "Raw RNA-Seq Counts", las = 2)
# With this "uneven" distribution of counts, how do our samples cluster?
# Are there any outliers before normalization?
# Dist calculates distances between rows, use t() to transpose
RawDist <- dist(t(ExpLog2), method = "euclidean")
plot(hclust(RawDist, method = "average"), xlab="Sample")
# Although most like samples are clustering together,
# we can see some samples not clustering by treatment.
# So, we will need to normalize our data!
Simple Normalization (Correction Factor)
# Simple normalization
SampleMedians <- apply(ExpLog2, 2, median) # Find the median value of each column
GrandMedian <- mean(SampleMedians) # Take the average of those
CorrectionFactors <- GrandMedian - SampleMedians # Calculate correction factor to apply to data
CorrectionFactors
## alpha-KG injection alpha-KG injection alpha-KG injection PBS injection
## -0.144586826 0.005127821 -0.152353512 0.196092553
## PBS injection PBS injection
## -0.038626870 0.134346833
# These correction factors will align the medians of all of our samples
# This changes our counts just enough to correct for variability between
# samples/groups without loosing any actual effects of changes in gene expression.
# Loop through each column (sample) and fill in medians adjusted with
# sample correction factor
ExpNorm <- ExpLog2
for(col in colnames(ExpNorm)){
ExpNorm[, col] <- ExpLog2[, col] + CorrectionFactors[col]
}
# Visualize normalized data
# We want to see the medians all in a straight line
# Boxplot of normalized counts
boxplot(ExpNorm, ylab = "log2 counts", main = "Normalized Counts", las = 2, col="lightsteelblue2")
# And do our treatments cluster together better now?
# Cluster dendrogram of normalized data
NormDist <- dist(t(ExpNorm), method = "euclidean")
plot(hclust(NormDist, method = "average"), xlab="Sample")
# Yes, yes they do
# PCA Plot
PCA <- prcomp(t(NormDist))
plot(PCA$x[ , 1], PCA$x[ , 2], pch = 16)
# But that doesn't tell you too much, does it?
# You can color the PCA plot by factor to visualize what specific samples
# may be outliers/driving variation in your data
# By donor
#plot(PCA$x[ , 1], PCA$x[ , 2], pch = 16, col = Donor,
# main = "Colored by Donor")
# You can even add a legend if you'd like
#legend("topleft", legend = unique(Donor), pch = 16, col = unique(Donor))
# By treatment
plot(PCA$x[ , 1], PCA$x[ , 2], pch = 16, cex = 1.5, col = treatment,
main = "PCA",
ylab="PC2",
xlab="PC1")
legend("topright", legend = unique(treatment), pch = 16, cex = 1.5, col = unique(treatment))
# Take home message:
# Always perform an exploratory analysis before your differential gene
# expression analysis!
# This helps you to:
# 1. Make sure your experiment worked
# 2. Find sources of variation within your data
# 3. Find samples that may be outliers in your data
# 4. Visualize your raw data before normalization and your normalized
# data before you perform differential expression analysis or any
# statistical tests
edgeR
# We'll start by creating a DGE object using our counts file
# We will set the "counts =" to our counts dataframe, remove any rows with zero
# counts, and set our gene names are the row names in our counts table.
# You can set the genes argument to an annotation table that contains multiple
# gene identifiers as well (ENSEMBL, ENTREZ, etc.)
DGE <- DGEList(counts = txi_counts, remove.zeros = TRUE, genes = rownames(txi_counts))
## Repeated column names found in count matrix
## Removing 14956 rows with all zero counts
DGE_df <- as.data.frame(DGE)
head(DGE_df)
## genes alpha-KG injection alpha-KG injection alpha-KG injection
## 1 ENSMUSG00000000001 860.058 752.006 891.041
## 2 ENSMUSG00000000028 48.001 51.000 36.002
## 3 ENSMUSG00000000037 82.992 72.957 106.194
## 4 ENSMUSG00000000049 1.000 0.000 1.000
## 5 ENSMUSG00000000056 1460.882 1033.427 1168.653
## 6 ENSMUSG00000000058 1193.490 1378.757 1342.572
## PBS injection PBS injection PBS injection
## 1 451.121 552.919 489.607
## 2 57.000 88.001 58.000
## 3 39.012 63.003 58.001
## 4 0.000 1.000 1.000
## 5 868.843 1077.856 855.352
## 6 1615.895 1665.718 1599.284
# CHANGE AS NEEDED
# Create design matrix
# This will inform our DGE object of which samples will be reference conditions
# for each factor we are interested in
design <- model.matrix(~ 0 + treatment)
rownames(design) <- colnames(DGE)
design
## treatmentPBS injection treatmentalpha-KG injection
## alpha-KG injection 0 1
## alpha-KG injection 0 1
## alpha-KG injection 0 1
## PBS injection 1 0
## PBS injection 1 0
## PBS injection 1 0
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$treatment
## [1] "contr.treatment"
# Sanity check removed low counts
head(DGE_remove_zeros_check, 10)
## genes alpha-KG injection alpha-KG injection alpha-KG injection
## 1 ENSMUSG00000000001 860.058 752.006 891.041
## 2 ENSMUSG00000000028 48.001 51.000 36.002
## 3 ENSMUSG00000000037 82.992 72.957 106.194
## 4 ENSMUSG00000000056 1460.882 1033.427 1168.653
## 5 ENSMUSG00000000058 1193.490 1378.757 1342.572
## 6 ENSMUSG00000000078 2415.239 2425.265 2241.423
## 7 ENSMUSG00000000085 1052.936 936.326 1162.540
## 8 ENSMUSG00000000088 4835.885 4754.421 4501.603
## 9 ENSMUSG00000000093 585.770 582.917 509.874
## 10 ENSMUSG00000000120 40.002 54.000 65.561
## PBS injection PBS injection PBS injection
## 1 451.121 552.919 489.607
## 2 57.000 88.001 58.000
## 3 39.012 63.003 58.001
## 4 868.843 1077.856 855.352
## 5 1615.895 1665.718 1599.284
## 6 1123.035 1168.220 1045.878
## 7 534.430 696.956 555.495
## 8 6093.690 6828.672 6119.166
## 9 470.307 533.065 484.544
## 10 34.000 41.000 34.000
Normalize DGE (TMM)
# Normalize library sizes
# So, instead of using the total library size,
# which is the sum of the reads to all of the genes, we will use TMM normalization.
# Trimmed mean of M-values, where M-values are the log fold change between each
# sample and a reference.
# TMM trims off the most highly variable genes and then calculates a normalization
# factor that is used to adjust the library size
# Normalize for different library sizes (TMM, trimmed mean of M-values)
DGE <- calcNormFactors(DGE) # Calculates normalization factors and adds them to our DGE object
DGE$samples
## group lib.size norm.factors
## 1 1 17352246 1.2064802
## 2 1 17055056 1.1115677
## 3 1 18332156 1.1570329
## 4 1 19606884 0.8455670
## 5 1 21337947 0.9121537
## 6 1 20674619 0.8355687
# Now on to dispersion
# Dispersion means biological coeffient of variation (BCV) squared (BCV^2).
# This is a measure of variability
# Ex: If gene expression typically differs from replicate to replicate by 20%
# its BCV is 0.2, and its dispersion is 0.04.
# Dispersion can also vary based on the biological model you are using.
# For example, human data could have a BCV of 0.4,
# while genetically identical model organisms could be 0.1.
# edgeR estimates dispersion from replicates using the quantile-adjusted
# conditional maximum likelihood method (qCML). The qCML method is designed to
# address composition bias, which occurs when the proportion of reads mapping
# to different genes varies across samples. This bias can arise due to factors
# such as differences in sequencing depths or varying transcript abundances.
# The qCML method uses 2 main types of dispersion:
# Common dispersion calculates a common dispersion value for all genes,
# while the tagwise method calculates gene-specific dispersions.
# Calculate overall dispersion
# estimateGLMCommonDisp calculates within-group and between-group variability.
DGE <- estimateGLMCommonDisp(DGE, design, verbose = TRUE) #Disp = 0.01293 , BCV = 0.1137
## Disp = 0.01293 , BCV = 0.1137
# For this data: Disp = 0.0971 , BCV = 0.3116.
# So on average, the true abundance for each gene can vary up or down by ~31%
# between replicates
# Calculate dispersion trend based on gene abundance
DGE <- estimateGLMTrendedDisp(DGE, design)
# Calculate separate dispersion for each gene
DGE <- estimateGLMTagwiseDisp(DGE, design)
# Visualize dispersion
# counts per million, BCV = variation, each dot is a gene
# more variation in beginning is better, so trend should look like this
plotBCV(DGE)
#Good Results:
#U-Shaped Trend: Ideally, you want to see a U-shaped trend where dispersion is higher at low counts, decreases, and then stabilizes at higher counts. This indicates that the variability in your data is behaving as expected.
#Smooth Trend Line: The trend line should be smooth without large, erratic jumps. This smoothness indicates that the model is fitting well to the data.
#Interpretation of Points:
# Each dot represents a gene.
#Genes with low CPM but high BCV are more variable and less reliable for differential expression analysis.
#Genes with high CPM and low BCV are more consistent and reliable.
# Let's check that our normalization worked by plotting the log2(CPM) PCA
CPM <- cpm(DGE, normalized.lib.sizes = TRUE, log = TRUE)
CPM <- as.data.frame(CPM)
# PCA CPM Log2
pca <- prcomp(t(CPM), scale.=TRUE)
percentVar <- pca$sdev^2 / sum(pca$sdev^2) * 100
pca_df <- as.data.frame(pca$x)
pca_df$group <- DGE$samples$group
ggplot(pca_df, aes(PC1, PC2, color=treatment)) +
geom_point(size=4) +
theme_minimal() +
labs(
title="PCA of RNA-seq Data",
x=paste0("PC1: ", round(percentVar[1], 2), "% variance"),
y=paste0("PC2: ", round(percentVar[2], 2), "% variance")
) +
scale_color_manual(values=c("PBS injection"="black", "alpha-KG injection"="lightsteelblue2"))
# Visualize normalized data
boxplot(CPM, las = 2, ylab = "log2 CPM", main = "Normalized Data", col="lightsteelblue2")
# Now for the main point of this session, differential gene expression!
# We first should consider what specific comparisons we would like to make
# What's our research question and what would I like to know?
# We can make a contrast matrix to specify specific comparisons of interest
# Condition1 - Condition2 will result in the log2FC in Condition1 using
# Condition2 as a reference
# If there are only two factor levels and the reference level doesn't appear
# in our design matrix, we only need to specify the name of the level
# Add as many specific comparisons as you would like
# Create the design matrix with valid column names
colnames(design) <- make.names(colnames(design))
print(colnames(design))
## [1] "treatmentPBS.injection" "treatmentalpha.KG.injection"
contrast.matrix <- makeContrasts(KG_vs_Unt = treatmentalpha.KG.injection - treatmentPBS.injection,
levels = design)
lrt Test
# Sanity check
topTags(lrt)
## Coefficient: -1*treatmentPBS.injection 1*treatmentalpha.KG.injection
## genes logFC logCPM LR PValue
## ENSMUSG00000019997 ENSMUSG00000019997 2.867284 6.632305 726.4103 5.404433e-160
## ENSMUSG00000026202 ENSMUSG00000026202 -2.681592 7.532976 660.6129 1.099103e-145
## ENSMUSG00000030137 ENSMUSG00000030137 -3.871167 5.923674 531.5220 1.317747e-117
## ENSMUSG00000036854 ENSMUSG00000036854 -2.708684 7.930601 503.5062 1.640950e-111
## ENSMUSG00000028766 ENSMUSG00000028766 -3.563794 4.311521 488.5635 2.926750e-108
## ENSMUSG00000027875 ENSMUSG00000027875 3.261380 7.686943 480.7559 1.463041e-106
## ENSMUSG00000043811 ENSMUSG00000043811 -5.273222 3.548905 474.6156 3.172370e-105
## ENSMUSG00000040350 ENSMUSG00000040350 -3.417233 5.427093 464.0827 6.214839e-103
## ENSMUSG00000051855 ENSMUSG00000051855 2.412457 7.919094 458.0054 1.306000e-101
## ENSMUSG00000045776 ENSMUSG00000045776 -2.765881 6.590464 422.6877 6.345084e-94
## FDR
## ENSMUSG00000019997 7.527295e-156
## ENSMUSG00000026202 7.654155e-142
## ENSMUSG00000030137 6.117861e-114
## ENSMUSG00000036854 5.713789e-108
## ENSMUSG00000028766 8.152756e-105
## ENSMUSG00000027875 3.396207e-103
## ENSMUSG00000043811 6.312110e-102
## ENSMUSG00000040350 1.082003e-99
## ENSMUSG00000051855 2.021107e-98
## ENSMUSG00000045776 8.837434e-91
Final Results (From lrt)
# Here are the prelim results! You get a table that has your genes, the logFC,
# logCPM, LR (likelihood ratio test statistics), p-value, and FDR-adjusted p-value
# Now see how many genes are either up- or down-regulated w/ treatment
de <- decideTests(lrt, adjust.method = "fdr")
summary(de)
## -1*treatmentPBS.injection 1*treatmentalpha.KG.injection
## Down 3308
## NotSig 7247
## Up 3373
# Save your results
Results <- as.data.frame(topTags(lrt, n = dim(DGE)[1]))
# Merge with annotations
Annot <- read.csv("AnnotTable.csv",
header = TRUE)
# Merge your added annotation information with your results
Results <- merge(Annot, Results, by.x = "ENSEMBL", by.y = "genes")
# Save your results in a .csv file that you can open with excel
write.csv(Results, file = "log2FC_KG_vs_PBS.csv", row.names = FALSE)
# CHANGE AS NEEDED
# Significance thresholds
# But how many of those are actually significant?
# Save the ones with p-value < 0.05 and log2FC > 1
# You can change these values to be whatever threshold you'd like
# Using the unadj. p-value and FC can produce more reproducible results than
# using FDR alone!
DE <- Results[Results$PValue < 0.05 & abs(Results$logFC) > 0.3, ]
# Save those gene names as "detags"
detags <- DE$ENSEMBL
detags2 <- DE$SYMBOL
# Make a volcano plot
volcano <- plot(-log10(Results$PValue) ~ Results$logFC,
xlab = "log2 fold change", ylab = "-log10 p-value",
main = "DE Genes with Treatment", xlim = c(-10, 10))
points(-log10(Results$PValue[Results$ENSEMBL %in% detags2]) ~
Results$logFC[Results$ENSEMBL %in% detags2],
col = "grey", pch = 16)
abline(v = c(-1, 1), col = "lightsteelblue2")
legend("topleft", legend = "DE genes (FC > 2, p < 0.05)", pch = 16,
col = "grey", bty = "n")
# Want to label the symbols of some genes of interest?
# Let's label the top X genes (highest absolute FC)
Results_arranged <- arrange(Results, PValue)
top4_Results <- Results_arranged[1:20, ]
text(-log10(top4_Results$PValue) ~ top4_Results$logFC,
labels = top4_Results$SYMBOL, cex = 0.9, font = 2)
# Heatmap
# Next, lets define our colors so we can visualize clustering
# You need a character vector to represent each sample's annotation color
# This vector must be in the same order + have the same # columns as your
# heatmap data!
Colors <- as.vector(colnames(CPM))
Colors[grepl("PBS injection", Colors)] <- "Blue"
Colors[grepl("alpha-KG injection", Colors)] <- "Red"
head(Colors)
## [1] "Red" "Red" "Red" "Blue" "Blue" "Blue"
heatmap.2(as.matrix(CPM[rownames(CPM) %in% detags, ]),
scale = "row", trace = "none", labRow = FALSE, srtCol=45,
margins = c(8, 6), ColSideColors = Colors)
# ColSideColors is what adds the annotation
Genes of Interest
# Assume lrt_results is your LRT results object
# Convert lrt results to data frame
lrt_df <- as.data.frame(Results)
# Example genes of interest (replace with your actual genes)
genes_of_interest <- c("Erbb4", "Erbb2", "Nrg1", "Acvrl1", "Ccn2")
# Filter results for these genes
filtered_lrt_df <- lrt_df[lrt_df$SYMBOL %in% genes_of_interest, ]
# Plot logFC values with bars for reference and treatment
ggplot(filtered_lrt_df, aes(x = SYMBOL, y = logFC)) +
geom_bar(stat = "identity", position = "dodge", fill = "lightsteelblue2") + # Dodge position to place bars side by side
labs(x = "Gene", y = "Log Fold Change (logFC)", title = "LogFC for Genes of Interest") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
GSEA and GO
# This module needs an object from the edgeR module
# Import gene expression table - from the EdgeR module
edgeRdata <- read.csv("log2FC_KG_vs_PBS.csv")
head(edgeRdata)
## ENSEMBL ENTREZID SYMBOL GENENAME
## 1 ENSMUSG00000000001 14679 Gnai3 G protein subunit alpha i3
## 2 ENSMUSG00000000028 12544 Cdc45 cell division cycle 45
## 3 ENSMUSG00000000037 107815 Scml2 Scm polycomb group protein like 2
## 4 ENSMUSG00000000056 67608 Narf nuclear prelamin A recognition factor
## 5 ENSMUSG00000000058 12390 Cav2 caveolin 2
## 6 ENSMUSG00000000078 23849 Klf6 Kruppel-like transcription factor 6
## UNIPROT logFC logCPM LR PValue FDR
## 1 A2AE36 0.5478569 5.113070 24.156352 8.882269e-07 5.262112e-06
## 2 F8WJ72 -0.7665309 1.622421 7.453704 6.330603e-03 1.572546e-02
## 3 B1AVB5 0.5159668 1.912651 4.101894 4.283520e-02 8.243867e-02
## 4 Q8BVW9 0.1895550 5.813677 2.580349 1.081978e-01 1.797232e-01
## 5 Q924U3 -0.5140279 6.288215 22.575200 2.020763e-06 1.124908e-05
## 6 Q3UGJ1 0.8901408 6.487324 75.306899 4.029493e-18 1.081364e-16
# Identify up-regulated differentially expressed (DE) genes (FDR < 0.05)
SigGenes <- edgeRdata[which(edgeRdata$FDR<0.05 & edgeRdata$logFC>0),"ENSEMBL"]
length(SigGenes) #1032 DE genes
## [1] 3373
# down-regulated
SigGenesdown <- edgeRdata[which(edgeRdata$FDR<0.05 & edgeRdata$logFC<0),"ENSEMBL"]
length(SigGenesdown) #1032 DE genes
## [1] 3308
# all genes detected
ExpGenes<-edgeRdata$ENSEMBL
# the Venn diagram helps you to visualize the number and construct the contingency table
venn(list("Up regulated"=SigGenes,
"Down Regulated"=SigGenesdown,
"Expressed"=ExpGenes))
Goana
# Top significant results
head(GoanaResults, 10)
## Term Ont N Up Down P.Up
## GO:0006066 alcohol metabolic process BP 253 69 69 1.787561e-01
## GO:0006629 lipid metabolic process BP 1004 254 267 3.116885e-01
## GO:0006807 nitrogen compound metabolic process BP 7362 1857 1694 3.741783e-02
## GO:0006810 transport BP 3304 860 810 1.604953e-02
## GO:0006869 lipid transport BP 303 87 67 5.559856e-02
## GO:0006897 endocytosis BP 577 166 126 1.103266e-02
## GO:0006898 receptor-mediated endocytosis BP 219 63 52 8.846764e-02
## GO:0008104 protein localization BP 2219 618 463 7.084678e-05
## GO:0008150 biological_process BP 13081 3234 3124 7.050006e-02
## GO:0008152 metabolic process BP 8499 2124 2016 9.584681e-02
## P.Down FDR.Up FDR.Down
## GO:0006066 0.11568157 0.887118342 1.0000000
## GO:0006629 0.02083256 1.000000000 0.6121146
## GO:0006807 0.99554438 0.420268954 1.0000000
## GO:0006810 0.16954802 0.244849219 1.0000000
## GO:0006869 0.78682501 0.522433779 1.0000000
## GO:0006897 0.89113313 0.197339587 1.0000000
## GO:0006898 0.54591560 0.662839983 1.0000000
## GO:0008104 0.99989873 0.003871295 1.0000000
## GO:0008150 0.53434332 0.576604316 1.0000000
## GO:0008152 0.72515036 0.704258068 1.0000000
# Less than 200 genes in enrichment result!
head(Goana_Results_top10_BP_less200,10)
## Term Ont N Up Down P.Up
## GO:0060038 cardiac muscle cell proliferation BP 55 35 5 9.683459e-10
## GO:0014855 striated muscle cell proliferation BP 75 41 7 2.252438e-08
## GO:0021915 neural tube development BP 162 72 24 2.399149e-08
## GO:0003205 cardiac chamber development BP 171 75 40 2.464652e-08
## GO:0003231 cardiac ventricle development BP 131 61 33 3.461409e-08
## GO:0034332 adherens junction organization BP 46 29 8 3.640241e-08
## GO:0003281 ventricular septum development BP 77 41 16 6.094599e-08
## GO:0055123 digestive system development BP 96 48 13 6.411563e-08
## GO:0001841 neural tube formation BP 114 54 20 1.025749e-07
## GO:0045216 cell-cell junction organization BP 152 67 29 1.071997e-07
## P.Down FDR.Up FDR.Down
## GO:0060038 0.9987252 2.710477e-07 1
## GO:0014855 0.9996755 4.095395e-06 1
## GO:0021915 0.9984371 4.288832e-06 1
## GO:0003205 0.5899627 4.369212e-06 1
## GO:0003231 0.3951325 5.664196e-06 1
## GO:0034332 0.8898378 5.866580e-06 1
## GO:0003281 0.7779253 9.132443e-06 1
## GO:0055123 0.9959800 9.537984e-06 1
## GO:0001841 0.9597549 1.416932e-05 1
## GO:0045216 0.9352366 1.461833e-05 1
GSEA
# GSEA Results
head(GSEA_GO_df,10)
## ID Description
## GO:0006119 GO:0006119 oxidative phosphorylation
## GO:0009060 GO:0009060 aerobic respiration
## GO:0045333 GO:0045333 cellular respiration
## GO:0022900 GO:0022900 electron transport chain
## GO:0009145 GO:0009145 purine nucleoside triphosphate biosynthetic process
## GO:0009206 GO:0009206 purine ribonucleoside triphosphate biosynthetic process
## GO:0032543 GO:0032543 mitochondrial translation
## GO:0009142 GO:0009142 nucleoside triphosphate biosynthetic process
## GO:0009201 GO:0009201 ribonucleoside triphosphate biosynthetic process
## GO:0140053 GO:0140053 mitochondrial gene expression
## setSize enrichmentScore NES pvalue p.adjust qvalue rank
## GO:0006119 136 -0.7546599 -3.075602 1e-10 3.78e-09 2.8e-09 2172
## GO:0009060 180 -0.7055015 -2.985293 1e-10 3.78e-09 2.8e-09 2172
## GO:0045333 228 -0.6647674 -2.908973 1e-10 3.78e-09 2.8e-09 2122
## GO:0022900 105 -0.7325576 -2.839381 1e-10 3.78e-09 2.8e-09 2172
## GO:0009145 107 -0.6866816 -2.670389 1e-10 3.78e-09 2.8e-09 2329
## GO:0009206 106 -0.6852837 -2.663010 1e-10 3.78e-09 2.8e-09 2329
## GO:0032543 130 -0.6536067 -2.651364 1e-10 3.78e-09 2.8e-09 2692
## GO:0009142 117 -0.6664917 -2.642250 1e-10 3.78e-09 2.8e-09 2329
## GO:0009201 111 -0.6696953 -2.626507 1e-10 3.78e-09 2.8e-09 2329
## GO:0140053 167 -0.6110423 -2.551284 1e-10 3.78e-09 2.8e-09 3132
## leading_edge
## GO:0006119 tags=76%, list=16%, signal=64%
## GO:0009060 tags=67%, list=16%, signal=57%
## GO:0045333 tags=60%, list=16%, signal=52%
## GO:0022900 tags=72%, list=16%, signal=61%
## GO:0009145 tags=67%, list=17%, signal=56%
## GO:0009206 tags=67%, list=17%, signal=56%
## GO:0032543 tags=69%, list=20%, signal=56%
## GO:0009142 tags=64%, list=17%, signal=54%
## GO:0009201 tags=66%, list=17%, signal=55%
## GO:0140053 tags=62%, list=23%, signal=49%
## core_enrichment
## GO:0006119 12866/17995/13382/227197/17708/12064/226646/11950/66072/12858/12861/67892/14004/107197/67003/66416/66576/214917/66152/17721/12857/105675/28080/67273/17706/68073/66218/68550/104130/17720/66108/66148/66052/67130/67914/17722/102060/66445/17709/66142/17705/66046/71679/57320/66043/68197/72900/66383/68194/68349/17719/17710/22273/17991/17716/67530/11957/17993/12859/17711/230075/66414/56312/12867/54405/67680/67184/57423/67264/12862/84682/12850/68375/102631912/66091/78330/17718/103172/66694/11958/70316/407785/66594/75406/69875/22272/68198/67267/67126/595136/68202/66377/215951/66916/69702/17717/225887/68342/13063/66495/12869/58805/12865
## GO:0009060 12866/17995/13382/227197/17708/12064/226646/11950/11429/66072/12858/12861/67892/14004/107197/67003/14194/66416/66576/78920/20916/170718/214917/56451/66152/70383/17721/12857/105675/12974/14533/17448/28080/67273/17706/68073/66218/68550/104130/17720/66108/66148/66052/17449/67130/67914/17722/102060/68263/67834/66445/17709/66142/17705/66046/71679/57320/66043/68197/235339/72900/66383/68194/68349/17719/17710/22273/17991/17716/67530/11957/17993/18293/12859/17711/230075/66414/56312/12867/54405/66128/63873/67680/67184/57423/67264/12862/84682/12850/68375/102631912/66091/78330/17718/103172/66694/11958/70316/407785/66594/75406/69875/22272/68198/67267/67126/595136/68202/66377/215951/66916/69702/17717/225887/68342/13063/66495/12869/58805/239017/12865
## GO:0045333 17995/13382/227197/17708/12064/226646/11950/11429/66072/12858/72170/12861/67892/246730/14004/107197/67003/14194/66416/13603/66576/78920/20916/170718/214917/56451/66152/70383/17721/12857/105675/12974/14533/17448/28080/67273/17706/68073/66218/68550/104130/100126824/17720/66108/66148/18124/66052/17449/67130/67914/17722/102060/68263/67834/66445/17709/20656/66142/17705/66046/71679/57320/66043/68197/235339/72900/66383/68194/68349/17719/17710/22273/17991/17716/67530/11957/17993/18293/12859/17711/93757/230075/110826/66414/56312/12867/54405/66128/63873/67680/67184/57423/67264/12862/84682/12850/68375/23960/102631912/66091/56857/78330/17718/103172/67885/66694/52637/11958/70316/407785/66594/75406/69875/22272/68198/67267/67126/67876/595136/75530/68202/66377/215951/68267/66916/69702/17717/225887/68342/13063/66495/12869/68002/14555/58805/239017/12865
## GO:0022900 12866/17995/13382/227197/17708/12064/226646/66072/12858/12861/67892/68165/14004/107197/67003/66416/66576/66152/17721/12857/67273/66218/100126824/17720/66148/66052/67914/17722/66445/17709/20656/66142/57320/68197/14148/72900/66383/68349/17719/17710/22273/17716/67530/17993/12859/17711/93757/230075/110826/66414/12867/67680/67264/12862/84682/12850/68375/56857/78330/17718/66694/407785/66594/75406/22272/68202/215951/68267/69702/17717/225887/13063/66495/12869/14555/12865
## GO:0009145 66114/18103/11674/18950/17995/227197/226646/11950/50873/107197/66416/214917/17721/69802/11534/28080/67273/17706/68073/66218/104130/17720/66108/27425/66052/67130/17722/18102/17705/66046/71679/66043/68197/72900/68194/68349/17719/17991/17716/11957/17993/230075/66414/228033/54405/67680/67184/57423/67264/68375/102631912/66091/78330/17718/11958/70316/407785/79059/75406/11951/69875/68198/67126/595136/68202/66377/66916/17717/225887/68342/66495/52815
## GO:0009206 66114/18103/11674/18950/17995/227197/226646/11950/50873/107197/66416/214917/17721/69802/28080/67273/17706/68073/66218/104130/17720/66108/27425/66052/67130/17722/18102/17705/66046/71679/66043/68197/72900/68194/68349/17719/17991/17716/11957/17993/230075/66414/228033/54405/67680/67184/57423/67264/68375/102631912/66091/78330/17718/11958/70316/407785/79059/75406/11951/69875/68198/67126/595136/68202/66377/66916/17717/225887/68342/66495/52815
## GO:0032543 72181/94067/118451/94061/66845/68572/77721/229487/107732/94064/69527/226539/101122/320806/64656/66258/218506/67994/56280/67840/102436/66973/68735/57312/75619/66407/18120/121022/52856/66223/60441/66416/28028/69163/24030/94063/66419/232536/71984/70120/69956/94062/76784/68537/14548/12261/107734/64655/216767/64660/68611/27398/94254/94066/27397/66047/66077/353242/50529/66242/102060/72416/74600/70207/66292/66448/107733/67036/224805/68499/66230/94065/64658/68463/66121/67270/52469/233870/66163/56284/67681/66493/67308/56282/66399/64659/67267/79044/76563/67212
## GO:0009142 66114/18103/11674/18950/69719/17995/227197/80914/226646/11950/50873/107197/66416/214917/17721/69802/11534/28080/67273/17706/68073/66218/104130/17720/66108/27425/66052/67130/17722/18102/11636/17705/66046/71679/66043/68197/72900/68194/68349/17719/17991/17716/11957/17993/230075/66414/228033/54405/67680/67184/57423/67264/68375/102631912/66091/78330/17718/11958/70316/407785/79059/75406/11951/69875/68198/67126/595136/68202/66377/66916/17717/225887/68342/66495/52815
## GO:0009201 66114/18103/11674/18950/69719/17995/227197/80914/226646/11950/50873/107197/66416/214917/17721/69802/28080/67273/17706/68073/66218/104130/17720/66108/27425/66052/67130/17722/18102/17705/66046/71679/66043/68197/72900/68194/68349/17719/17991/17716/11957/17993/230075/66414/228033/54405/67680/67184/57423/67264/68375/102631912/66091/78330/17718/11958/70316/407785/79059/75406/11951/69875/68198/67126/595136/68202/66377/66916/17717/225887/68342/66495/52815
## GO:0140053 320720/72181/94067/118451/94061/66845/68572/224481/64384/77721/229487/107732/94064/21780/69527/226539/101122/320806/64656/66258/218506/67994/56280/67840/71701/102436/66973/68735/57312/75619/66407/216151/18120/74238/121022/52856/66223/67851/60441/66416/28028/69163/24030/94063/66419/232536/71984/70120/69956/94062/76784/68537/14548/12261/107734/64655/216767/64660/68611/27398/94254/94066/27397/66047/68550/66077/353242/50529/66242/66587/102060/338359/72416/74600/70207/66292/66448/107733/67036/170826/224805/68499/66230/94065/64658/68463/66121/67270/52469/233870/21379/66163/56284/67681/103172/66493/67308/56282/66399/64659/67267/79044/76563/67212
# Assuming GSEA_GO is your GSEA result object
# Filter for the top 10 results by NES (normalized enrichment score) or p-value
top10_results <- GSEA_GO@result %>%
dplyr::arrange(desc(NES)) %>%
dplyr::slice(1:10)
top10_results
## ID Description setSize
## GO:0030177 GO:0030177 positive regulation of Wnt signaling pathway 118
## GO:0003206 GO:0003206 cardiac chamber morphogenesis 126
## GO:0000070 GO:0000070 mitotic sister chromatid segregation 172
## GO:0010466 GO:0010466 negative regulation of peptidase activity 118
## GO:0140014 GO:0140014 mitotic nuclear division 240
## GO:0030111 GO:0030111 regulation of Wnt signaling pathway 267
## GO:0060070 GO:0060070 canonical Wnt signaling pathway 250
## GO:0198738 GO:0198738 cell-cell signaling by wnt 367
## GO:0016055 GO:0016055 Wnt signaling pathway 366
## GO:0007292 GO:0007292 female gamete generation 116
## enrichmentScore NES pvalue p.adjust qvalue rank
## GO:0030177 0.5257707 2.100636 5.784822e-08 1.977225e-06 1.464611e-06 1734
## GO:0003206 0.4770743 1.923782 3.385255e-06 8.724725e-05 6.462759e-05 1886
## GO:0000070 0.4475769 1.886465 3.089279e-06 8.341054e-05 6.178559e-05 4351
## GO:0010466 0.4646511 1.856442 2.009894e-05 4.469058e-04 3.310414e-04 600
## GO:0140014 0.4205424 1.840917 4.050287e-07 1.208691e-05 8.953266e-06 4357
## GO:0030111 0.4123730 1.826167 2.654904e-07 8.136923e-06 6.027351e-06 3132
## GO:0060070 0.4092492 1.811428 1.234530e-06 3.499892e-05 2.592513e-05 3660
## GO:0198738 0.3871508 1.802579 5.928189e-08 1.977225e-06 1.464611e-06 3660
## GO:0016055 0.3864180 1.800150 1.905235e-07 6.001491e-06 4.445549e-06 3660
## GO:0007292 0.4475454 1.784904 1.897903e-04 3.416225e-03 2.530537e-03 3738
## leading_edge
## GO:0030177 tags=31%, list=13%, signal=27%
## GO:0003206 tags=33%, list=14%, signal=29%
## GO:0000070 tags=53%, list=32%, signal=37%
## GO:0010466 tags=18%, list=4%, signal=17%
## GO:0140014 tags=51%, list=32%, signal=35%
## GO:0030111 tags=35%, list=23%, signal=27%
## GO:0060070 tags=41%, list=27%, signal=31%
## GO:0198738 tags=40%, list=27%, signal=30%
## GO:0016055 tags=40%, list=27%, signal=30%
## GO:0007292 tags=48%, list=27%, signal=35%
## core_enrichment
## GO:0030177 21752/330938/12560/20677/11538/20356/22402/26564/20377/58198/14296/20319/71461/14734/13607/14183/58799/20779/12316/70122/218581/64297/11865/20729/21951/140577/21833/78560/59036/232560/66241/19664/107515/244667/27373/72114
## GO:0003206 21952/140781/12814/18741/12162/268902/20564/20677/215798/20319/67451/15214/14465/20563/211323/20666/19876/14183/21814/13857/17684/22762/116701/11504/14613/57265/15111/21388/57246/13617/109620/20181/17300/19664/319757/380718/21813/22411/63958/15251/21808/93840
## GO:0000070 56742/51944/229841/12235/381293/193385/12236/268697/230936/73804/72119/74393/26934/70099/107995/70218/100502766/11799/76464/233406/108907/76789/16551/71819/20843/22137/18817/208718/215387/66371/19357/228421/66468/20479/76707/68612/52276/18519/110033/73420/54141/101706/60411/18005/66977/13006/20842/74201/209334/14211/218914/26554/56317/219114/11920/106042/235661/224023/67052/16580/67064/102920/52696/213582/93759/12534/71175/70385/54392/16563/99412/76044/108989/101994/16319/24061/68014/72124/67141/70799/20877/78798/59008/66442/69654/105513/234852/16569/14841/67070/67849/68501
## GO:0010466 20715/20716/100034684/20861/18054/20862/213945/17002/17395/20708/278507/215001/11625/20319/110595/235505/14734/94227/21859/232345/20779
## GO:0140014 21802/12162/16002/56742/12223/51944/229841/16001/12235/381293/13860/14183/193385/12236/268697/230936/73804/434130/72119/17345/74393/26934/70099/107995/70218/12576/16337/100502766/11799/76464/233406/108907/76789/66515/16551/71819/20843/22137/18817/208718/215387/218294/66371/83946/19357/228421/66468/20479/76707/68612/52276/16000/18519/110033/73420/54141/101706/60411/18005/66977/13006/108912/20842/74201/209334/14211/71943/218914/17295/26554/56317/219114/11920/106042/235661/20878/224023/67052/71678/16580/67064/102920/54369/52696/213582/93759/12534/71782/71175/70385/54392/16563/99412/76044/108989/101994/16319/77744/24061/68014/50850/72124/67141/70799/20877/78309/68053/18596/78267/78798/59008/67655/66442/69654/105513/234852/229776/16569/14841/67070/67849/68501/105837
## GO:0030111 21752/330938/230899/12560/50781/20677/11538/20356/22402/26564/56484/20377/58198/14296/20319/71461/14734/13607/74499/14183/58799/20779/68588/66277/12316/70122/19277/76365/14783/17240/219134/218581/64297/11865/330096/12558/83675/432940/20729/29806/21951/140577/21833/78560/59036/24117/228357/232560/66241/19664/18647/107515/244667/232334/380718/27373/72114/71890/22411/213326/21885/74775/93840/14634/494504/17128/81004/229473/207181/72345/621976/56274/74493/214897/13205/70425/106042/80707/21983/21416/240725/16400/56458/260305/14682/12005/68031/12388/407821/17242/65114/67974/18128
## GO:0060070 330938/230899/12560/50781/20677/11538/20356/26564/56484/20377/14296/20319/71461/13618/14734/14367/13607/74499/21414/14183/20779/68588/12316/70122/19277/76365/14365/64297/57265/11865/330096/12558/83675/432940/20969/29806/21951/140577/21833/13617/78560/59036/228357/232560/66241/19664/18647/107515/244667/20437/380718/27373/72114/71890/22411/21415/213326/21885/74775/20479/14634/17128/81004/72135/229473/207181/72345/621976/56274/74493/214897/13205/70425/106042/21983/21416/56458/260305/14682/12005/16480/68031/12388/407821/17242/65114/67974/18128/107351/16973/73668/104318/74996/22419/80288/18073/74025/19211/12387/58231/14370/76441/13649
## GO:0198738 21752/18741/330938/230899/12560/50781/20677/11538/20356/22402/26564/17869/56484/20377/58198/14296/20319/71461/13618/14734/14367/13607/74499/21414/14183/67488/73389/58799/20779/68588/66277/12316/70122/19277/76365/14783/17240/219134/218581/14365/64297/57265/11865/330096/12558/83675/432940/20969/20729/29806/21951/140577/21833/13617/78560/59036/24117/75723/13194/228357/232560/66241/108079/103583/19664/18647/107515/26563/244667/232334/20437/380718/27373/72114/71890/22411/21415/56805/213326/21885/74775/93840/20479/56332/14634/268980/494504/17128/81004/72135/53627/229473/207181/72345/621976/56274/16825/74493/214897/13205/26554/70425/106042/80707/14773/21983/21416/240725/16400/56458/260305/14682/12005/208846/16480/68031/12388/240756/407821/17242/212398/65114/67974/18128/107351/16973/73668/104318/74996/13000/104831/218952/22419/80288/76281/18073/74025/24088/76246/19211/12387/58231/12614/14370/76441/13649
## GO:0016055 21752/18741/330938/230899/12560/50781/20677/11538/20356/22402/26564/17869/56484/20377/58198/14296/20319/71461/13618/14734/14367/13607/74499/21414/14183/67488/73389/58799/20779/68588/66277/12316/70122/19277/76365/14783/17240/219134/218581/14365/64297/57265/11865/330096/12558/83675/432940/20969/20729/29806/21951/140577/21833/13617/78560/59036/24117/75723/13194/228357/232560/66241/108079/103583/19664/18647/107515/26563/244667/232334/20437/380718/27373/72114/71890/22411/21415/56805/213326/21885/74775/93840/20479/56332/14634/268980/494504/17128/81004/72135/53627/229473/207181/72345/621976/56274/16825/74493/214897/13205/26554/70425/106042/80707/14773/21983/21416/240725/16400/56458/260305/14682/12005/208846/16480/68031/12388/240756/407821/17242/212398/65114/67974/18128/107351/16973/73668/104318/74996/13000/218952/22419/80288/76281/18073/74025/24088/76246/19211/12387/58231/12614/14370/76441/13649
## GO:0007292 19242/19229/210510/56484/12442/229841/20779/268697/12316/230594/381404/67121/71567/11504/54004/20729/68767/15081/13617/230126/236266/13194/68549/21973/22137/18817/215387/223593/72634/13121/16000/110147/54611/21930/16324/242202/15212/11920/20112/20878/67052/54204/22187/22697/17242/100121/93759/26920/328365/11651/19211/12387/268903/56739/16323/213541
# Create a dot plot of the top 10 GSEA results
ggplot(top10_results, aes(x = reorder(Description, NES), y = NES, color = p.adjust)) +
geom_point(size = 4) +
coord_flip() +
labs(x = "GO Terms", y = "Normalized Enrichment Score (NES)", title = "Top 10 GSEA Results") +
scale_color_gradient(low = "blue", high = "red") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
# GSEA plot
gseaplot2(GSEA_GO,geneSetID = "GO:0016055", pvalue_table = TRUE)
gseaplot2(GSEA_GO,geneSetID = "GO:0003007", pvalue_table = TRUE)
gseaplot2(GSEA_GO,geneSetID = "GO:0019827", pvalue_table = TRUE)
# Bar plot for the top 10 GO terms
ggplot(top10_results, aes(x = reorder(Description, NES), y = NES, fill = p.adjust)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "GO Terms", y = "Normalized Enrichment Score (NES)", title = "Top 10 GSEA Results (Bar Plot)") +
scale_fill_gradient(low = "blue", high = "red") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
# BIOLOGICAL PROCESS ONLY
# website of code https://rdrr.io/bioc/rWikiPathways/f/vignettes/Pathway-Analysis.Rmd
# Up-regulated only BP
egobp_up <- clusterProfiler::enrichGO(
gene = up.genes.entrez[[2]],
universe = bkgd.genes.entrez[[2]],
OrgDb = org.Mm.eg.db,
ont = "BP",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05,
readable = TRUE)
# Make results a data frame and save to CSV
BP_enrichment_results_upreg <- as.data.frame(egobp_up)
head(BP_enrichment_results_upreg,10)
## ID Description GeneRatio BgRatio
## GO:0048598 GO:0048598 embryonic morphogenesis 187/3241 488/13113
## GO:0060038 GO:0060038 cardiac muscle cell proliferation 35/3241 55/13113
## GO:0003007 GO:0003007 heart morphogenesis 102/3241 243/13113
## GO:0198738 GO:0198738 cell-cell signaling by wnt 141/3241 367/13113
## GO:0016055 GO:0016055 Wnt signaling pathway 140/3241 366/13113
## GO:0048568 GO:0048568 embryonic organ development 144/3241 379/13113
## GO:0060070 GO:0060070 canonical Wnt signaling pathway 102/3241 250/13113
## GO:0014855 GO:0014855 striated muscle cell proliferation 41/3241 75/13113
## GO:0021915 GO:0021915 neural tube development 72/3241 162/13113
## GO:0003205 GO:0003205 cardiac chamber development 75/3241 171/13113
## pvalue p.adjust qvalue
## GO:0048598 8.289006e-12 4.692407e-08 4.004899e-08
## GO:0060038 1.090223e-09 3.085875e-06 2.633748e-06
## GO:0003007 1.963513e-09 3.705149e-06 3.162289e-06
## GO:0198738 2.641642e-09 3.738585e-06 3.190826e-06
## GO:0016055 4.202838e-09 4.005897e-06 3.418973e-06
## GO:0048568 4.245783e-09 4.005897e-06 3.418973e-06
## GO:0060070 1.200302e-08 9.707013e-06 8.284790e-06
## GO:0014855 2.552632e-08 1.674226e-05 1.428927e-05
## GO:0021915 2.866443e-08 1.674226e-05 1.428927e-05
## GO:0003205 2.957474e-08 1.674226e-05 1.428927e-05
## geneID
## GO:0048598 Gna12/Itga5/Tcf7/Mmp14/Sec24b/Sp1/Smo/Prkra/Cul3/Prkacb/Trim28/Gata6/Ctnnb1/Arid1a/Fzd3/Hipk1/Bmp7/Kcnq1/Bcr/Prox1/Lmbr1/Sulf1/Ift52/Tbx5/Kif3a/Pcgf2/Apaf1/Frs2/Cobl/Wdpcp/Rock2/Prkar1a/Id2/Llgl2/Ntn1/Hif1a/Zfp36l1/Gli3/Ror2/Ptch1/Bmpr1a/Gata4/Nipbl/Gdnf/Stk3/Myc/Scrib/Ptk7/Foxp4/Ift140/Axin1/Mib1/Brd2/Smad4/Smad2/Fbn2/Gnaq/Kif20b/Lrp5/Tcf7l2/Ldb1/Sufu/Bbs4/Ric8a/Col5a2/Fn1/Epb41l5/Rnf2/Vangl2/Pou2f1/Itga8/Tsc1/Acvr1/Col5a1/Brd3/Notch1/Nckap1/Itgav/Sp3/Traf6/Ext2/Lrp4/Rtf1/Bcl2l11/Vcam1/Col11a1/Sfrp2/Pitx2/Mmp16/Lmo4/Stil/Syf2/Ski/Dnajb6/Pdgfra/Rest/Glmn/Shroom3/Kdm2b/Exoc4/Tead2/Fgfr2/Dlc1/Dusp4/Sfrp1/Ednra/Casp3/Sall1/Smad1/Irx3/Irx5/Tbx20/Col12a1/Tbx18/Tgfbr2/Mdfi/Ift57/Phldb2/Rpgrip1l/Specc1l/Mks1/Scx/Pbx2/Fras1/Hectd1/Ripor2/Prickle1/Plxnb2/Aff3/Kdm6a/Wdr19/Rara/Cdon/Tbc1d32/Hand2/Tctn1/Kif16b/Rbpj/Tgfb2/Cc2d2a/Cplane1/Cited2/Ift88/Hs2st1/Ofd1/Mthfd1l/Erf/Dact1/Setd2/Fzd5/Megf8/Tcf21/Tshz1/Dync2h1/C2cd3/Tenm4/Igf2/Fzd2/Foxc1/Gja1/Acvr2a/Pbx1/Ndst1/Cthrc1/Tead1/Gpc3/Gja5/Ece1/T2/Arhgap35/Map2k5/Rnf207/Deaf1/Mfap2/Hipk2/Tgif2/Mapk3/Mapk1/Sox11/Clasp1/Bpnt2/Zbtb16/Bmpr2/Zfp568/Spint2/Pdzd7/Ttbk2
## GO:0060038 Gata6/Ctnnb1/Pten/Rxra/Tbx5/Hey2/Cdk1/Prkar1a/Myh10/Bmpr1a/Gata4/Zfpm2/Cxadr/Pim1/Notch1/Tgfbr3/Foxp1/Vgll4/Fgfr2/Smad1/Tbx20/Tgfbr2/Arid2/Rbpj/Tgfb2/Ncam1/Rxrb/Cited2/Ccnb1/Tenm4/Foxc1/Gja1/Erbb4/Nrg1/Mapk1
## GO:0003007 Sec24b/Smo/Insr/Gata6/Crkl/Epor/Tek/Ctnnb1/Zmiz1/Fgfrl1/Bmp7/Mrtfb/Prox1/Rxra/Ift52/Tbx5/Kif3a/Hey2/Rtn4/Id2/Hif1a/Ptch1/Ptcd2/Cert1/Bmpr1a/Gata4/Nipbl/Zfpm2/Angpt1/Robo1/Adamts1/Pim1/Mib1/Smad4/Smad2/Ankrd1/Pdcd4/Sufu/Tnni1/Vangl2/Acvr1/Col5a1/Notch1/Atf2/Col11a1/Sfrp2/Pitx2/Stil/Ube4b/Tgfbr3/Ephb4/Abcc9/Parva/Tead2/Fgfr2/Dlc1/Slit2/Ednra/Tbx20/Nedd4/Tgfbr2/Ift57/Arid2/Mks1/Prickle1/Spry1/Kdm6a/Hand2/Adgrg6/Rbpj/Tgfb2/Synpo2l/Cited2/Ift88/Thbs1/Ppp1r13l/Pkp2/Pbrm1/Asxl1/Megf8/S1pr1/C2cd3/Flrt2/Zfpm1/Fzd2/Foxc1/Gja1/Robo2/Jun/Myh7/Slc8a1/Dsp/Tead1/Slit3/Gja5/Nsd2/Rnf207/Nrg1/Sox11/Bmpr2/Sox4/Xirp1
## GO:0198738 Tcf7/Ddx3x/Sost/Itga3/Jup/Hbp1/Grk5/Cul3/Wwox/Ccn4/Crbn/Nid1/Ctnnb1/Fzd3/Tle1/Gprc5b/Vps4b/Amotl1/Pten/Bicc1/Sulf1/Grb10/Wif1/Fbxw11/Kremen1/Sdc1/Ppm1a/Gli3/Spin1/Ror2/Tert/Ednrb/Sema5a/Stk3/Myc/Csnk1e/Gsk3b/Lmbr1l/Calcoco1/Ptk7/Strn/Axin1/Ccny/Cdh2/Smad4/Csnk1a1/Gnaq/Ddb1/Tnks2/Lrp5/Tcf7l2/Csnk1d/Ldb1/Limd1/Tmem9/Klhl12/Vangl2/Notch1/Lrp4/Rtf1/Tbl1xr1/Src/Sfrp2/Pitx2/Mllt3/Prkaa2/Cdk14/Ctnnbip1/Mad2l2/Klf15/Ptpro/Caprin2/Vgll4/Dkk3/Fgfr2/Zranb1/Adgra2/Tnks/Sfrp1/Ednra/Sall1/Vps35/Amfr/Dixdc1/Csnk1g1/Tbx18/Amotl2/Ryk/Kank1/Mdfi/Tmem131l/Aspm/Ctnnd1/Mks1/Daam1/Pygo1/Ror1/Ptpn23/Prickle1/Siah1a/Fermt2/Rnf146/Rbpj/Rbms3/Tax1bp3/Ankrd6/Daam2/Nppa/Zbed3/Disc1/Ppp2r3a/Tmem64/Foxo1/Dact1/Fzd5/Otulin/Csnk2a2/Zbtb33/Foxo3/Fzd4/Shisa3/Lgr4/Fzd2/Amer1/Adnp/Chd8/Cthrc1/Bmal1/Gpc3/Tcf7l1/Thra/Usp47/Eda/Cdh3/Bcl9l/Frat1/Scyl2/Apcdd1/Csnk1g3/Sox4/Tmem170b
## GO:0016055 Tcf7/Ddx3x/Sost/Itga3/Jup/Hbp1/Grk5/Cul3/Wwox/Ccn4/Crbn/Nid1/Ctnnb1/Fzd3/Tle1/Gprc5b/Vps4b/Amotl1/Pten/Bicc1/Sulf1/Grb10/Wif1/Fbxw11/Kremen1/Sdc1/Ppm1a/Gli3/Spin1/Ror2/Tert/Ednrb/Sema5a/Stk3/Myc/Csnk1e/Gsk3b/Lmbr1l/Calcoco1/Ptk7/Strn/Axin1/Ccny/Cdh2/Smad4/Csnk1a1/Gnaq/Ddb1/Tnks2/Lrp5/Tcf7l2/Csnk1d/Ldb1/Limd1/Tmem9/Klhl12/Vangl2/Notch1/Lrp4/Rtf1/Tbl1xr1/Src/Sfrp2/Pitx2/Mllt3/Prkaa2/Cdk14/Ctnnbip1/Mad2l2/Klf15/Ptpro/Caprin2/Vgll4/Dkk3/Fgfr2/Zranb1/Adgra2/Tnks/Sfrp1/Ednra/Sall1/Vps35/Amfr/Dixdc1/Csnk1g1/Tbx18/Amotl2/Ryk/Kank1/Mdfi/Tmem131l/Aspm/Ctnnd1/Mks1/Daam1/Pygo1/Ror1/Prickle1/Siah1a/Fermt2/Rnf146/Rbpj/Rbms3/Tax1bp3/Ankrd6/Daam2/Nppa/Zbed3/Disc1/Ppp2r3a/Tmem64/Foxo1/Dact1/Fzd5/Otulin/Csnk2a2/Zbtb33/Foxo3/Fzd4/Shisa3/Lgr4/Fzd2/Amer1/Adnp/Chd8/Cthrc1/Bmal1/Gpc3/Tcf7l1/Thra/Usp47/Eda/Cdh3/Bcl9l/Frat1/Scyl2/Apcdd1/Csnk1g3/Sox4/Tmem170b
## GO:0048568 Tcf7/Mmp14/Sec24b/Sp1/Col18a1/Akt1/Smo/Tead3/Prkra/Trim28/Kit/Ctnnb1/Arid1a/Fzd3/Rpl10/Hipk1/Bmp7/Kcnq1/Mrtfb/Bcr/Prox1/Arnt/Ift52/Kif3a/Pcgf2/Sp2/Hey2/Frs2/Cobl/Wdpcp/Id2/Llgl2/Ntn1/Plcd3/Hif1a/Zfp36l1/Vash1/Gli3/Ror2/Ptch1/Bmpr1a/Gata4/Nipbl/Gdnf/Zfpm2/Stk3/Myc/Scrib/Tbc1d23/Ptk7/Ift140/Axin1/Mib1/Smad2/Tcf7l2/Ldb1/Sufu/Vangl2/Hsd17b7/Itga8/Acvr1/Notch1/Nckap1/Itgav/Sp3/Vcam1/Col11a1/Pitx2/Mmp16/Stil/Syf2/Dnajb6/Pdgfra/Rest/Kdm2b/Plxna4/Tead2/Fgfr2/Ednra/Sall1/Gab1/Irx5/Tbx20/Tbx18/Tgfbr2/Wdr48/Mdfi/Ift57/Arid2/Mks1/Thoc5/Pbx2/Hectd1/Ripor2/Prickle1/Kdm1a/Kdm6a/Vash2/Wdr19/Rara/Hand2/Ncoa6/Pbx3/Rbpj/Tgfb2/Cc2d2a/Cited2/Bptf/Mthfd1l/Ppp1r13l/Erf/Setd2/Fzd5/Megf8/Tubb2b/Tcf21/Tshz1/C2cd3/Igf2/Krt8/Zfpm1/Fzd2/Foxc1/Gja1/Cited1/Pbx1/Nrk/Ndst1/Cthrc1/Slc8a1/Tead1/Gja5/Ece1/T2/Rnf207/Mfap2/Hipk2/Mapk3/Mapk1/Sox11/Zfp568/Spint2/Pdzd7/Ttbk2
## GO:0060070 Tcf7/Ddx3x/Sost/Jup/Nid1/Ctnnb1/Fzd3/Tle1/Gprc5b/Vps4b/Pten/Bicc1/Kremen1/Sdc1/Ppm1a/Gli3/Ror2/Ednrb/Sema5a/Stk3/Csnk1e/Gsk3b/Lmbr1l/Ptk7/Axin1/Ccny/Cdh2/Smad4/Csnk1a1/Gnaq/Tnks2/Lrp5/Tcf7l2/Csnk1d/Limd1/Tmem9/Notch1/Lrp4/Tbl1xr1/Src/Sfrp2/Mllt3/Cdk14/Ctnnbip1/Mad2l2/Ptpro/Caprin2/Dkk3/Fgfr2/Adgra2/Tnks/Sfrp1/Ednra/Vps35/Amfr/Dixdc1/Csnk1g1/Tbx18/Kank1/Tmem131l/Aspm/Ctnnd1/Mks1/Pygo1/Prickle1/Siah1a/Rnf146/Rbpj/Rbms3/Ankrd6/Daam2/Nppa/Zbed3/Disc1/Ppp2r3a/Tmem64/Foxo1/Dact1/Fzd5/Otulin/Foxo3/Fzd4/Shisa3/Lgr4/Fzd2/Amer1/Adnp/Chd8/Cthrc1/Bmal1/Gpc3/Tcf7l1/Thra/Usp47/Eda/Cdh3/Bcl9l/Frat1/Scyl2/Csnk1g3/Sox4/Tmem170b
## GO:0014855 Gata6/Ctnnb1/Pten/Rxra/Tbx5/Hey2/Cdk1/Sirt1/Prkar1a/Myh10/Bmpr1a/Gata4/Sugt1/Zfpm2/Angpt1/Cxadr/Pim1/Ndc80/Jak2/Notch1/Dsn1/Tgfbr3/Foxp1/Vgll4/Fgfr2/Smad1/Tbx20/Tgfbr2/Arid2/Rbpj/Tgfb2/Ncam1/Rxrb/Cited2/Ccnb1/Tenm4/Foxc1/Gja1/Erbb4/Nrg1/Mapk1
## GO:0021915 Tcf7/Sec24b/Smo/Prkacb/Arid1a/Fzd3/Bmp7/Ift52/Kif3a/Apaf1/Cobl/Wdpcp/Rock2/Hif1a/Zfp36l1/Gli3/Ptch1/Stk3/Scrib/Ptk7/Ift140/Mib1/Brd2/Kif20b/Tcf7l2/Sufu/Bbs4/Vangl2/Tsc1/Notch1/Nckap1/Traf6/Sfrp2/Lmo4/Stil/Ski/Glmn/Shroom3/Kdm2b/Tead2/Dlc1/Sfrp1/Casp3/Sall1/Ift57/Rpgrip1l/Specc1l/Mks1/Hectd1/Prickle1/Plxnb2/Kdm6a/Dzip1l/Wdr19/Rara/Tbc1d32/Tctn1/Tgfb2/Nup133/Cc2d2a/Cited2/Mthfd1l/Dact1/Setd2/C2cd3/Fzd2/Cthrc1/Arhgap35/Deaf1/Ssbp3/Spint2/Ttbk2
## GO:0003205 Smo/Ltbp1/Gata6/Epor/Tek/Ctnnb1/Arid1a/Fgfrl1/Bmp7/Prox1/Rxra/Tbx5/Hey2/Frs2/Myocd/Id2/Myh10/Hif1a/Ptcd2/Bmpr1a/Gata4/Zfpm2/Robo1/Adamts1/Ptk7/Smad4/Sufu/Tnni1/Vangl2/Acvr1/Notch1/Col11a1/Sfrp2/Pitx2/Lmo4/Ube4b/Tgfbr3/Parva/Fgfr2/Slit2/Ednra/Sall1/Tbx20/Tgfbr2/Mks1/Ap2b1/Hectd1/Matr3/Hand2/Adgrg6/Rbpj/Tgfb2/Cplane1/Cited2/Ift88/Ppp1r13l/Pkp2/Greb1l/Adamts6/Zfpm1/Fzd2/Foxc1/Maml1/Robo2/Myh7/Ndst1/Dsp/Slit3/Cntrl/Gja5/Nsd2/Nrg1/Sox11/Bmpr2/Sox4
## Count
## GO:0048598 187
## GO:0060038 35
## GO:0003007 102
## GO:0198738 141
## GO:0016055 140
## GO:0048568 144
## GO:0060070 102
## GO:0014855 41
## GO:0021915 72
## GO:0003205 75
# Plot 1:20 all terms
go_terms_all_upreg <- ggplot(egobp_up[1:20], aes(x=reorder(Description, -pvalue), y=Count, fill=-p.adjust)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_continuous(low="#b3b3b3", high="lightsteelblue2") +
labs(x = "", y = "", fill = "p.adjust") +
theme(axis.text=element_text(size=11),
colour = "black") +
theme_minimal() +
ggtitle(label = "Up-regulated",
subtitle = "Biological Process") +
ylab("Hit Proteins") +
theme(axis.text.y=element_text(colour = "black"))
go_terms_all_upreg
# Filter the enrichment results for the GO term "angiogenesis"
go_term_data <- BP_enrichment_results_upreg[grep("angiogenesis", BP_enrichment_results_upreg$Description, ignore.case = TRUE), ]
# Extract the gene list associated with this GO term
if (nrow(go_term_data) > 0) {
gene_list <- strsplit(go_term_data$geneID, "/")[[1]]
# View the gene list
print(gene_list)
# Save the gene list to a CSV file
write.csv(gene_list, "angiogenesis_gene_list.csv", row.names = FALSE)
} else {
print("No GO term related to angiogenesis found in the results.")
}
## [1] "Itga5" "Sp1" "Efnb2" "Col18a1" "Jup" "Akt1"
## [7] "Hif3a" "Hdac9" "Klf5" "Hmox1" "Gata6" "Tek"
## [13] "Ctnnb1" "Elk3" "Hdac5" "Vav2" "Vps4b" "Amotl1"
## [19] "Pten" "Sash1" "Rxra" "Sulf1" "Plcg1" "Mmp9"
## [25] "Akt3" "Ccn2" "Sirt1" "Rtn4" "Rock2" "Gna13"
## [31] "Pxdn" "Map3k3" "Plcd3" "Hif1a" "Vash1" "Syk"
## [37] "Smad5" "Tert" "Plk2" "Bmpr1a" "Gata4" "Npr3"
## [43] "Sema5a" "Mtdh" "Angpt1" "Robo1" "Adamts1" "Nr4a1"
## [49] "Thbs2" "Cyp1b1" "C3" "Cemip2" "Hgs" "Fn1"
## [55] "Ncl" "Cfh" "Glul" "Vangl2" "Acvr1" "Notch1"
## [61] "Itgav" "Atf2" "Tspan18" "Pik3ca" "Nras" "Sfrp2"
## [67] "Pitx2" "B4galt1" "Jak1" "Wasf2" "Add1" "Emilin1"
## [73] "Pdgfra" "Tgfbr3" "Ephb4" "Ptn" "Tgfa" "Adipor2"
## [79] "Tjp1" "Pak4" "Mfge8" "Pik3c2a" "Pde3b" "Parva"
## [85] "Fgfr2" "Srpx2" "Mecp2" "Adgra2" "Slit2" "Ednra"
## [91] "Smad1" "Gab1" "Mmp2" "Cdh5" "Nfatc3" "Tbx20"
## [97] "Mcam" "Rora" "Tgfbr2" "Ngp" "Amotl2" "Thsd7a"
## [103] "Flcn" "Vav3" "Jcad" "Ramp1" "Grn" "Prkx"
## [109] "Btg1" "Ago2" "Vash2" "Hand2" "Rbpj" "Tgfb2"
## [115] "Htatip2" "Thbs1" "Ppp1r15a" "Cldn5" "Ago1" "C1galt1"
## [121] "Shc1" "Foxo4" "Setd2" "Shb" "Fzd5" "S1pr1"
## [127] "Tcf21" "Otulin" "Igf2" "Agtr1a" "Fzd4" "Foxc1"
## [133] "Lgals3" "Jun" "Adam12" "Mia3" "Lepr" "Map2k5"
## [139] "Bcas3" "Ptk2b" "Gtf2i" "Prcp" "Hipk2" "Nrxn3"
## [145] "Apold1"
# Wnt signaling pathway
go_term_data_wnt <- BP_enrichment_results_upreg[grep("Wnt", BP_enrichment_results_upreg$Description, ignore.case = TRUE), ]
# Extract the gene list associated with this GO term
if (nrow(go_term_data_wnt) > 0) {
gene_list_wnt <- strsplit(go_term_data_wnt$geneID, "/")[[1]]
# View the gene list
print(gene_list_wnt)
# Save the gene list to a CSV file
write.csv(gene_list_wnt, "wnt_gene_list.csv", row.names = FALSE)
} else {
print("No GO term related to angiogenesis found in the results.")
}
## [1] "Tcf7" "Ddx3x" "Sost" "Itga3" "Jup" "Hbp1"
## [7] "Grk5" "Cul3" "Wwox" "Ccn4" "Crbn" "Nid1"
## [13] "Ctnnb1" "Fzd3" "Tle1" "Gprc5b" "Vps4b" "Amotl1"
## [19] "Pten" "Bicc1" "Sulf1" "Grb10" "Wif1" "Fbxw11"
## [25] "Kremen1" "Sdc1" "Ppm1a" "Gli3" "Spin1" "Ror2"
## [31] "Tert" "Ednrb" "Sema5a" "Stk3" "Myc" "Csnk1e"
## [37] "Gsk3b" "Lmbr1l" "Calcoco1" "Ptk7" "Strn" "Axin1"
## [43] "Ccny" "Cdh2" "Smad4" "Csnk1a1" "Gnaq" "Ddb1"
## [49] "Tnks2" "Lrp5" "Tcf7l2" "Csnk1d" "Ldb1" "Limd1"
## [55] "Tmem9" "Klhl12" "Vangl2" "Notch1" "Lrp4" "Rtf1"
## [61] "Tbl1xr1" "Src" "Sfrp2" "Pitx2" "Mllt3" "Prkaa2"
## [67] "Cdk14" "Ctnnbip1" "Mad2l2" "Klf15" "Ptpro" "Caprin2"
## [73] "Vgll4" "Dkk3" "Fgfr2" "Zranb1" "Adgra2" "Tnks"
## [79] "Sfrp1" "Ednra" "Sall1" "Vps35" "Amfr" "Dixdc1"
## [85] "Csnk1g1" "Tbx18" "Amotl2" "Ryk" "Kank1" "Mdfi"
## [91] "Tmem131l" "Aspm" "Ctnnd1" "Mks1" "Daam1" "Pygo1"
## [97] "Ror1" "Ptpn23" "Prickle1" "Siah1a" "Fermt2" "Rnf146"
## [103] "Rbpj" "Rbms3" "Tax1bp3" "Ankrd6" "Daam2" "Nppa"
## [109] "Zbed3" "Disc1" "Ppp2r3a" "Tmem64" "Foxo1" "Dact1"
## [115] "Fzd5" "Otulin" "Csnk2a2" "Zbtb33" "Foxo3" "Fzd4"
## [121] "Shisa3" "Lgr4" "Fzd2" "Amer1" "Adnp" "Chd8"
## [127] "Cthrc1" "Bmal1" "Gpc3" "Tcf7l1" "Thra" "Usp47"
## [133] "Eda" "Cdh3" "Bcl9l" "Frat1" "Scyl2" "Apcdd1"
## [139] "Csnk1g3" "Sox4" "Tmem170b"
# BMP signaling pathway
go_term_data_bmp <- BP_enrichment_results_upreg[grep("Bmp", BP_enrichment_results_upreg$Description, ignore.case = TRUE), ]
# Extract the gene list associated with this GO term
if (nrow(go_term_data_bmp) > 0) {
gene_list_bmp <- strsplit(go_term_data_bmp$geneID, "/")[[1]]
# View the gene list
print(gene_list_bmp)
# Save the gene list to a CSV file
write.csv(gene_list_bmp, "bmp_gene_list.csv", row.names = FALSE)
} else {
print("No GO term related to angiogenesis found in the results.")
}
## [1] "Sost" "Itga3" "Gata6" "Bmp7" "Sulf1" "Smurf2" "Ppm1a"
## [8] "Hivep1" "Ror2" "Smad5" "Bmpr1a" "Gata4" "Fstl1" "Smad4"
## [15] "Smad2" "Pdcd4" "Tcf7l2" "Acvr1" "Notch1" "Atf2" "Fbn1"
## [22] "Sfrp2" "Ski" "Tgfbr3" "Usp9x" "Sfrp1" "Smad1" "Comp"
## [29] "Cdh5" "Tbx20" "Neo1" "Trim33" "Ctdspl2" "Scx" "Fstl4"
## [36] "Rbpj" "Megf8" "Acvr2a" "Gpc3" "Hipk2" "Mapk3" "Sox11"
## [43] "Numa1" "Bmpr2"
# Genes from GO-terms of interest (bmp)
# Filter results for these genes
filtered_lrt_df_bmp <- lrt_df[lrt_df$SYMBOL %in% gene_list_bmp[1:10], ]
# Plot logFC values with bars for reference and treatment
ggplot(filtered_lrt_df_bmp, aes(x = SYMBOL, y = logFC)) +
geom_bar(stat = "identity", position = "dodge", fill = "lightsteelblue2") + # Dodge position to place bars side by side
labs(x = "Gene", y = "Log Fold Change (logFC)", title = "Bmp Terms") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
# Filter results for these genes
filtered_lrt_df_bmp_2 <- lrt_df[lrt_df$SYMBOL %in% gene_list_bmp[11:20], ]
# Plot logFC values with bars for reference and treatment
ggplot(filtered_lrt_df_bmp_2, aes(x = SYMBOL, y = logFC)) +
geom_bar(stat = "identity", position = "dodge", fill = "lightsteelblue2") + # Dodge position to place bars side by side
labs(x = "Gene", y = "Log Fold Change (logFC)", title = "Bmp Terms") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
# Genes from GO-terms of interest (wnt)
# Filter results for these genes
filtered_lrt_df_wnt <- lrt_df[lrt_df$SYMBOL %in% gene_list_wnt[1:10], ]
# Plot logFC values with bars for reference and treatment
ggplot(filtered_lrt_df_wnt, aes(x = SYMBOL, y = logFC)) +
geom_bar(stat = "identity", position = "dodge", fill = "lightsteelblue2") + # Dodge position to place bars side by side
labs(x = "Gene", y = "Log Fold Change (logFC)", title = "Wnt Terms") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
# Genes from GO-terms of interest (angiogenesis)
# Filter results for these genes
filtered_lrt_df_angio <- lrt_df[lrt_df$SYMBOL %in% gene_list[1:10], ]
# Plot logFC values with bars for reference and treatment
ggplot(filtered_lrt_df_angio, aes(x = SYMBOL, y = logFC)) +
geom_bar(stat = "identity", position = "dodge", fill = "lightsteelblue2") + # Dodge position to place bars side by side
labs(x = "Gene", y = "Log Fold Change (logFC)", title = "Angiogenesis Terms") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
axis.title.x = element_text(face = "italic"),
axis.title.y = element_text(face = "italic"),
axis.text = element_text(size = 12)
)
# Plot 1:20 all terms
go_terms_all_downreg <- ggplot(egobp_down[1:20], aes(x=reorder(Description, -pvalue), y=Count, fill=-p.adjust)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_continuous(low="#b3b3b3", high="lightsteelblue2") +
labs(x = "", y = "", fill = "p.adjust") +
theme(axis.text=element_text(size=11),
colour = "black") +
theme_minimal() +
ggtitle(label = "Down-regulated",
subtitle = "Biological Process") +
ylab("Hit Proteins") +
theme(axis.text.y=element_text(colour = "black"))
go_terms_all_downreg
# Package versions
sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] clusterProfiler_4.12.2 enrichplot_1.24.2 ggridges_0.5.6
## [4] gplots_3.1.3.1 lubridate_1.9.3 forcats_1.0.0
## [7] stringr_1.5.1 purrr_1.0.2 tidyr_1.3.1
## [10] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
## [13] edgeR_4.2.1 limma_3.60.4 org.Hs.eg.db_3.19.1
## [16] biomaRt_2.60.1 org.Mm.eg.db_3.19.1 dplyr_1.1.4
## [19] GenomicFeatures_1.56.0 AnnotationDbi_1.66.0 Biobase_2.64.0
## [22] GenomicRanges_1.56.1 GenomeInfoDb_1.40.1 IRanges_2.38.1
## [25] S4Vectors_0.42.1 BiocGenerics_0.50.0 readr_2.1.5
## [28] tximport_1.32.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.4.1 BiocIO_1.14.0
## [3] bitops_1.0-8 ggplotify_0.1.2
## [5] filelock_1.0.3 polyclip_1.10-7
## [7] XML_3.99-0.17 lifecycle_1.0.4
## [9] httr2_1.0.2 vroom_1.6.5
## [11] lattice_0.22-6 MASS_7.3-60.2
## [13] magrittr_2.0.3 sass_0.4.9
## [15] rmarkdown_2.27 jquerylib_0.1.4
## [17] yaml_2.3.10 cowplot_1.1.3
## [19] DBI_1.2.3 RColorBrewer_1.1-3
## [21] abind_1.4-5 zlibbioc_1.50.0
## [23] ggraph_2.2.1 RCurl_1.98-1.16
## [25] yulab.utils_0.1.5 tweenr_2.0.3
## [27] rappdirs_0.3.3 GenomeInfoDbData_1.2.12
## [29] ggrepel_0.9.5 tidytree_0.4.6
## [31] codetools_0.2-20 DelayedArray_0.30.1
## [33] DOSE_3.30.2 xml2_1.3.6
## [35] ggforce_0.4.2 tidyselect_1.2.1
## [37] aplot_0.2.3 UCSC.utils_1.0.0
## [39] farver_2.1.2 viridis_0.6.5
## [41] matrixStats_1.3.0 BiocFileCache_2.12.0
## [43] GenomicAlignments_1.40.0 jsonlite_1.8.8
## [45] tidygraph_1.3.1 tools_4.4.1
## [47] progress_1.2.3 treeio_1.28.0
## [49] snow_0.4-4 Rcpp_1.0.13
## [51] glue_1.7.0 gridExtra_2.3
## [53] SparseArray_1.4.8 xfun_0.46
## [55] qvalue_2.36.0 MatrixGenerics_1.16.0
## [57] withr_3.0.1 fastmap_1.2.0
## [59] fansi_1.0.6 caTools_1.18.2
## [61] digest_0.6.36 timechange_0.3.0
## [63] R6_2.5.1 gridGraphics_0.5-1
## [65] colorspace_2.1-1 GO.db_3.19.1
## [67] gtools_3.9.5 RSQLite_2.3.7
## [69] utf8_1.2.4 generics_0.1.3
## [71] data.table_1.15.4 rtracklayer_1.64.0
## [73] prettyunits_1.2.0 graphlayouts_1.1.1
## [75] httr_1.4.7 S4Arrays_1.4.1
## [77] scatterpie_0.2.3 pkgconfig_2.0.3
## [79] gtable_0.3.5 blob_1.2.4
## [81] XVector_0.44.0 shadowtext_0.1.4
## [83] htmltools_0.5.8.1 fgsea_1.30.0
## [85] scales_1.3.0 png_0.1-8
## [87] ggfun_0.1.5 knitr_1.48
## [89] rstudioapi_0.16.0 tzdb_0.4.0
## [91] reshape2_1.4.4 rjson_0.2.21
## [93] nlme_3.1-164 curl_5.2.1
## [95] cachem_1.1.0 KernSmooth_2.23-24
## [97] parallel_4.4.1 HDO.db_0.99.1
## [99] restfulr_0.0.15 pillar_1.9.0
## [101] grid_4.4.1 vctrs_0.6.5
## [103] dbplyr_2.5.0 evaluate_0.24.0
## [105] cli_3.6.3 locfit_1.5-9.10
## [107] compiler_4.4.1 Rsamtools_2.20.0
## [109] rlang_1.1.4 crayon_1.5.3
## [111] labeling_0.4.3 plyr_1.8.9
## [113] fs_1.6.4 stringi_1.8.4
## [115] viridisLite_0.4.2 BiocParallel_1.38.0
## [117] munsell_0.5.1 Biostrings_2.72.1
## [119] lazyeval_0.2.2 GOSemSim_2.30.0
## [121] Matrix_1.7-0 hms_1.1.3
## [123] patchwork_1.2.0 bit64_4.0.5
## [125] KEGGREST_1.44.1 statmod_1.5.0
## [127] highr_0.11 SummarizedExperiment_1.34.0
## [129] igraph_2.0.3 memoise_2.0.1
## [131] bslib_0.8.0 ggtree_3.12.0
## [133] fastmatch_1.1-4 bit_4.0.5
## [135] gson_0.1.0 ape_5.8