生信分析之R语言绘图整理

#####################################
##              第二节            ##
#####################################
data() #查看数据库
class(rivers) #"numeric"，
plot(rivers) #散点图
hist(rivers) #直方图
class(state.abb) #"character"字符串
pie(table(state.abb))
pie(table(state.region))
yanse <- c("red","green","blue","yellow")

pie(table(state.region),col = yanse) #颜色
labs <- c("NC","S","N","W")
pie(table(state.region),labels = labs) #标签
class(state.division) #"factor" 因子-分组
plot(mtcars$cyl,mtcars$disp) #"numeric" #散点图
plot(as.factor(mtcars$disp),mtcars$cyl) #箱线图
class(state.x77) #"matrix" "array"
heatmap(state.x77) #矩阵-》热图
heatmap(as.data.frame(state.x77)) #报错，'x' must be a numeric matrix
fit <- lm(weight~height,data=women)
class(fit) #"lm"
plot(fit)




#####################################
##              第三节             ##
#####################################
example("barplot") #案例
example("heatmap") #案例
library(ggplot2)
example("qplot")
help("barplot")
??barplot


#####################################
##          条形图 barplot         ##
#####################################
#输入的数据是向量
getwd()
setwd("/Users/marscolono/Desktop/Macmini 工作数据包/R代码整理/R语言绘图配套资料")
m <- read.csv("Rdata/homo_length.csv",header = T)
m
class(m) #"data.frame"
head(m)
x <- m[1:24,]
x
class(x$length) #"integer"
barplot(height = x$length) #height
barplot(height = x$length,names.arg = x$chr,las=3) #设置横坐标,las=3设置方向
yanse <- sample(colours(),24,replace = F) #随机抽取24个颜色
barplot(height = x$length,names.arg = x$chr,col = yanse) #设置颜色
barplot(height = x$length,names.arg = x$chr,col = rainbow(7)) #彩虹颜色
barplot(height = x$length,names.arg = x$chr,col = rainbow(7),horiz =T) # horiz 设置方向
barplot(height = x$length,names.arg = x$chr,col = rainbow(7),border = F) #是否加边框
barplot(height = x$length,names.arg = x$chr,col = rainbow(7),border = F,
        main = "Human chromosome length distribution",xlab = "Chromosome Name",ylab = "Chromosome Length")
#main：标题，xlab，x轴标签，ylab ，y轴标签



#####################################
##             分组条形图         ##
#####################################
#输入的数据是矩阵
x <- read.csv("Rdata/sv_distrubution.csv",header = T,row.names = 1) #row.names=1,第一列作为行名
x #"data.frame"
barplot(x)
#Error in barplot.default(x) : 'height' must be a vector or a matrix
barplot(as.matrix(x)) #必须矩阵 #堆砌条形图
barplot(t(as.matrix(x))) #t来转换，将矩阵转置
barplot(t(as.matrix(x)),col = rainbow(4)) #加颜色，不同的突变显示不同的颜色
barplot(t(as.matrix(x)),col = rainbow(4),beside = T) #beside分组条形图
barplot(t(as.matrix(x)),col = rainbow(4),legend.text = colnames(x)) #添加图例 legend.text = colnames(x) #列名
barplot(t(as.matrix(x)),col = rainbow(4),legend.text = colnames(x),ylim = c(0,800)) #设置y坐标范围，记得加c
barplot(t(as.matrix(x)),col = rainbow(4),legend.text = colnames(x),ylim = c(0,800),
        main = "SV Distribution",xlab="Chromosome Number",ylab="SV Numbers")


#####################################
##          直方图            ##
#####################################
#基因组注释文件
x <- read.table("Rdata/H37Rv.gff",sep = "\t",header = F,skip = 7,quote = "") 
#sep = "\t"，#分隔
#skip = 7，#跳过前面7行
#quote = ""，表示里面的字符串不用引号括起来

x <- x[x$V3=="gene",]

x <- abs(x$V5-x$V4+1) #第5列减去第四列
length(x)
range(x) #看一下最小值，最大值
hist(x) #直方图
hist(x,breaks = 80) #breaks，设置条目的数目
hist(x,breaks = c(0,500,1000,1500,2000,2500,15000))
hist(x,breaks = 80,freq = F) #freq，根据概率密度画图
hist(x,breaks = 80,density = T) #density表示是否绘制斜线
hist(rivers,density = T,breaks = 10)
?hist

