genexsea

[genexsea]

Gene Expression and Enrichment Analysis (R code)

library(edgeR) library(limma) library(Glimma) library(gplots) library(org.Mm.eg.db) library(RColorBrewer) library(RnaSeqGeneEdgeRQL) library(statmod) library(sva) library(bladderbatch) library(pamr) library(limma) library(ggplot2) library(gridExtra) library(plotly) library(tidyverse) library(gridExtra) library(grid) library(png) library(downloader) library(grDevices) library(cowplot) library(ggplot2)

library(snpStats) getwd()

#Load the Sample information, Batches and the conditions. This will correspond to the Phenotype data. #The variables are stored in the phenotype data slot and can be obtained as follows:

sampleinfo <- read.delim("Desktop/Abhilash/study/Post-Doc/Lab_notebook/eQTL/Input_3/targets_red_blue.txt",stringsAsFactors = default.stringsAsFactors(),row.names = "samples")

#Load the expression Data / Reading in count-data #The expression data can be obtained from count matrix as follows: seqdata_all <- read.delim("Desktop/Abhilash/study/Post-Doc/Lab_notebook/eQTL/Input_3/FP_TR_red_blue.txt",                          row.names="ORF19_ID",stringsAsFactors = FALSE)

#seqdata_all <- seqdata_all[,c(1:36)]

#Removing the 'Gene Length' column from the Expression dataset

seqdata_all <- seqdata_all[,-c(1)]

#Check if the row names of 'sampleinfo' is the same as the column names of 'seqdata' stopifnot(all(rownames(sampleinfo) == colnames(seqdata_all)))

group_isolate <- factor(sampleinfo$isolate, levels = unique(sampleinfo$isolate)) group_isolate

group_type <- factor(sampleinfo$culture_type, levels = unique(sampleinfo$culture_type)) group_type

group_plate <- factor(sampleinfo$plate_info, levels = unique(sampleinfo$plate_info)) group_plate

group_flowcell <- factor(sampleinfo$Flowcell, levels = unique(sampleinfo$Flowcell)) group_flowcell

#Converting the read data (Expression data) into a DGEList-object using the DGEList function.

x <- DGEList(seqdata_all) class(x) dim(x) x colnames(x)

# Updating the group and batch information to the samples in DGElist x$samples$group_isolate <- group_isolate x$samples$group_type <- group_type x$samples$group_plate <- group_plate x$samples$group_flowcell <- group_flowcell

x

################################################################################################

#####################################Adding annotations to the DGElist############################### # load the file Containing the Gene names and Gene Ids of all the annotated Candida genes

IDs <- read.delim("Desktop/Abhilash/study/Post-Doc/Lab_notebook/eQTL/Input_3/C_albicans_genes_allids_desc.txt", row.names="ORF19_ID") #Merging the ID dataframe to DGElist dataframe containing the gene details x$genes <- merge(x$genes, IDs, by="row.names", all=T)

#Convert the values in a "orf19" column into row names in DGElist existing data frame containing gene information

samp2 <- x$genes[,-1]

x$genes <- samp2

head(x$genes) dim(x)

###########################################################################################

#################################Pre-processing the data################

#Transforming the raw read counts to counts per million (cpm) #This transformation accounts for the differences in the library size #CPM and log-CPM transformations do not account for gene length differences as RPKM and FPKM values do #raw counts are converted to CPM and log-CPM values using the cpm function in edgeR.

cpm_x <- cpm(x) lcpm_x <- cpm(x, log=TRUE, prior.count=2) head(cpm_x)[,1] head(lcpm_x)[,1] dim(lcpm_x)

L_x <- mean(x$samples$lib.size) * 1e-6 M_x <- median(x$samples$lib.size) * 1e-6 c(L_x, M_x)

# minimum log-CPM value for each sample becomes log2(2/3.923991) = -3.18 summary(lcpm_x)

#######################################################################

######################Removing lowly expressed Candida genes##################################

table(rowSums(x$counts==0)==20)

#############################################################################################

#########################Filtering low-expressed genes######################################

#The **filterByExpr** function in the edgeR package provides an automatic way to filter genes, while keeping as many genes as possible with worthwhile counts. keep.exprs_x <- filterByExpr(x, group=group_isolate) x <- x[keep.exprs_x,, keep.lib.sizes=FALSE] dim(x)

