
rm(list = ls())    #清空数据

setwd("C:/Users/86193/Desktop/GSE198710")   #设置工作路径

#GSE198710的注释
# 清空暂存数据
rm(list=ls())
#载入原始数据
raw = read.csv("dat.csv",
                 header = T,   #第一行表头  T F
                 row.names=1,   #第一列  1 0 
                 check.names = FALSE)     
head(raw)
dim(raw)   #
# 画箱式图，比较数据分布情况
# par(mfrow=c(1,2))
boxplot(raw,col="#9400D3",outline = F) 
dat=raw
#dat = log2(raw+1)
#boxplot(dat,col="#ff4757",outline = F) 

#ids = read.csv("GSE122709symbol.csv",header=TRUE,check.names = FALSE)
#colnames(ids)=c('probe_id','symbol')  
#ids=ids[ids$symbol != '',]
#ids=ids[ids$probe_id %in%  rownames(dat),]

#dat[1:4,1:4]   
#dat=dat[ids$probe_id,] 
#ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
#ids=ids[order(ids$median,ids$symbol,decreasing = T),]
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
#ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
#dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
#rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
#dat[1:4,1:4]  #保留每个基因ID第一次出现的信息(median最大的那个)

dat <- dat[which(rowSums(dat) > 0),]   #去除一行全0


library(limma)
#install.packages('limma')


#分位数标准化
dat = normalizeBetweenArrays(dat)  
par(mfrow=c(1,1))    #调整画布
boxplot(data.frame(dat),col="#546de5",outline = F) 

write.csv(dat,'198710bzh.csv')






# 后续的分析基于标准化后的数据做
rm(list = ls())    #清空数据
dat = read.csv("198710bzh.csv",
               header = T,   #第一行表头  T F
               row.names=1,   #第一列  1 0 
               check.names = FALSE)     
#dat = dat[,c(1:6,8:15)]


#山峰图  判断正态分布
library(ggplot2)
library(ggridges)
library(reshape)
library(ggprism) 

# 清空暂存数据

df <- as.data.frame(dat)
df1 <- melt(df) 
show(df1)

ggplot(df1,aes(x=value,y=variable,fill=variable))+
  geom_density_ridges(stat = "density_ridges")+
  theme(legend.position = "none")+scale_fill_prism(palette = "candy_bright") 
####详细参数
ggplot(df1, aes(x = value, y = variable, fill = variable)) +#数据
  geom_density_ridges(stat = "density_ridges") +#绘制山脊图
  theme_prism(palette = "candy_bright",#主题样式
              base_fontface = "plain", # 字体样式
              base_family = "serif", # 字体格式
              base_size = 12,  #字体大小
              base_line_size = 0.5)+ #坐标轴粗细
  scale_fill_prism(palette = "summer")+#颜色diazo/candy_soft/candy_bright/summer/pastels/winter_soft/shades_of_gray
  theme(legend.position = "none",#图例去除
        axis.line.y = element_blank(),#去除Y轴轴线
        axis.ticks.y = element_blank())+#去除Y轴刻度
  scale_y_discrete(expand = c(0.05, 0))#调整起始图形距离X轴距离

#相关性
library(Hmisc)
res2 <- rcorr(as.matrix(df))
res2
library(corrplot)
#corrplot(res2$r, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
corrplot(res2$r, method = "number")
corrplot(res2$r, method="circle", #square方形，ellipse, 椭圆形，number数值，shade阴影，color颜色，pie饼图
         #title = "pearson",   #指定标题
         type="full",  #full完全(默认)，lower下三角，upper上三角
         #col=c("#FF6666", "white", "#0066CC"), #指定图形展示的颜色，默认以均匀的颜色展示。支持grDevices包中的调色板，也支持RColorBrewer包中调色板。
         outline = T,  #是否添加圆形、方形或椭圆形的外边框，默认为FALSE。
         diag = TRUE,  #是否展示对角线上的结果，默认为TRUE
         mar = c(0,0,0,0), #设置图形的四边间距。数字分别对应(bottom, left, top, right)。
         bg="white", #指定背景颜色
         add = FALSE, #表示是否添加到已经存在的plot中。默认FALSE生成新plot。
         is.corr = TRUE, #是否为相关系数绘图，默认为TRUE,FALSE则可将其它数字矩阵进行可视化。
         addgrid.col = "darkgray", #设置网格线颜色，当指定method参数为color或shade时， 默认的网格线颜色为白色，其它method则默认为灰色，也可以自定义颜色。
         addCoef.col = NULL, #设置相关系数值的颜色，只有当method不是number时才有效
         addCoefasPercent = FALSE, #是否将相关系数转化为百分比形式，以节省空间，默认为FALSE。
         order = "original", #指定相关系数排序的方法, 可以是original原始顺序，AOE特征向量角序，FPC第一主成分顺序，hclust层次聚类顺序，alphabet字母顺序。
         hclust.method = "complete", # 指定hclust中细分的方法，只有当指定order参数为hclust时有效。有7种可选：complete,ward,single,average,mcquitty,median,centroid。
         addrect = NULL, #是否添加矩形框，只有当指定order参数为hclust时有效， 默认不添加， 用整数指定即可添加。
         rect.col = "black", #指定矩形框的颜色。
         rect.lwd = 2, #指定矩形框的线宽。
         tl.pos = NULL,  #指定文本标签(变量名称)相对绘图区域的位置，为"lt"(左侧和顶部),"ld"(左侧和对角线),"td"(顶部和对角线),"d"(对角线),"n"(无);当type="full"时默认"lt"。当type="lower"时默认"ld"。当type="upper"时默认"td"。
         tl.cex = 1,  #设置文本标签的大小
         tl.col = "black", #设置文本标签的颜色。
         cl.pos = NULL #设置图例位置，为"r"(右边),"b"(底部),"n"(无)之一。当type="full"/"upper"时，默认"r"; 当type="lower"时，默认"b"。
         #addshade = c("negative", "positive", "all"), # 表示给增加阴影，只有当method="shade"时有效。#为"negative"(对负相关系数增加阴影135度)；"positive"(对正相关系数增加阴影45度)；"all"(对所有相关系数增加阴影)。
         #shade.lwd = 1,  #指定阴影线宽。
         #shade.col = "white",  #指定阴影线的颜色。    
         #p.mat= res1$p,sig.level= 0.01,insig= "pch", pch.col= "blue", pch.cex= 3,#只有指定矩阵的P值，sig.level，pch等参数才有效。只有当insig = "pch"时，pch.col和pch.cex参数才有效
)