#添加正太分布的曲线
#正太函数的值保存为一个变量
h=hist(x,nclass=80,col="pink",xlab="Gene Length (bp)",main="Histogram of Gene Length");
h
rug(x) #绘制分布的密度
xfit<-seq(min(x),max(x),length=100) #将基因长度平均分割成100份
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x)) #生成基因长度的正态分布
yfit <- yfit*diff(h$mids[1:2])*length(x) #计算正态分布的值
lines(xfit, yfit, col="blue", lwd=2) #绘制

#####################################
##              散点图             ##
#####################################
m <- read.table("Rdata/prok_representative.txt",sep="\t")
genome_size <- m[,2]
gene_size <- m[,4]
plot(genome_size,gene_size)
plot(genome_size,gene_size,pch=16,cex=0.5,xlab="Genome Size",ylab="Genes") #pch设置形状
#添加回归曲线
fit <- lm(gene_size ~ genome_size) #回归曲线
abline( fit,col="blue",lwd=1.8 )
summary(fit)
text(12,10000,labels = 'y=843.691x+286.598\nR2=0.97') #人为添加公式

rr <- round( summary(fit)$adj.r.squared,2) #r方
intercept <- round( summary(fit)$coefficients[1],2) #截距值
slope <- round( summary(fit)$coefficients[2],2) #系数项
eq <- bquote( atop( "y = " * .(slope) * " x + " * .(intercept), R^2 == .(rr) ) ) #组合
text(12,6e3,eq)



#####################################
##            饼图               ##
#####################################
x <- read.csv("Rdata/homo_length.csv",header = T)
#比较染色体大小
x <- x[1:24,]
barplot(height = x$length,names.arg = x$chr)
pie(x$length/sum(x$length))

m <- read.table("Rdata/Species.txt")
m  #"data.frame"
x <- m[,3]
pie(x)
pie(x,col=rainbow(length(x)))
lbls <- paste(m[,1],m[,2],"\n",m[,3],"%" )
pie(x,col=rainbow(length(x)),labels = lbls)
pie(x,col=rainbow(length(x)),labels = lbls)
pie(x,col=rainbow(length(x)),labels = lbls,radius = 1)
pie(x,col=rainbow(length(x)),labels = lbls,radius = 1,cex=0.8)

#3D饼图
install.packages("plotrix")
library(plotrix)
pie3D(x,col=rainbow(length(x)),labels = lbls)
pie3D(x,col=rainbow(length(x)),labels = lbls,cex=0.8)
pieplot <- pie3D(x,col=rainbow(length(x)),radius = 1,explode = 0.1) # explode 分开
pie3D.labels(pieplot,labels = lbls,labelcex = 0.8,height = 0.1,labelrad = 1.75)

#扇形图
fan.plot(x,col=rainbow(length(x)),labels = lbls,cex=0.8,radius = 1)



#####################################
##        箱线图            ##
#####################################
boxplot(mpg ~cyl,data = mtcars) #数值向量～分组的因子
#RNAseq表达数据，每个值表示基因的表达量
m <- read.csv("Rdata/gene_expression.csv",header = T,row.names = 1)
head(m)
boxplot(m)
boxplot(m,outline=F) #把离群的点去掉
install.packages("reshape2")
library(reshape2)
x <- melt(m) #修改
head(x)
boxplot(value ~ variable,data = x)
boxplot(value ~ variable,data = x,outline=F)

boxplot(len ~ dose, data = ToothGrowth)
boxplot(len ~ dose:supp, data = ToothGrowth) #分组
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"))#颜色
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"),notch=T) #notch添加缺口
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"),width=c(1,0.5,1,0.5,1,0.5))#width设置每个盒子的款度
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"),varwidth=T) #宽度可变
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"),boxwex=0.5) #缩放
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"),staplewex=0.5) #设置盒子上下盖子线的宽度
boxplot(len ~ dose:supp, data = ToothGrowth,col = c("orange", "red"),sep = ":",lex.order = T) #sep = ":"设置连接，lex.order排列顺序


###########################
##        par颜色       ##
###########################
#par参数集合
opar <- par(no.readonly = T)
?par
par() #集合#列表
x <- par()
x$pch
par(pch=16)
opar <- par(no.readonly = T) #保存一份
opar
par(opar) #恢复

plot(women$height,women$weight)
plot(women$height,women$weight,type = "b",pch=16,col="red",lty=2,lwd=2,
     main = "Main Title",sub = "Sub Title",xlab = "Heigth",
     ylab = "Weight",xlim = c(50,100),ylim = c(100,200),cex=1.5,las=3,
     adj=0.3,bty="c",fg="blue",tck=-0.03,col.main="green",cex.lab=2)



