R语言：热图，火山图

> 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

交流学习