> library(data.table)
> library(tidyverse)
> library(ggsignif)
> library(RColorBrewer)
#install.packages("gplots")
> library(gplots)
> library(limma)
> library(ggplot2)
> library(ggrepel)
> library(Rcpp)
> rt=read.table("diffGeneExp.txt",sep="\t",header=T,check.names=F)
> rt=as.matrix(rt)
id logFC logCPM PValue adj.P.Val CKM "CKM" "-8.576279957" " 5.459229724" " 0.000000e+00" " 0.000000e+00" ACTA1 "ACTA1" "-7.228533360" " 6.535926507" " 0.000000e+00" " 0.000000e+00" MYLPF "MYLPF" "-7.209575984" " 2.618654657" " 0.000000e+00" " 0.000000e+00" PYGM "PYGM" "-7.202706200" " 4.074830196" " 0.000000e+00" " 0.000000e+00" SLN "SLN" "-6.733717683" " 1.961201804" " 0.000000e+00" " 0.000000e+00" KLHL41 "KLHL41" "-6.656135657" " 3.030242156" " 0.000000e+00" " 0.000000e+00"
> rownames(rt)=rt[,1]
> exp=rt[,2:ncol(rt)]
> dimnames=list(rownames(exp),colnames(exp))
> data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
> diffExpLevel=avereps(data)
> hmExp=log10(diffExpLevel+0.00001)
> hmMat=as.matrix(hmExp)
> pdf(file="heatmap.pdf",height=7,width=7)
> par(oma=c(3,3,3,5))
> heatmap.2(hmMat,col='greenred',trace="none",cexCol=1)
> dev.off()
> rt=read.table("alldiff.txt",sep="\t",header=T,check.names=F)
> rt=as.matrix(rt)
> rownames(rt)=rt[,1]
> exp=rt[,2:ncol(rt)]
> dimnames=list(rownames(exp),colnames(exp))
> data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
> allDiff=avereps(data)
> allDiff=as.data.frame(allDiff)
> pdf(file="vol.pdf")
#不太好看。
> xMax=max(-log10(allDiff$adj.P.Val+0.00001))
> yMax=max(abs(allDiff$logFC))
> plot(-log10(allDiff$adj.P.Val), allDiff$logFC, xlab="-log10(adj.P.Val)",ylab="logFC",main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
> diffSub=subset(allDiff, allDiff$adj.P.Val<0.05 & abs(allDiff$logFC)>1)
> points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.4)
> abline(h=0,lty=2,lwd=3)
> dev.off()
> data = read.table("alldiff.txt", header=TRUE)
> head(data,10)
id logFC logCPM PValue adj.P.Val 1 CKM -8.576280 5.4592297 0 0 2 ACTA1 -7.228533 6.5359265 0 0 3 MYLPF -7.209576 2.6186547 0 0 4 PYGM -7.202706 4.0748302 0 0 5 SLN -6.733718 1.9612018 0 0 6 KLHL41 -6.656136 3.0302422 0 0 7 MYOT -6.612275 0.8828837 0 0 8 TNNC2 -6.600796 3.1647841 0 0 9 ACTN3 -6.583526 0.9872158 0 0 10 NEB -6.500059 5.1964626 0 0
> logFC_cutoff <- with(data,mean(abs(logFC)) + 2*sd(abs(logFC)))
> logFC_cutoff
> data$sig = as.factor(ifelse(data$adj.P.Val < 0.05 & abs(data$logFC) > logFC_cutoff,ifelse(data$logFC > logFC_cutoff ,'UP','DOWN'),'NOT'))
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),'\nThe number of up gene is ',nrow(data[data$sig =='UP',]),'\nThe number of down gene is ',nrow(data[data$sig =='DOWN',]))
> g = ggplot(data=data, aes(x=logFC, y=-log10(adj.P.Val), color=sig)) +
geom_point(alpha=0.4, size=2.5) +
theme_bw(base_size=15)+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))+
xlab("logFC") + ylab("-Lgp") +
ggtitle( this_tile ) +
theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red'))
g
> g2<- g+ geom_hline(yintercept=-log10(0.05),colour="black", linetype="dashed") +
geom_vline(xintercept=c(-1,1),colour="black", linetype="dashed")
> g2
交流学习