#####################################
##            颜色             ##
#####################################
colors()
sample(colours(),24)
rgb(0, 1, 0) #将数值转换成RBG
hsv(0,1,0) #色调/饱和度/值
hcl(0,1,0) #色调/   /亮度

rainbow(7)
pie(1:7,col=rainbow(7))
heat.colors(7) # 暖色调
pie(1:7,col=heat.colors(7))
terrain.colors(7) # 
pie(1:7,col=terrain.colors(7))
topo.colors(7)
pie(1:7,col=topo.colors(7))
cm.colors(7)
pie(1:7,col=cm.colors(7))
gray.colors(7)
pie(1:7,col=gray.colors(7))

install.packages("RColorBrewer")
library(RColorBrewer)
display.brewer.all() #所有颜色列表
display.brewer.pal(name="Set1",n=9) #Set1的调色板
brewer.pal.info
brewer.pal(name = "Set1",n=7)

#########################
##       热图          ##
#########################
#1 heatmap()
#2 gplots ?? heatmap.2()
#3 pheatmap?? pheatmap
#4 ComplexHeatmap??

#输入数据必须为数值矩阵
example("heatmap") #内置heatmap
m <- read.csv("Rdata/heatmap.csv",header = T, row.names = 1) #名字都列好
class(m) #"data.frame"
x <- as.matrix(m) #转换格式
heatmap(x)
heatmap(t(x)) #互换行和列
heatmap(x,col=c("green","red")) #修改颜色
yanse <- colorRampPalette(c("red","black","green"))(nrow(x)) #colorRampPalette随机生成渐变色
heatmap(x,col=yanse)
heatmap(x,col=yanse,RowSideColors = yanse)  #设置颜色条RowSideColors
heatmap(x,col=yanse,ColSideColors = colorRampPalette(c("red","black","green"))(ncol(x))) #设置颜色条ColSideColors
heatmap(x,col=yanse,Rowv = NA) #Rowv=NA 按行聚类取消
heatmap(x,col=yanse,Rowv = NA,Colv = NA)
heatmap(state.x77)
heatmap(state.x77,scale = "col") #按列标准化

#####################################
##        热图2           ##
#####################################
#install.packages("gplots")
library(gplots)
heatmap.2(x)
#example(heatmap.2) #绘图风格，多了图例
heatmap.2(x)
heatmap.2(x,key = F) #是否显示图例key
heatmap.2(x,symkey = F) #颜色是否对chen
heatmap.2(x,symkey = T,density.info = "none")
heatmap.2(x,symkey = T,trace = "none") #是否添加追踪的线
heatmap.2(x,symkey = T,tracecol = "black") #设置追踪线的颜色
heatmap.2(x,dendrogram = "none")
heatmap.2(x,dendrogram = "row")
heatmap.2(x,dendrogram = "col")

#install.packages("pheatmap")
library(pheatmap)
example("pheatmap")
pheatmap(x)
annotation_col <- data.frame(CellType=factor(rep(c("N1","T1","N2","T2"),each=5)))
rownames(annotation_col) <- colnames(x)
pheatmap(x,annotation_col = annotation_col)
pheatmap(x,annotation_col = annotation_col,display_numbers = T) #在热土上显示数字
pheatmap(x,annotation_col = annotation_col,display_numbers = T,number_format = "%.2f") #设置数字格式
pheatmap(x,annotation_col = annotation_col,display_numbers = T,number_format = "%.1f",number_color = "black") #

#####################################
##        韦恩图           ##
#####################################
#install.packages("VennDiagram")
setwd("/Users/marscolono/Desktop/Macmini 工作数据包/R代码整理/R语言绘图配套资料")
listA <- read.csv("Rdata/genes_list_A.txt",header=FALSE)
A <- listA$V1
listB <- read.csv("Rdata/genes_list_B.txt",header=FALSE)
B <- listB$V1
listC <- read.csv("Rdata/genes_list_C.txt",header=FALSE)
C <- listC$V1
listD <- read.csv("Rdata/genes_list_D.txt",header=FALSE)
D <- listD$V1
listE <- read.csv("Rdata/genes_list_E.txt",header=FALSE)
E <- listE$V1
length(A);length(B);length(C);length(D);length(E)
library(VennDiagram)
venn.diagram(list(C = C, D = D),fill = c("yellow","cyan"), cex = 1.5,filename = "venn2.png")
venn.diagram(list(A = A, C = C, D = D), fill = c("yellow","red","cyan"), cex = 1.5,filename="venn3.png")
venn.diagram(list(A = A, B = B, C = C, D = D), fill = c("yellow","red","cyan","forestgreen"), cex = 1.5,filename="venn4.png")
venn.diagram(list(A = A, B = B, C = C, D = D , E = E ), fill = c("yellow","red","cyan","forestgreen","lightblue"), cex = 1.5,filename="venn5.png")



