#### Example rm(list = ls()) library(ggplot2) library(dplyr) set.seed(2012) data <- data.frame( id = paste0("ca",1:20), p = runif(20,0,0.05), # runif(n, min, max): n = number of observations, min/max = distribution range change = sample(c("up","down"),20,replace = T) ) head(data) data <- data %>% mutate(p2 = ifelse(change == "up", -log10(p), log10(p))) %>% arrange(change,p2) data$id = factor(data$id,levels = unique(data$id),ordered = T) head(data) tmp = with(data, labeling::extended(range(p2)[1], range(p2)[2], m = 5));tmp lm = tmp[c(1,length(tmp))];lm ggplot(data, aes(x=id, y=p2)) + geom_bar(stat='identity', aes(fill=change), width=.7)+ # fill = 'pink' # change fill color; width = 0.9 # change bar width coord_flip()+ # flip cartesian coordinates to make horizontal bars vertical theme_light() + # See: https://zhuanlan.zhihu.com/p/463041897 # theme_grey() default background with light gray and white grid lines, no border # theme_bw() similar to default but with white background and light gray grid lines, no border # theme_linedraw() white background with black grid lines and borders # theme_light() white background with light gray grid lines and light gray borders ylim(lm)+ scale_y_continuous(breaks = tmp, labels = abs(tmp))+ theme( legend.position = "none", # legend.position parameter places legend in chart area panel.border = element_blank() # remove border ) #### Practical Application rm(list = ls()) library(ggplot2) library(dplyr) input <- read_tsv("F:/Decoding Helix Courses/GEOdatasets_spsis1/DEG_SPSIS/RRP/all_signi_result_P_F.tsv") head(input) input_up <- input[1:95,] input_down <- input[96:164,] gene.df_up <- bitr(input_up$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", # drop = TRUE, # remove NA and unmatched gene names OrgDb = org.Hs.eg.db) gene.df_down <- bitr(input_down$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", # drop = TRUE, # remove NA and unmatched gene names OrgDb = org.Hs.eg.db) ### GO Enrichment Analysis ---- ego_up <- enrichGO(gene = gene.df_up$ENTREZID, OrgDb = org.Hs.eg.db, ont = "ALL", # default is "MF", ALL means all categories pvalueCutoff =0.05, qvalueCutoff =0.05, readable = TRUE) %>% mutate(change='up') ego_up_data <- as.data.frame(ego_up) ego_down <- enrichGO(gene = gene.df_down$ENTREZID, OrgDb = org.Hs.eg.db, ont = "ALL", # default is "MF", ALL means all categories pvalueCutoff =0.05, qvalueCutoff =0.05, readable = TRUE) %>% mutate(change='down') ego_down_data <- as.data.frame(ego_down) ego_merge_data <- bind_rows(ego_up_data,ego_down_data) ego_merge_data1 <- ego_merge_data[c(1:3,7,11)] table(ego_merge_data1$ONTOLOGY) ## Extract GO results for plotting ego_merge_BP <- ego_merge_data1[ego_merge_data1$ONTOLOGY =="BP",] %>% mutate(p2 = ifelse(change == "up", -log10(p.adjust), log10(p.adjust))) %>% arrange(change,p2) ego_merge_CC <- ego_merge_data1[ego_merge_data1$ONTOLOGY =="CC",] %>% mutate(p2 = ifelse(change == "up", -log10(p.adjust), log10(p.adjust))) %>% arrange(change,p2) ego_merge_MF <- ego_merge_data1[ego_merge_data1$ONTOLOGY =="MF",] %>% mutate(p2 = ifelse(change == "up", -log10(p.adjust), log10(p.adjust))) %>% arrange(change,p2) ## BP Bidirectional Histogram ---- # Determine BP plot scale ego_merge_BP$Description = factor(ego_merge_BP$Description, levels = unique(ego_merge_BP$Description), ordered = T) head(ego_merge_BP) tmp = with(ego_merge_BP, labeling::extended(range(p2)[1], range(p2)[2], m = 10));tmp lm = tmp[c(1,length(tmp))];lm # Determine BP plot data data <- ego_merge_BP[c(1:5,(nrow(ego_merge_BP)-4):nrow(ego_merge_BP)),] id <- data$Description p2 <- data$p2 # BP Plot pdf("./GO_KEGG/Bidirectional histogram/GO_BP_Bi_barplot.pdf") ggplot(data, aes(x=id, y=p2)) + geom_bar(stat='identity', aes(fill=change), width=.7)+ coord_flip()+ theme_light() + ylim(lm)+ scale_y_continuous(breaks = tmp, labels = abs(tmp))+ theme( legend.position = "right" #, panel.border = element_blank() # "none", "left", "right", "bottom", "top": # "left" left, "right" right, "bottom" bottom, "top" top # panel.border = element_blank() removes the outer border (default has border) )+ scale_fill_manual(values=c("#2874C5", "#f87669"))+ facet_grid(ONTOLOGY~., scale='free') # Add ONTOLOGY type dev.off() ## CC Bidirectional Histogram ---- # Determine CC plot scale ego_merge_CC$Description = factor(ego_merge_CC$Description, levels = unique(ego_merge_CC$Description), ordered = T) head(ego_merge_CC) tmp = with(ego_merge_CC, labeling::extended(range(p2)[1], range(p2)[2], m = 12));tmp lm = tmp[c(1,length(tmp))];lm # Determine CC plot data data <- ego_merge_CC[c(1:5,(nrow(ego_merge_CC)-4):nrow(ego_merge_CC)),] id <- data$Description p2 <- data$p2 # CC Plot pdf("./GO_KEGG/Bidirectional histogram/GO_CC_Bi_barplot.pdf") ggplot(data, aes(x=id, y=p2)) + geom_bar(stat='identity', aes(fill=change), width=.7)+ coord_flip()+ theme_light() + ylim(lm)+ scale_y_continuous(breaks = tmp, labels = abs(tmp))+ theme( legend.position = "right" #, panel.border = element_blank() # "none", "left", "right", "bottom", "top": # "left" left, "right" right, "bottom" bottom, "top" top # panel.border = element_blank() removes the outer border (default has border) )+ scale_fill_manual(values=c("#2874C5", "#f87669"))+ facet_grid(ONTOLOGY~., scale='free') # Add ONTOLOGY type dev.off() ## MF Bidirectional Histogram ---- # Determine MF plot scale ego_merge_MF$Description = factor(ego_merge_MF$Description, levels = unique(ego_merge_MF$Description), ordered = T) head(ego_merge_MF) tmp = with(ego_merge_MF, labeling::extended(range(p2)[1], range(p2)[2], m = 10));tmp lm = tmp[c(1,length(tmp))];lm # Determine MF plot data data <- ego_merge_MF[c(1:5,(nrow(ego_merge_MF)-4):nrow(ego_merge_MF)),] id <- data$Description p2 <- data$p2 # MF Plot pdf("./GO_KEGG/Bidirectional histogram/GO_MF_Bi_barplot.pdf") ggplot(data, aes(x=id, y=p2)) + geom_bar(stat='identity', aes(fill=change), width=.7)+ coord_flip()+ theme_light() + ylim(lm)+ scale_y_continuous(breaks = tmp, labels = abs(tmp))+ theme( legend.position = "right" #, panel.border = element_blank() # "none", "left", "right", "bottom", "top": # "left" left, "right" right, "bottom" bottom, "top" top # panel.border = element_blank() removes the outer border (default has border) )+ scale_fill_manual(values=c("#2874C5", "#f87669"))+ facet_grid(ONTOLOGY~., scale='free') # Add ONTOLOGY type dev.off() ### KEGG Enrichment Analysis ---- kk_up <- enrichKEGG(gene = gene.df_up$ENTREZID, organism = "hsa", # keyType = "kegg", pvalueCutoff =0.9999, qvalueCutoff =0.9999) %>% mutate(change='up') kk_up_data <- as.data.frame(kk_up) kk_down <- enrichKEGG(gene = gene.df_down$ENTREZID, organism = "hsa", # keyType = "kegg", pvalueCutoff =0.9999, qvalueCutoff =0.9999) %>% mutate(change='down') kk_down_data <- as.data.frame(kk_down) kk_merge_data <- bind_rows(kk_up_data,kk_down_data) kk_merge_data1 <- kk_merge_data[c(1,2,6,10)] kk_merge <- kk_merge_data1 %>% mutate(p2 = ifelse(change == "up", -log10(p.adjust), log10(p.adjust))) %>% arrange(change,p2) ## KEGG Bidirectional Histogram ---- # Determine KEGG plot scale kk_merge$Description = factor(kk_merge$Description, levels = unique(kk_merge$Description), ordered = T) head(kk_merge) tmp = with(kk_merge, labeling::extended(range(p2)[1], range(p2)[2], m = 13));tmp lm = tmp[c(1,length(tmp))];lm # Determine KEGG plot data data <- kk_merge[c(1:5,(nrow(kk_merge)-4):nrow(kk_merge)),] id <- data$Description p2 <- data$p2 # KEGG Plot pdf("./GO_KEGG/Bidirectional histogram/kegg_Bi_barplot.pdf") ggplot(data, aes(x=id, y=p2)) + geom_bar(stat='identity', aes(fill=change), width=.7)+ coord_flip()+ theme_light() + ylim(lm)+ scale_y_continuous(breaks = tmp, labels = abs(tmp))+ theme( legend.position = "right" #, panel.border = element_blank() # "none", "left", "right", "bottom", "top": # "left" left, "right" right, "bottom" bottom, "top" top # panel.border = element_blank() removes the outer border (default has border) )+ scale_fill_manual(values=c("#2874C5", "#f87669")) dev.off() ########################################## # miRNA-mRNA and mRNA-miRNA sections ########################################## # Read original mRNA-miRNA interaction data mir_mRNA=read.table("starbase/mRNA_miRNA_interaction.txt", sep="\t", stringsAsFactors=F, header=T) # Extract only miRNAname and geneName columns and remove duplicate rows mir_mRNA=unique(mir_mRNA[,c("miRNAname","geneName")]) # View first six rows head(mir_mRNA) # Convert dataframe to list with gene names as list names mRNA_mir_list=unstack(mir_mRNA,miRNAname~geneName) # Convert dataframe to list with miRNA names as list names mir_mRNA_list=unstack(mir_mRNA,geneName~miRNAname) # Set random seed for reproducibility set.seed(123) # Randomly select 100 miRNAs to simulate candidate miRNAs mir_candidate1=sample(names(mir_mRNA_list),100) # Get subset of the list to get miRNA-mRNA interactions, then convert to dataframe mir_target1=stack(mir_mRNA_list[mir_candidate1]) # Rename dataframe columns names(mir_target1)=c("gene","miRNA") # Randomly select 100 mRNAs to simulate candidate mRNAs mRNA_candidate=sample(names(mRNA_mir_list),100) # Get subset of the list to get mRNA-miRNA interactions, then convert to dataframe mRNA_target=stack(mRNA_mir_list[mRNA_candidate]) # Rename dataframe columns names(mRNA_target)=c("miRNA","gene") ########################################## # miRNA-lncRNA and lncRNA-miRNA sections ########################################## # Read original lncRNA-miRNA interaction data mir_lnc=read.table("starbase/lncRNA_miRNA_interaction.txt", sep="\t", stringsAsFactors=F, header=T) # Extract only miRNAname and geneName columns and remove duplicate rows mir_lnc=unique(mir_lnc[,c("miRNAname","geneName")]) # View first six rows head(mir_lnc) # Convert dataframe to list with lncRNA names as list names lnc_mir_list=unstack(mir_lnc,miRNAname~geneName) # Convert dataframe to list with miRNA names as list names mir_lnc_list=unstack(mir_lnc,geneName~miRNAname) # Set random seed for reproducibility set.seed(123) # Randomly select 100 miRNAs to simulate candidate miRNAs mir_candidate2=sample(names(mir_lnc_list),100) # Get subset of the list to get miRNA-lncRNA interactions, then convert to dataframe mir_target2=stack(mir_lnc_list[mir_candidate2]) # Rename dataframe columns names(mir_target2)=c("lncRNA","miRNA") # Randomly select 100 lncRNAs to simulate candidate lncRNAs lnc_candidate=sample(names(lnc_mir_list),100) # Get subset of the list to get lncRNA-miRNA interactions, then convert to dataframe lnc_target=stack(lnc_mir_list[lnc_candidate]) # Rename dataframe columns names(lnc_target)=c("miRNA","lncRNA")