rm(list = ls())
library(openxlsx)
library(dplyr)
##
## 载入程辑包：'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggdendro)
library(ggplot2)
###############################input data
dir_path <- "C:\\Users\\liyix\\OneDrive\\Desktop\\"
dir_path_name <- list.files(pattern = ".*xlsx",dir_path,full.names = T, recursive = F)
#dir_path_name
#########################################################sig_stress_related_genes
#getSheetNames(grep("data.xlsx",dir_path_name, value = T))
data_1 <- read.xlsx(grep("data.xlsx",dir_path_name, value = T), sheet = 3,startRow = 1)
#dim(data_1) #[1] 532 31
l <- LETTERS
letter_1 <- unlist(lapply(2:5, function(n) combn(l, n, paste, collapse="")))
#length(letter_1)
data_1$Accession <- letter_1[1:nrow(data_1)]
colnames(data_1)[29] <- "pvalue"
#min(data_1$pvalue) #3.842396e-16
data_1 <- data_1[data_1$pvalue > min(data_1$pvalue), ]
#dim(data_1) #[1] 383 31
##########################################################cutoff = 2
colnames(data_1)[28] <- "fold"
data_1 <- data_1[data_1$fold >= 2.2 |
data_1$fold <= 1/2.2,] ##cutoff = 2
##########################################################cutoff = 2
#dim(data_1) #[1] 249 31
#View(data_1)
s1 <- data_1[grep("mitochondrion",data_1$Cellular.Component, ignore.case=TRUE, value = F), ]
#View(s1)
s2 <- s1[, c(2, 28, 29)]
#dim(s2) #[1] 12 3
#head(s2)
s2_up <- s2[s2$fold >= 1.5, ]
s2_dow <- s2[s2$fold <= 1.5, ]
#head(s2_up)
#dim(s2_up); dim(s2_dow) #[1] 57 3 [1] 10 3
s3 <- rbind(s2_up, s2_dow)
#head(s3)
#write.csv(s3, paste0(dir_path,Sys.Date(),"-","mitochondrion_15_fold_list.csv"),row.names = FALSE,na = "")
######################################################all data
data_2 <- data_1
#dim(data_2) #[1] 7047 31
data_2 <- data_2[!is.na(data_2$pvalue),]
#dim(data_2) #[1] 6828 31
#colnames(data_2)
data_2 <- data_2[data_2$pvalue > min(data_2$pvalue), ]
#################################################cutoff 2__-reminder
data_2 <- data_2[data_2$fold >= 2.2 |
data_2$fold <= 1/2.2,] ##cutoff 2
#####################################################
data_2 <- data_2[, c(2, 30, 31)]
#head(data_2)
colnames(data_2)<- c("Accession", "TR", "Con")
#dim(data_2) #[1] 49 3
#head(data_2)
rownames(data_2) <- data_2$Accession
data_2$Accession <- NULL
#dim(data_2) #[1] 6679 2
#View(data_2)
###########################################################dendro_1_LEFT
#data_2[is.na(data_2)] <- 0
dd <- dist(scale(data_2), method = "euclidean")
hc <- hclust(dd, method = "complete")
#ggdendrogram(hc, rotate = T, size = 0.1)
######
dend <- as.dendrogram(hc)
dendata <- dendro_data(dend,type = "rectangle")
#dim(data_2)#[1] 6679 2
#dim(dendata)
#add column family to labs dataframe (look like VLOOKUP in excel)
labs <- label(dendata)
####plot
data_plot <- segment(dendata)
#dim(data_plot)
#head(data_plot)
write.csv(data_plot, paste0(dir_path,Sys.Date(),"-","segment_dendrogram.csv"),row.names = FALSE,na = "")
write.csv(labs, paste0(dir_path,Sys.Date(),"-","labs_dendrogram.csv"),row.names = FALSE,na = "")
######dendro_1
labs$y <- -0.1
#head(labs)
head(data_plot)
## x y xend yend
## 1 19.182617 4.64154854 5.693359 4.64154854
## 2 5.693359 4.64154854 5.693359 1.10158625
## 3 5.693359 1.10158625 1.500000 1.10158625
## 4 1.500000 1.10158625 1.500000 0.05656794
## 5 1.500000 0.05656794 1.000000 0.05656794
## 6 1.000000 0.05656794 1.000000 0.00000000
#dim(data_2)
p <- ggplot(data_plot) +
geom_segment(aes(x=x, y=y, xend=xend, yend=yend, color=""),size = .1) +
coord_flip() +
#geom_text(data=labs,size= 1,
# aes(label=label, x=x, y= y),vjust = 0.5,hjust = 0, color = "black")+
geom_segment(data=labs,aes(x=x, y=y, xend=x, yend=0, color=""),size = .1) +
scale_colour_manual(values = c("black")) +
scale_x_continuous(expand = c(0,0),limits = c(0.5,nrow(data_2)+0.5)) +
scale_y_reverse(expand = c(0, 0))
p_dendro <- p + theme(legend.position = "right")+
theme(axis.line = element_blank(), plot.background = element_rect(fill="green"))
p_dendro
##################################################2___HM-1
data_s3 <- s3
#head(data_s3)
data_s3$label[data_s3$fold > 1.5] <- "red"
data_s3$label[data_s3$fold< 1.5] <- "green"
data_s3 <- data_s3[, c(1,4)]
#head(data_s3)
data_label <- labs
colnames(data_label)[3] <- "Accession"
#head(data_label); head(data_s3)
data_12 <- merge(data_label, data_s3, by = "Accession", all.x = T, sort = F)
#dim(data_12) #[1] 49 4
data_12$label[is.na(data_12$label)] <- "gray"
#table(data_12$label)
#head(data_12)
data_12$xend <- data_12$x
data_12$yend <- 2
data_12$y <- 1
#head(data_plot)
data_4 <- merge(data_plot[, c("x", "xend")], data_12, by = c("x", "xend"), all = F)
#dim(data_4) #[1] 49 8
#head(data_4)
######plot
data_4$label <- factor(data_4$label, levels = c("green","red","gray"))
p_h1 <- ggplot(data_4) +
geom_segment(aes(x=x, y=y, xend=xend, yend=yend, color = label),size = .1) +
coord_flip() +
#geom_text(data=data_4,size= 1,
# aes(label=Accession, x=x, y= y),vjust = 0.5,hjust = 0, color = "black")+
#geom_segment(data=labs,aes(x=x, y=y, xend=x, yend=0, color=""),size = .1) +
scale_colour_manual(values = c("#1a53ff","darkred","#cce5ff"),
labels = c("DM","UM","Other"),
name = "Category") +
scale_x_continuous(expand = c(0,0),limits = c(0.5,nrow(data_2)+0.5)) +
scale_y_reverse(expand = c(0, 0)) +
theme(legend.position = "right")+
theme(axis.line =element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title =element_blank(),
panel.grid=element_blank(),
plot.background = element_rect(fill="green"),
legend.key = element_rect(colour = NA, fill = NA)) +
labs(x = "", y = "", title = "")
p_h1
#################################################################3____heatmap_2
data_21 <- data_2
#View(data_21)
#head(data_21)
data_21$num <- row.names(data_21)
data_hm <- gather(data_21, key = key, value = value, -num)
#head(data_hm)
data_hm$num <- factor(data_hm$num, levels = unique(labs$label))
data_hm$inver <- cut(data_hm$value, breaks = seq(0, 200, 50), include.lowest = T)
#table(data_hm$inver)
#str(data_hm)
#data_hm$A
#unique(data_hm$key)
library(RColorBrewer)
#brewer.pal.info
head(data_hm)
## num key value inver
## 1 CM TR 159.5 (150,200]
## 2 CN TR 158.8 (150,200]
## 3 CO TR 157.9 (150,200]
## 4 CP TR 154.9 (150,200]
## 5 CQ TR 150.1 (150,200]
## 6 CR TR 149.6 (100,150]
p2 <- ggplot() +
geom_tile(data = data_hm,aes(x = key,
y= num, height = 1,
fill = inver,width = 1),size= 0.2) +
ylab("") + xlab("") +
theme(legend.position = "right",
plot.background = element_rect(fill="green")) +
scale_y_discrete(expand = c(0, 0),position = "right") +
scale_x_discrete(expand = c(0, 0)) +
scale_fill_manual(values = brewer.pal(n = length(unique(data_hm$inver)), name = "Pastel1"),
name = "Abundance")
#scale_fill_gradient2(low="blue", high="red", name = expression('Value'), na.value="gray90")
p_h2 <- p2 + guides(colour = guide_legend("Missing value", override.aes = list(fill="gray90"))) +
theme(axis.text = element_text(angle = 0, vjust = 0.5, hjust = 0.5, color = "black", size = 10)) +
scale_colour_manual(values= "gray80")
#theme(axis.ticks = element_blank(),
#axis.text.y = element_blank(),
#panel.background = element_blank(),
# legend.key.height = unit(0.8,"line"),
#legend.key.width = unit(0.7,"line"))
p_h2
###################################################add label _4
head(data_4)
## x xend Accession y label yend
## 1 1 1 AFP 1 gray 2
## 2 10 10 AFX 1 gray 2
## 3 11 11 AFM 1 gray 2
## 4 12 12 AFR 1 gray 2
## 5 13 13 AFY 1 gray 2
## 6 14 14 AFU 1 gray 2
library(ggrepel)
label_1 <- data_4[data_4$label != "gray", ]$Accession
p_label <- ggplot(data_4,
aes(x=0, y= x, label = ifelse(Accession %in% label_1, Accession, ""))) +
geom_point(color = "red") +
geom_text_repel(
point.padding = 0.001,
size= 3,
nudge_x = .15,
force_pull = 0,
nudge_y = .5,
segment.curvature = -1e-20,
segment.size = 0.1,
arrow = arrow(length = unit(0, "npc"))
) + theme(legend.position = "right") +
scale_y_continuous(expand = c(0, 0),limits = c(0.5,nrow(data_2)+0.5)) +
scale_x_continuous(expand = c(0, 0), limits = c(0,0.25))
p_label
##################################################6___MERGE
library(patchwork)
p_dendro + p_h1 + p_h2 + p_label +
plot_layout(widths = c(1,0.5,2,0.5), guides = 'collect')
#################################################################
(p_dendro + theme(axis.text = element_blank(),
legend.position= "",
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.ticks.y = element_line(color= NA),
panel.background=element_rect(fill=NA),
plot.background = element_rect(fill=NA),
panel.grid=element_blank(),
plot.margin = unit(c(0.2,0,0,0), "cm"))) +
(p_h1 + theme(axis.text = element_blank(),
#legend.position=c(1,0.8),
#legend.direction = "vertical",
axis.title = element_blank(),
axis.ticks = element_blank(),
panel.background=element_rect(fill=NA),
plot.background = element_rect(fill=NA),
panel.grid=element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.length.y = unit(0, "mm"),
plot.margin = unit(c(0.2,0,0,0), "cm"))) +
theme(legend.position=c(30,0.92), legend.direction = "vertical") +
(p_h2 + theme(axis.text.y = element_blank(),
#legend.position=c(1.2,0.5),
#legend.direction = "vertical",
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.length.y = unit(0, "mm"),
panel.background=element_rect(fill=NA),
plot.background = element_rect(fill=NA),
plot.margin = unit(c(0.2,0,0,0), "cm"))) +
theme(legend.position=c(1.5,0.6), legend.direction = "vertical") +
(p_label + theme(panel.background=element_rect(fill=NA),
plot.background = element_rect(fill=NA),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
panel.grid=element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.length.y = unit(0, "mm"),
plot.margin = unit(c(0.2,3,0,0), "cm"))) +
plot_layout(widths = c(1,0.1,2,0.6), guides = 'auto')
####################################################
ggsave(filename = paste0(Sys.Date(),"heatmap.tif"),
plot = last_plot(), device = "tiff", path = dir_path,
width = 15, height = 18, units = "cm",
dpi = 300, limitsize = TRUE)
ggsave(filename = paste0(Sys.Date(),"heatmap.pdf"),
plot = last_plot(), device = "pdf", path = dir_path,
width = 15, height = 18, units = "cm",
dpi = 300, limitsize = TRUE)
## Warning: 'mode(bg)' 新的值和旧的值不一样
## ==> NOT changing 'bg'