#############################
##        曼哈顿图         ##
#############################
#全基因组关联分析 / 大数据可视化
setwd("/Users/marscolono/Desktop/Macmini 工作数据包/R代码整理/R语言绘图配套资料")
#install.packages("qqman")
library(qqman)
library(RColorBrewer)
str(gwasResults) #自带模拟数据
head(gwasResults)
manhattan(gwasResults)
manhattan(gwasResults, main = "Manhattan Plot", ylim = c(0, 6), cex = 0.6,cex.axis = 0.9, col = c("blue4", "orange3"), suggestiveline =
            F, genomewideline = F,chrlabs = c(1:20, "P", "Q"))

unique(gwasResults$CHR)
number <- length(unique(gwasResults$CHR))
yanse <- brewer.pal(n = 4,name = "Set1")
manhattan(gwasResults, main = "Manhattan Plot", ylim = c(0, 6), cex = 0.6,cex.axis = 0.9, col = yanse, suggestiveline =
            F, genomewideline = F,chrlabs = c(1:20, "P", "Q"))


manhattan(subset(gwasResults,CHR==3)) #单独看一下想要的染色体
#高亮显示部分SNP
snpsOfInterest
manhattan(gwasResults, highlight = snpsOfInterest)#高亮

#注释SNP结果
manhattan(gwasResults, annotatePval = 0.001) #注释
manhattan(gwasResults, annotatePval = 0.001, annotateTop = FALSE) #top最显著的点
#???????ݿ??Բ鿴manhattan??qqman?İ????ĵ?
help("manhattan")
vignette("qqman")

##############################
##        火山图            ##
##############################
#基因表达分析中
m <- read.csv("Rdata/Deseq2.csv",header = T,row.names = 1)
head(m)
m <- na.omit(m) #去除包含缺失值的行
plot(m$log2FoldChange,m$padj)
plot(m$log2FoldChange,-1*log10(m$padj))
plot(m$log2FoldChange,-1*log10(m$padj),xlim = c(-10,10),ylim=c(0,100))

m <- transform(m,padj=-1*log10(m$padj)) #对数据框的列进行修改

#将数据分成三组
down <- m[m$log2FoldChange<=-1,] 
up <- m[m$log2FoldChange>=1,]
no <- m[m$log2FoldChange>-1 & m$log2FoldChange <1,] 
#分别画图
plot(no$log2FoldChange,no$padj,xlim = c(-10,10),ylim=c(0,100),col="blue",pch=16,cex=0.8,main = "Gene Expression",xlab = "log2FoldChange",ylab="-log10(Qvalue)")
points(up$log2FoldChange,up$padj,col="red",pch=16,cex=0.8) #其余的点加上去
points(down$log2FoldChange,down$padj,col="green",pch=16,cex=0.8) #其余的点加上去


#############################
##        GCdepth图        ##
#############################
#评估拼接的基因组效果，散点图+直方图
opar <- par(no.readonly = T)
par(mfrow=c(2,3)) #布局，2行3列
plot(pressure)
m <- read.table("Rdata/GC-depth.txt");
head(m)

#布局layout
nf <- layout(matrix(c(0,2,0,0,1,3),2,3,byrow=T),c(0.5,3,1),c(1,3,0.5),TRUE);
par(mar=c(5,5,0.5,0.5))
layout.show(3) #分割好了
par("mar")
par(mar=c(5,5,0.5,0.5)) #反复调整之后得到的，mar设置
x <- m[,1]
y <- m[,2]
#散点图
plot(x,y,xlab='GC Content(%)',ylab='Depth',pch=16,col="#FF000077",xlim=c(0,100),ylim=c(0,max(y)))
#hist(x,breaks = 100) #直方图

xhist <- hist(x,breaks=100,plot=FALSE) # 保存到一个变量里
yhist <- hist(y,breaks=floor(max(y)-0),plot=FALSE) # 保存到一个变量里

par(mar=c(0,5,1,1))#重新生成一个布局的内容
barplot(xhist$counts,space=0,xlim=c(0,100) )#最上侧的条形图