# 分位数PCA  主成分分析  放文章
library(factoextra)
library(FactoMineR)
qlog = t(df)
type = c(rep("c",5),rep("s",5))     #改
type <- as.data.frame(type)
qlogtype <- data.frame(qlog,type)

qlogtype.pca <- PCA(qlogtype[,-35182],graph = FALSE)    #改
qlogtype$type <- factor(qlogtype$type)

fviz_pca_ind(qlogtype.pca,
             habillage = qlogtype$type, # color by groups
             palette = c("#546de5","#ff4757","green"),
             addEllipses = FALSE, # Concentration ellipses
             ellipse.level = 0.95, # Specify the confidence level
             repel = FALSE,
             invisible = 'quali'
) + 
  ggforce::geom_mark_ellipse(aes(fill = Groups, color = Groups)) +
  theme(legend.position = 'bottom')



write.csv(dat,"198710_.csv")






#到这

rm(list=ls())
df = read.csv("198710_.csv",header=TRUE,row.names=1,check.names = FALSE)  
########################
dfs <- df[,c(1:10)]

#差异分析
library('limma')
#分组信息
type = c(rep("s",5),rep("c",5))   #改

#####
group <- as.data.frame(type)
#定义比较组，按分组数目修改
group <- group[,1] 
#把group设置成一个model matrix
design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(dfs)

#计算
fit <- lmFit(dfs,design)
cont.matrix<-makeContrasts(paste0(unique(group),collapse = "-"),levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
tempOutput = topTable(fit2,coef=1,n=Inf,adjust="BH") 
nrDEG = na.omit(tempOutput) 
#####

#阈值
logfoldchange = 1
padj = 0.05
#输出
allDiff <- nrDEG
diff=allDiff
write.csv(diff, "limmaOut.csv")   #输出全部结果
######



#火山图
library(ggplot2)#声明画图包
library(ggrepel)
# 清空暂存数据
rm(list=ls())

dataset = read.csv("limmaOut.csv",header = T,check.names = FALSE)

cut_off_pvalue = 0.05
cut_off_logFC = 0.5

# 根据阈值分别为上调基因设置‘up’，下调基因设置‘Down’，无差异设置‘Stable’，保存到change列
# 这里的change列用来设置火山图点的颜色
dataset$change = ifelse(dataset$P.Value < cut_off_pvalue & abs(dataset$logFC) >= cut_off_logFC, 
                        ifelse(dataset$logFC> cut_off_logFC ,'Up','Down'),'Stable')
# 绘制火山图====================================
p <- ggplot(
  #设置数据
  dataset, 
  aes(x = logFC, 
      y = -log10(P.Value), 
      colour=change)) +
  geom_point(alpha = 0.75 , size = 3) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  
  # 辅助线
  geom_vline(xintercept=c(-cut_off_logFC,cut_off_logFC),lty=4,col="black",lwd=0.5) +
  geom_hline(yintercept = -log10(cut_off_pvalue),lty=4,col="black",lwd=0.5) +
  
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  
  # 图例
  theme(plot.title = element_text(hjust = 1), 
        legend.position="right", 
        legend.title = element_blank()
  )

p

#标注
###########################################################
dataset$label = ifelse(dataset$P.Value < cut_off_pvalue & abs(dataset$logFC) >= 4, as.character(dataset[,1]),"")
p + geom_text_repel(data = dataset, aes(x = logFC, 
                                        y = -log10(P.Value), 
                                        label = label),
                    size = 3,box.padding = unit(0.5, "lines"),
                    point.padding = unit(0.8, "lines"), 
                    segment.color = "black", 
                    show.legend = FALSE)
###########################################################


#热图
# 清空暂存数据
rm(list=ls())
library(pheatmap)

df = read.csv("198710_.csv",header=TRUE,row.names=1,check.names = FALSE) 

interest <- read.table('retu.txt')
dataset <- df[which(rownames(df)%in%interest[,1]),]
ds <- dataset



exp_ds = ds  #热图的行\列数量
cell_list = c(rep("s",5),rep("c",5))  
annotation_c <- data.frame(cell_list)
rownames(annotation_c) <- colnames(exp_ds)

pheatmap((exp_ds), #表达数据    #可以改
         cluster_rows = F,#行聚类
         cluster_cols = F,#列聚类
         annotation_col = annotation_c, #样本分类数据
         annotation_legend = T, # 显示样本分类
         show_rownames = T,# 显示行名
         show_colnames = T,# 显示列名
         scale = "row", #对行标准化    #可以改
         #legend_breaks = c(0:6),  #图例范围
         border = FALSE,   #边框线
         display_numbers = F,   #显示表达量
         color = colorRampPalette(c("#546de5", "#d2dae2","#ff4757"))(100) # 热图基准颜色
)

