#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma")

#install.packages("pheatmap")


#引用包
library(limma)
library(pheatmap)

logFCfilter=1        #logFC过滤阈值
adjPfilter=0.05      #矫正后p值阈值
gene="GADD45B"       #基因的名称(根据研究的目标基因进行修改)
expFile="merge.normalzie.txt"      #表达数据文件
setwd("C:\\Users\\lexb\\Desktop\\geoGene\\19.groupDiff")      #设置工作目录

#读取输入文件
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
rt=avereps(rt)

#去除对照组的样品
Type=gsub("(.*)\\_(.*)\\_(.*)", "\\3", colnames(rt))
rt=rt[,Type=="Treat"]

#按照目标基因表达的中位值对样品分组
low=rt[gene,]<=median(rt[gene,])
high=rt[gene,]>median(rt[gene,])
lowRT=rt[,low]
highRT=rt[,high]
conNum=ncol(lowRT)        #低表达组样品数目
treatNum=ncol(highRT)     #高表达组样品数目
rt=cbind(lowRT,highRT)

#差异分析
Type=c(rep("con",conNum),rep("treat",treatNum))
design <- model.matrix(~0+factor(Type))
colnames(design) <- c("con","treat")
fit <- lmFit(rt,design)
cont.matrix<-makeContrasts(treat-con,levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
allDiff=topTable(fit2,adjust='fdr',number=200000)

#输出所有基因的差异情况
write.table(allDiff,file="all.txt",sep="\t",quote=F)

#输出显著的差异基因
diffSig=allDiff[with(allDiff, (abs(logFC)>logFCfilter & adj.P.Val < adjPfilter )), ]
diffSigOut=rbind(id=colnames(diffSig),diffSig)
write.table(diffSigOut,file="diff.txt",sep="\t",quote=F,col.names=F)

#绘制差异基因热图
geneNum=20      #设置展示基因的数据
diffUp=diffSig[diffSig$logFC>0,]
diffDown=diffSig[diffSig$logFC<0,]
geneUp=row.names(diffUp)
geneDown=row.names(diffDown)
if(nrow(diffUp)>geneNum){geneUp=row.names(diffUp)[1:geneNum]}
if(nrow(diffDown)>geneNum){geneDown=row.names(diffDown)[1:geneNum]}
hmExp=rt[c(geneUp,geneDown),]
#注释文件
Type=c(rep("Low",conNum),rep("High",treatNum))
Type=factor(Type, levels=c("Low", "High"))
names(Type)=colnames(rt)
Type=as.data.frame(Type)
colnames(Type)[1]=gene
#输出图形
pdf(file="heatmap.pdf", width=7.5, height=5)
pheatmap(hmExp, 
         annotation=Type, 
         color = colorRampPalette(c("blue", "white", "red"))(50),
         cluster_cols =F,
         show_colnames = F,
         scale="row",
         fontsize = 8,
         fontsize_row=6,
         fontsize_col=8)
dev.off()

#输出相关性分析的数据
hmExp=rt[unique(c(gene,geneUp,geneDown)),]
hmExpOut=rbind(id=colnames(hmExp),hmExp)
write.table(hmExpOut,file="corInput.txt",sep="\t",quote=F,col.names=F)

#绘制火山图
pdf(file="vol.pdf",width=5,height=5)
xMax=4      #x轴范围
yMax=10     #y轴范围
plot(allDiff$logFC, -log10(allDiff$adj.P.Val), ylab="-log10(adj.P.Val)",xlab="logFC", col="grey",
     main="Volcano", ylim=c(0,yMax),xlim=c(-xMax,xMax),yaxs="i",pch=20, cex=1)
diffSub=subset(allDiff, adj.P.Val<adjPfilter & logFC>logFCfilter)
points(diffSub$logFC, -log10(diffSub$adj.P.Val), pch=20, col="red",cex=1.2)
diffSub=subset(allDiff, adj.P.Val<adjPfilter & logFC<(-logFCfilter))
points(diffSub$logFC, -log10(diffSub$adj.P.Val), pch=20, col="green",cex=1.2)
abline(v=0,lty=2,lwd=3)
dev.off()


######生信自学网: https://www.biowolf.cn/
######课程链接1: https://shop119322454.taobao.com
######课程链接2: https://ke.biowolf.cn
######课程链接3: https://ke.biowolf.cn/mobile
######光俊老师邮箱: seqbio@foxmail.com
######光俊老师微信: seqBio