par(mar=c(5,0,1,1));#最右侧的条形图
barplot(yhist$counts,space=0,horiz=TRUE,ylim=c(0,max(y) ) );


#####################################
##        COG功能注释图         ##
#####################################
#蛋白质
m <- read.table("Rdata/cog.class.annot.txt",head=T,sep="\t");
head(m)
layout(matrix(c(1,2),nr=1),widths=c(20,13));#进行布局
layout.show(2)
par( mar=c(3,4,4,1)+0.1 );#调整一下边界值

#创建一个class
class <- c("J","A","K","L","B","D","Y","V","T","M","N","Z","W","U","O","C","G","E","F","H","I","P","Q","R","S");
t <- factor( as.character(m$Code),levels=class )
m <- m[order(t),] #对数据进行排序

#对第三列进行绘图
barplot(m$Gene.Number,space=F,col=rainbow(25),ylab="Number of genes",names.arg = m$Code)

#将分组信息添加上：
l <- c(0,5,15,23,25);
id<- c("INFORMATION STORAGE\nAND PROCESSING","CELLULAR PROCESSES\nAND SIGNALING","METABOLISM","POORLY\nCHARACTERIZED")
abline( v = l[c(-1,-5)]) #绘制斜线

#将分类的信息添加到图中
for( i in 2:length(l) ){
  text( (l[i-1]+l[i])/2,max(m[,3])*1.1,id[i-1],cex=0.8,xpd=T)
}

par(mar=c(2,0,2,1)+0.1 ); #设置边界
#绘制空的图，设定xy轴的范围
plot(0,0,type="n",xlim=c(0,1),ylim=c(0,26),bty="n",axes=F,xlab="",ylab="") #axes=F不显示坐标轴

for( i in 1:length(class) ){
  text(0,26-i+0.5,paste(m$Code[i],m$Functional.Categories[i]),pos=4,cex=1,pty=T)
}
title(main = "COG function classification"); #添加标题




#####################################
##        GO功能注释条目图          ##
#####################################
#富集分析

library(ggplot2)

go <- read.csv("Rdata/go.csv",header = T)
head(go)
go_sort <- go[order(go$Ontology,-go$Percentage),] #对数据框进行排序
m <- go_sort[go_sort$Ontology=="Molecular function",][1:10,] #前10个term
c <- go_sort[go_sort$Ontology=="Cellular component",][1:10,]
b <- go_sort[go_sort$Ontology=="Biological process",][1:10,]
slimgo <- rbind(b,c,m) #合并

#将Tremת转换为因子
slimgo$Term=factor(slimgo$Term,levels=slimgo$Term) #转化为因子

colnames(slimgo) #输出列名
pp <- ggplot(data = slimgo, mapping = aes(x=Term,y=Percentage,fill=Ontology)) #aes映射
pp #输出映射的内容
pp+geom_bar(stat="identity") #分组式的条目图
pp+geom_bar(stat="identity")+coord_flip() #coord_flip()翻转坐标轴
pp+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(slimgo$Term))) #从大到小排序 scale_x_discrete
pp+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(slimgo$Term)))+guides(fill=FALSE) #去掉图例
pp+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(slimgo$Term)))+guides(fill=FALSE)+theme_bw() #切换成黑白背景主题


#####################################
##        KEGG功能注释图          ##
#####################################
#得到基因差异的信号通路
library(ggplot2)

pathway <-  read.csv("Rdata/kegg.csv",header=T)
head(pathway)
colnames(pathway)#设置列名

pp <-  ggplot(data=pathway,mapping = aes(x=richFactor,y=Pathway))
pp #完成映射
pp + geom_point() #绘制点图
pp + geom_point(aes(size=AvsB)) #设置点的大小
pp + geom_point(aes(size=AvsB,color=Qvalue)) #设置颜色
pp + geom_point(aes(size=AvsB,color=Qvalue)) + scale_colour_gradient(low="green",high="red")#修改默认颜色
#修改标题
pp + geom_point(aes(size=AvsB,color=Qvalue)) + scale_colour_gradient(low="green",high="red")+labs(title="Top20 of pathway enrichment",x="Rich factor",y="Pathway Name",color="-log10(Qvalue)",size="Gene Numbers")
pp + geom_point(aes(size=AvsB,color=Qvalue)) + scale_colour_gradient(low="green",high="red")+labs(title="Top20 of pathway enrichment",x="Rich factor",y="Pathway Name",color="-log10(Qvalue)",size="Gene Numbers")+theme_bw()