rm(list = ls())  
options(stringsAsFactors = F)


#################差异分析----
library(readxl)
data <- read_excel("Data.exp.xls")
colnames(data)<-c("Tumor.1","Tumor.2","Tumor.3","Tumor.4","Tumor.5",
                  "Normal.1","Normal.2","Normal.3","Normal.4")

#因子排序
group <- c(rep("treat",5),rep("con",4)) 
group <- factor(group,levels = c("treat","con"),ordered = F)


library(limma)
library(edgeR)

#分组矩阵design构建
design <- model.matrix(~0+factor(group)) #构建分组矩阵
colnames(design) <- levels(factor(group))
rownames(design) <- colnames(data)
contrast.matrix<-makeContrasts(paste0(unique(group),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较

##step1
fit <- lmFit(data,design)#fit linear model
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)

write.csv(nrDEG,"22.06.09.cir.limma.raw.csv")

##筛选
diff <- data %>% filter(P.Value <0.01 & abs(logFC)>2)

##################热图----
##########画图
library(pheatmap)
annotation_col = data.frame(type = factor(rep(c("Tumor","Normal"),c(5,4))))
rownames(annotation_col) = colnames(df)

ann_colors = list(type = c(Normal = "blue", Tumor = "red"))

pheatmap(df,cellwidth =16,
         cellheight = 0.2,
         fontsize = 8,
         method="spearman", #计算gene或sample之间的相关性的方法，可选"pearson" (default), "kendall", or "spearman"
         scale="row", #为基因做scale
         cluster_rows=T,#为基因做聚类FALSE
         cluster_cols=T,#为sample做聚类
         color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
         show_colnames=F,show_rownames =F,
         annotation_col = annotation_col,
         annotation_colors = ann_colors,
         treeheight_row = "0",treeheight_col = "0",#不画树
         border_color = "NA",
         filename = "heatmap.diff.pdf")

############火山图
library(ggplot2)
library(ggrepel)


data33 <- mutate(data2, log10pvalue = -log10(P.Value), 
                signif. = (ifelse(P.Value > 0.01, "No.signif.",
                                  ifelse(logFC > 2, "Up.Reg.",
                                         ifelse(logFC < -2, "Down.Reg.", "No.signif.")))))

p=ggplot(data33, aes(x = logFC, y = log10pvalue, 
                    color = factor(signif., levels = c("Up.Reg.", "Down.Reg.", "No.signif."))))+
  geom_point(size = 2, alpha = 0.7)+theme_bw()+
  labs(color = "Signifance")+
  ylab("-log10(PValue)")+
  xlab("log2(FC)")+
  scale_color_manual(values = c("#e6550d", "#3182bd", "gray60"))+
  # scale_color_d3(palette = "category10", )+
  theme(axis.text.x = element_text(size = 15, face = "plain", hjust=0.5, colour="black", family = "ArialMT"),
        axis.text.y = element_text(size = 15, face = "plain", hjust=1, colour="black", family = "ArialMT"),
        axis.title = element_text(size = 15, face = "plain", colour="black", family = "ArialMT"),
        panel.grid = element_blank(),
        legend.position = "bottom", 
        legend.title = element_text(size = 12, face = "plain", colour="black", family = "ArialMT"),
        legend.text = element_text(size = 12, face = "plain", colour="black", family = "ArialMT"))+
  geom_hline(aes(yintercept = 2), linetype = "dashed", color = "#d62728")+
  geom_vline(aes(xintercept = 2), linetype = "dashed", color = "#d62728")+
  geom_vline(aes(xintercept = -2), linetype = "dashed", color = "#d62728")+ 
  ggtitle("Up.Reg=728 Down.Reg=228")+ 
  scale_x_continuous(breaks=seq(-6, 6, 2))+   ## X 轴每隔单位显示一个刻度+
  scale_y_continuous(breaks=seq(0, 8, 2))+
theme(plot.title = element_text(hjust = 0.45))  #标题居中


table(data33$signif.)

p

ggsave("vocanol.pdf",width = 7.09, height =6,dpi = 300)   #保存成pdf

heatmap.png

vol.png