#Plotting the distribution log-CPM values shows that a sizeable proportion of genes within each sample are either unexpressed or lowly-expressed with log-CPM values that are small or negative lcpm.cutoff_x <- log2(10/M_x + 2/L_x) library(RColorBrewer) nsamples_x <- ncol(x) col_x <- brewer.pal(nsamples_x, "Paired") par(mfrow=c(1,2)) plot(density(lcpm_x[,1]), col=col_x[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title(main="A. Raw data", xlab="Log-cpm") abline(v=lcpm.cutoff_x, lty=) for (i in 2:nsamples_x){  den <- density(lcpm_x[,i])  lines(den$x, den$y, col=col_x[i], lwd=2) } legend("topright", colnames(x), text.col=col_x, bty="n", cex = 0.3, ncol=2)

#The density of log-CPM values for post-filtered data are shown for each sample

lcpm_x <- cpm(x, log=TRUE) head(lcpm_x)[,1] plot(density(lcpm_x[,1]), col=col_x[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title(main="B. Filtered data", xlab="Log-cpm") abline(v=lcpm.cutoff_x, lty=3) for (i in 2:nsamples_x){  den <- density(lcpm_x[,i])  lines(den$x, den$y, col=col_x[i], lwd=2) } legend("topright", colnames(x), text.col=col_x, bty="n", cex = 0.3, ncol=2)

##################################################################################################

######################Qulaity control steps##################################################

#### Determining the sample library sizes #First, we can check how many reads we have for each sample in **x**(DGElist-object) x$samples$lib.size

#We can also plot the library sizes as a barplot to see whether there are any major discrepancies between the samples more easily. # The names argument tells the barplot to use the sample names on the x-axis # The las argument rotates the axis names

par(mfrow=c(1,1)) barplot(x$samples$lib.size,names=colnames(x),las=2, cex.axis = 0.8,cex.names = 0.3) # Add a title to the plot title("Barplot of library sizes")

#### Normalising gene expression distributions for compositional Bias x <- calcNormFactors(x, method = "TMM") x$samples$norm.factors x

#####################Observing the effects of normalization##########################

# x2 <- x # x2$samples$norm.factors <- 1 # x2$counts[,1] <- ceiling(x2$counts[,1]*0.05) # x2$counts[,2] <- x2$counts[,2]*5 # # lcpm_x <- cpm(x2, log=TRUE,prior.count=2) # head(lcpm_x)[,1] # par(mfrow=c(1,2)) # boxplot(lcpm_x, las=2, col=col_x, main="") # title(main="A. Example: Unnormalised data",ylab="Log-cpm") # x2 <- calcNormFactors(x2)   # x2$samples$norm.factors # # lcpm_x <- cpm(x2, log=TRUE,prior.count=2) # head(lcpm_x)[,1] # boxplot(lcpm_x, las=2, col=col_x, main="") # title(main="B. Example: Normalised data",ylab="Log-cpm")

lcpm_x <- cpm(x, log=TRUE, prior.count = 2) head(lcpm_x)[,1]

############Adjusting for the batch effects###########################33

#Simple way

mod = model.matrix(~as.factor(group_type) + as.factor(group_flowcell)) fit = lm.fit(mod, t(lcpm_x)) hist(fit$coefficients[2,], col=2,breaks=100)

table(group_isolate, group_flowcell)

# using combat

batch = group_flowcell modcombat = model.matrix(~1, data=sampleinfo) modtype = model.matrix(~0+group_type) combat_lcpm = ComBat(dat=lcpm_x, batch = batch, mod=modcombat,par.prior = T, prior.plots = FALSE) combat_fit = lm.fit(modtype, t(combat_lcpm)) hist(combat_fit$coefficients[2,], col=2, breaks = 100)

#plotting coefficients from the original liner model vs Combat fit

plot(fit$coefficients[2,],combat_fit$coefficients[2,], col=2, xlab="Linear Model",     ylab="Combat", xlim=c(-5,5), ylim=c(-5,5))

abline(c(0,1),col=1, lwd=3)

par(mfrow=c(1,1))

#Using SVA package to infer the batch effects

mod=model.matrix(~0 + group_type, data=sampleinfo) # Build model matrix with the outcome you care about mod0=model.matrix(~1, data=sampleinfo)  # Null model n.sv = num.sv(lcpm_x,mod)

sva1=sva(lcpm_x,mod,mod0) names(sva1) dim(sva1$sv)

#correlating batch effects with the observed with actual observed batch variable

summary(lm(sva1$sv ~ sampleinfo$Flowcell)) boxplot(sva1$sv[,4] ~ sampleinfo$Flowcell)

points(sva1$sv[,4] ~ jitter(as.numeric(sampleinfo$Flowcell)), col=as.numeric(sampleinfo$Flowcell)) modsv = cbind(mod,sva1$sv)

fitsv = lm.fit(modsv, t(lcpm_x))

# Comparing SVA and Combat par(mfrow=c(1,2))

plot(fitsv$coefficients[2,],combat_fit$coefficients[2,],col=2,     xlab="SVA", ylab="Combat", xlim=c(-5,5), ylim=c(-5,5)) abline(c(0,1), col=1, lwd=3)

# Comparing SVA and liner adjusted model

plot(fitsv$coefficients[2,],fit$coefficients[2,],col=2,     xlab="SVA", ylab="linear model", xlim=c(-5,5), ylim=c(-5,5)) abline(c(0,1), col=1, lwd=3)

cleaningP = function(x, mod, sva1,  P=ncol(mod)) {  X=cbind(mod,sva1$sv)  Hat=solve(t(X)%*%X)%*%t(X)  beta=(Hat%*%t(x))  cleany=x-t(as.matrix(X[,-c(1:P)])%*%beta[-c(1:P),])  return(cleany) }

sva_lcpm = cleaningP(lcpm_x, mod, sva1) sva_lcpm

#write.csv(sva_lcpm, "Desktop/Abhilash/study/Post-Doc/Lab_notebook/eQTL/Input_3/sva_lcpm.csv")

################Checking for the variation between sample with PCA plot

par(mfrow=c(1,2))

pca_lcpm_x = prcomp(t(lcpm_x), scale=FALSE) plot(pca_lcpm_x$x[,c(1,2)],col=sampleinfo$isolate,pch=16) text(pca_lcpm_x$x[,c(1,2)],labels=sampleinfo$isolate,pos=3,cex=0.5) title("PC1 vs PC2- No correction for batch effect")

pca_combat_lcpm = prcomp(t(combat_lcpm), scale=FALSE) plot(pca_combat_lcpm$x[,c(1,2)],col=sampleinfo$isolate,pch=16) text(pca_combat_lcpm$x[,c(1,2)],labels=sampleinfo$isolate,pos=3,cex=0.5) title("PC1 vs PC2- Correction with Combat")

pca_sva_lcpm = prcomp(t(sva_lcpm), scale=FALSE) plot(pca_sva_lcpm$x[,c(1,2)],col=sampleinfo$isolate,pch=16) text(pca_sva_lcpm$x[,c(1,2)],labels=sampleinfo$isolate,pos=3,cex=0.5) title("PC1 vs PC2 - Correction with SVA")

par(mfrow=c(1,2))

hist(lcpm_x[,2], col=2, breaks = 100)

hist(sva_lcpm[,2], col=2, breaks = 100)

plot(density(lcpm_x[,2]))

lines(density(sva_lcpm[,2]),col="green")

mod_type=model.matrix(~0+group_type) v <- voom(x, mod_type, plot=TRUE) v col.group <- group_type nb.cols <- 2 levels(col.group) <-  colorRampPalette(brewer.pal(nlevels(col.group), "Set1"))(nb.cols) col.group <- as.character(col.group) #pch <- c(0,0,1,1,2,2,3,3,3,3,3,3,4,4,5,5,6,6,7,7,8,8,8,8,9,9,10,10,11,11,11,11,11,11,11)

plotMDS(v$E, labels=group_type, col=col.group, cex = 0.5) #plotMDS(v$E, col=col.group, pch=pch, cex = 0.4)

col.group <- group_type nb.cols <- 2 levels(col.group) <-  colorRampPalette(brewer.pal(nlevels(col.group), "Set1"))(nb.cols) col.group <- as.character(col.group)

plotMDS(v$E, labels=group_type, col=col.group, cex = 0.5)

#Using SVA package to infer the batch effects

mod=model.matrix(~0+group_type) # Build model matrix with the outcome you care about mod0=model.matrix(~1, data=sampleinfo)  # Null model n.sv = num.sv(v$E,mod)

sva1=sva(v$E,mod,mod0) names(sva1) dim(sva1$sv)

#correlating batch effects with the observed with actual observed batch variable

summary(lm(sva1$sv ~ sampleinfo$Flowcell)) boxplot(sva1$sv[,1] ~ sampleinfo$Flowcell)

points(sva1$sv[,1] ~ jitter(as.numeric(sampleinfo$Flowcell)), col=as.numeric(sampleinfo$Flowcell)) modsv = cbind(mod,sva1$sv)

fitsv = lm.fit(modsv, t(v$E))

cleaningP = function(x, mod, sva1,  P=ncol(mod)) {  X=cbind(mod,sva1$sv)  Hat=solve(t(X)%*%X)%*%t(X)  beta=(Hat%*%t(x))  cleany=x-t(as.matrix(X[,-c(1:P)])%*%beta[-c(1:P),])  return(cleany) }

counts_sva = cleaningP(v$E, mod, sva1) counts_sva

#combat_v_E = ComBat(dat=v$E, batch = batch, mod=modcombat,par.prior = T, prior.plots = FALSE) #combat_fit = lm.fit(modisolate, t(combat_lcpm))

################Checking for the variation between sample with PCA plot

par(mfrow=c(1,2))

pca_V_E = prcomp(t(v$E), scale=FALSE) plot(pca_V_E$x[,c(1,2)],col=sampleinfo$isolate,pch=16) text(pca_V_E$x[,c(1,2)],labels=sampleinfo$isolate,pos=3,cex=0.5) title("PC1 vs PC2- No correction for batch effect")

pca_counts_sva = prcomp(t(counts_sva), scale=FALSE) plot(pca_counts_sva$x[,c(1,2)],col=sampleinfo$isolate,pch=16) text(pca_counts_sva$x[,c(1,2)],labels=sampleinfo$isolate,pos=3,cex=0.5) title("PC1 vs PC2 - Correction with SVA")

# pca_combat_v_E = prcomp(t(combat_v_E), scale=FALSE) # plot(pca_combat_v_E$x[,c(1,2)],col=sampleinfo$isolate,pch=16) # text(pca_combat_v_E$x[,c(1,2)],labels=sampleinfo$isolate,pos=3,cex=0.5) # title("PC1 vs PC2 - Correction with Combat")

par(mfrow=c(1,2))

plot(density(lcpm_x[,1]), col=col_x[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title(main="A. Notcorrected for Batch effect", xlab="Log-cpm") abline(v=lcpm.cutoff_x, lty=) for (i in 2:nsamples_x){  den <- density(lcpm_x[,i])  lines(den$x, den$y, col=col_x[i], lwd=2) }

legend("topright", colnames(x), text.col=col_x, bty="n", cex = 0.3, ncol=2)

plot(density(v$E[,1]), col=col_x[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title(main="B. Batch_corrected", xlab="Voom-transformed(Log-cpm)") abline(v=lcpm.cutoff_x, lty=) for (i in 2:nsamples_x){  den <- density(v$E[,i])  lines(den$x, den$y, col=col_x[i], lwd=2) }

legend("topright", colnames(x), text.col=col_x, bty="n", cex = 0.3, ncol=2)

plot(density(counts_sva[,1]), col=col_x[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title(main="B. Batch_corrected", xlab="Voom-transformed(Log-cpm)") abline(v=lcpm.cutoff_x, lty=) for (i in 2:nsamples_x){  den <- density(counts_sva[,i])  lines(den$x, den$y, col=col_x[i], lwd=2) }

legend("topright", colnames(x), text.col=col_x, bty="n", cex = 0.3, ncol=2)

colnames(mod_type) <- gsub("group_type", "", colnames(mod_type)) mod_type

#Creating contrast matrix contr.matrix_type <- makeContrasts(  MonoculturevsCoculture =Monoculture-Coculture,  levels = colnames(mod_type)) contr.matrix_type

#Replacing the corrected expression values to the voom objects

v[["E"]] <- counts_sva

vfit_nosv <- lmFit(v, mod_type) vfit_nosv <- contrasts.fit(vfit_nosv, contrasts=contr.matrix_type) efit_nosv <- eBayes(vfit_nosv) plotSA(efit_nosv, main="Model after batch effect removal: Mean-variance trend") de_nosv <- decideTests(efit_nosv,lfc = 0.3, genewise.p.value = 0.05) summary(decideTests(efit_nosv))

tfit_nosv <- treat(vfit_nosv, lfc=1) dt_nosv <- decideTests(tfit_nosv) summary(dt_nosv)

FPvsTR <- topTable(efit_nosv, coef=1, number=nrow(v$E))

write.csv(FPvsTR, "Desktop/Abhilash/study/Post-Doc/Presentations/TIM/redvsblue/results/FPvsTR.csv")

glMDPlot(efit_nosv, coef=1, status=de_nosv, main=colnames(efit_nosv)[1],         side.main="Genename", counts=v$E, groups=group_type, launch=TRUE,ylim=c(--10,13))

library(edgeR) library(limma) library(Glimma) library(gplots) library(org.Mm.eg.db) library(RColorBrewer) library(RnaSeqGeneEdgeRQL) library(statmod) library(sva) library(bladderbatch) library(pamr) library(limma) library(ggplot2) library(gridExtra) library(plotly) library(tidyverse) library(gridExtra) library(grid) library(png) library(downloader) library(grDevices) library(cowplot) library(ggplot2) library(DOSE)

library(AnnotationHub) hub <- AnnotationHub()

query(hub, "Candida albicans") Candida_alb <- hub[["AH73663"]] Candida_alb

library(clusterProfiler) FPvsTR_enrich_sig_enrich = read.csv("Desktop/Abhilash/study/Post-Doc/Presentations/TIM/redvsblue/results/FPvsTR.csv") FPvsTR_enrich_sig <- FPvsTR_enrich_sig_enrich %>% filter(FPvsTR_enrich_sig_enrich$adj.P.Val < 0.05) FPvsTR_enrich_sig_down <- FPvsTR_enrich_sig %>% filter(FPvsTR_enrich_sig$logFC < -0.5) dim(FPvsTR_enrich_sig_down)

#Top_upregulated #FPvsTR_significant_Upregulated (FDR<0.05 and Lfc > 0.5) FPvsTR_enrich_sig_up <- FPvsTR_enrich_sig %>% filter(FPvsTR_enrich_sig$logFC > 0.5) dim(FPvsTR_enrich_sig_up)

#Significant_downregulated #FPvsTR_Downregulated_all (FDR<0.05 and Lfc <0) FPvsTR_enrich_sig_down_all <- FPvsTR_enrich_sig %>% filter(FPvsTR_enrich_sig$logFC < 0) dim(FPvsTR_enrich_sig_down_all)

#Significant_upregulated #FPvsTR_Upregulated_all (FDR<0.05 and Lfc > 0) FPvsTR_enrich_sig_up_all <- FPvsTR_enrich_sig %>% filter(FPvsTR_enrich_sig$logFC > 0) dim(FPvsTR_enrich_sig_up_all)

####code for producing the genelists

#Creating genelists for different comaparisons #1. FPvsTR geneList_FPvsTR_enrich_sig_logFC <- FPvsTR_enrich_sig$logFC names(geneList_FPvsTR_enrich_sig_logFC) <- FPvsTR_enrich_sig$EntrezID length(geneList_FPvsTR_enrich_sig_logFC)

### Gene Set enrichment analysis

#all terms together #CEC3609_3hrsvsCEC3609_6hrs #all significant (FDR < 0.05) gse_FPvsTR_enrich_sig_all <- gseGO(geneList=sort(geneList_FPvsTR_enrich_sig_logFC,decreasing = T),ont = "ALL", OrgDb=Candida_alb, verbose=F,pvalueCutoff = 1, nPerm = 10000) head(summary(gse_FPvsTR_enrich_sig_all)) write.csv(gse_FPvsTR_enrich_sig_all@result,file = "Desktop/Abhilash/study/Post-Doc/Presentations/TIM/redvsblue/results/gse_FPvsTR_enrich_sig_all.csv")

#significant Biological process (FDR < 0.05)

gse_FPvsTR_enrich_sig_BP <- gseGO(geneList=sort(geneList_FPvsTR_enrich_sig_logFC,decreasing = T),ont = "BP", OrgDb=Candida_alb, verbose=F,pvalueCutoff = 1, nPerm = 10000) head(summary(gse_FPvsTR_enrich_sig_BP)) write.csv(gse_FPvsTR_enrich_sig_BP@result,file = "Desktop/Abhilash/study/Post-Doc/Presentations/TIM/redvsblue/results/gse_FPvsTR_enrich_sig_BP.csv")

gse_FPvsTR_enrich_sig_MF <- gseGO(geneList=sort(geneList_FPvsTR_enrich_sig_logFC,decreasing = T),ont = "MF", OrgDb=Candida_alb, verbose=F,pvalueCutoff = 1, nPerm = 10000) head(summary(gse_FPvsTR_enrich_sig_MF)) write.csv(gse_FPvsTR_enrich_sig_MF@result,file = "Desktop/Abhilash/study/Post-Doc/Presentations/TIM/redvsblue/results/gse_FPvsTR_enrich_sig_MF.csv")

gse_FPvsTR_enrich_sig_CC <- gseGO(geneList=sort(geneList_FPvsTR_enrich_sig_logFC,decreasing = T),ont = "CC", OrgDb=Candida_alb, verbose=F,pvalueCutoff = 1, nPerm = 10000) head(summary(gse_FPvsTR_enrich_sig_CC)) write.csv(gse_FPvsTR_enrich_sig_CC@result,file = "Desktop/Abhilash/study/Post-Doc/Presentations/TIM/redvsblue/results/gse_FPvsTR_enrich_sig_CC.csv")

### Visualization of Enrichment results

# ridgeplot will visualize expression distributions of core enriched genes for GSEA enriched categories. It helps users to interpret up/down-regulated pathways.

#FPvsTR #all significant (FDR < 0.05) #all GO ridgeplot(gse_FPvsTR_enrich_sig_all, fill = "p.adjust", showCategory = 10) #BP ridgeplot(gse_FPvsTR_enrich_sig_BP, fill = "p.adjust", showCategory = 10) #MF ridgeplot(gse_FPvsTR_enrich_sig_MF, fill = "p.adjust", showCategory = 10) #CC ridgeplot(gse_FPvsTR_enrich_sig_CC, fill = "p.adjust", showCategory = 10)

### GESEAplot

BiocManager::install("enrichplot") library(enrichplot)

#all significant (FDR < 0.05) gseaplot2(gse_FPvsTR_enrich_sig_all, geneSetID = 1:3, pvalue_table = TRUE, color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot") #BP gseaplot2(gse_FPvsTR_enrich_sig_BP, geneSetID = 1:3, pvalue_table = TRUE, color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot") #MF gseaplot2(gse_FPvsTR_enrich_sig_MF, geneSetID = 1:3, pvalue_table = TRUE, color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot") #CC gseaplot2(gse_FPvsTR_enrich_sig_CC, geneSetID = 1:3, pvalue_table = TRUE, color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot")

library(ggplot2) library(cowplot)

#CEC3609_3hrsvsCEC3609_6hrs pp <- lapply(1:20, function(i) {  anno <- gse_FPvsTR_enrich_sig_all[i, c("NES", "pvalue", "p.adjust")]  lab <- paste0(names(anno), "=",  round(anno, 3), collapse="\n")

 gsearank(gse_FPvsTR_enrich_sig_all, i, gse_FPvsTR_enrich_sig_all[i, 2]) + xlab(NULL) +ylab(NULL) +    annotate("text", 0, gse_FPvsTR_enrich_sig_all[i, "enrichmentScore"] * .9, label = lab, hjust=0, vjust=0) }) plot_grid(plotlist=pp, ncol=4, nrow = 5) #BP gseaplot2(gse_FPvsTR_enrich_sig_BP, geneSetID = c("GO:0007155","GO:0030447","GO:0006811","GO:0006696"), pvalue_table = T, color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot")