在多组学联合分析中，需要用得到九象限图来对两个组学获得得基因结果进行可视化，例如下面这样得，因此这两天主要是对这个内容进行整理。
代码解析
代码语言：javascript复制#===========================================================
# 1. 读入两个差异分析表（分别是转录组和翻译组）并合并；
##==========================================================
RNA =read.csv("rna.csv")
Ribo=read.csv("lribo.csv")
dim(RNA)
dim(Ribo)
##合并两个表格。选择表头为"id"的列进行合并；
combine= merge(RNA,Ribo,
by.x="id",
by.y="id",
suffixes = c(".RNA",".Ribo") ,
all.x=F,
all.y=F)
dim(combine)
#保存合并后的数据；
write.table(combine,"merge.txt",sep="t",row.names=FALSE)
#====================
# 2. 数据准备
#====================
dat = read.table("merge.txt",header=T,sep="t")
# 查看表头;
colnames(dat)
#取RNA的差异倍数，将作为X轴;
#取Ribo的差异倍数，将作为Y轴;
x=dat$RNA_log2FC
y=dat$H3K27me3_log2FC
#计算两个组学差异倍数的相关性，并取2位小数;
cor = round(cor(x,y),2)
#准备作为图形的标题;
main = paste("correlation=",cor,sep="")
main
#数据筛选：
#RNAdown:转录组显著下调的行;
RNAdown<-intersect(which(dat$RNA_log2FC <= -2),
which(dat$RNA_Pvalue <= 0.05))
#RNAup：转录组显著上调的行;
RNAup<-intersect(which(dat$RNA_log2FC >= 2),
which(dat$RNA_Pvalue <=0.05))
# RNA显著差异的行；
RNA_diff = intersect(which(abs(dat$RNA_log2FC) >= 2),
which(dat$RNA_Pvalue <=0.05))
# RNA不显著差异的行；
RNA_nodiff = union(which(abs(dat$RNA_log2FC) < 2),
which(dat$RNA_Pvalue > 0.05))
#Ribo显著下调的行 ；
Ribdown<-intersect(which(dat$Ribo_log2FC <= -2),
which(dat$RNA_Pvalue <= 0.05))
#Ribo显著上调的行；
Ribup<-intersect(which(dat$H3K27me3_log2FC >= 2),
which(dat$H3K27me3_Pvalue <=0.05))
# Ribo显著差异的行；
Ribo_diff = intersect(which(abs(dat$H3K27me3_log2FC) >= 2),
which(dat$H3K27me3_Pvalue <=0.05))
# Ribo不显著差异的行；
Ribo_nodiff = union(which(abs(dat$H3K27me3_log2FC) < 2),
which(dat$H3K27me3_Pvalue > 0.05))
# 两个组学都下调；
samedown<-intersect(RNAdown,Ribdown)
# 两个组学都上调；
sameup<-intersect(RNAup,Ribup)
# RNA下调，Ribo上调；
diff1<-intersect(RNAdown,Ribup)
# RNA上调，Ribo下调；
diff2<-intersect(RNAup,Ribdown)
#差异方向相同的基因数量；
homo=length(union(sameup,samedown))
#差异方向相反的基因数量；
opp=length(union(diff1,diff2))
# 仅仅RNA有差异的基因数量；
RNAonly=length(intersect(Ribo_nodiff,RNA_diff))
# 仅仅Ribo有差异的基因数量；
Riboonly=length(intersect(Ribo_diff,RNA_nodiff))
#生成图例所用标签；
RNA_change=paste("Transcription(",RNAonly,")",sep="")
Ribo_change=paste("H3K27me3(",Riboonly,")",sep="")
homo_change=paste("Homodirection(",homo,")",sep="")
opp_change=paste("Opposite(",opp,")",sep="")
#给图例准备颜色
id<-c(RNA_change,Ribo_change,homo_change,opp_change)
co<-c("red","blue","brown","yellow")
#=====================
# 3. 图形绘制
#=====================
# 绘制所有点散点图，注意可以修改X轴和y轴范围；
# 可利用type = "n",记得改坐标轴范围大小
plot(x,y,pch=16,cex = 1,
xlim=c(-9,6),ylim=c(-8,8),
cex.lab=1.2,main=main ,col="grey80",
xlab="gene expression",
ylab="ribo expression")
#利用次级函数绘制不同象限的图，并加入图例；
points(x[RNA_diff],y[RNA_diff],pch=16,cex=0.3,col="red")
points(x[Ribo_diff],y[Ribo_diff],pch=16,cex=0.3,col="blue")
points(x[c(sameup,samedown)],y[c(sameup,samedown)],pch=16,cex=0.3,col="#brown")
points(x[c(diff1,diff2)],y[c(diff1,diff2)],pch=16,cex=0.3,col="yellow")
legend("topright",legend=id,col=co,pch=16,bty = "n",cex = 0.8,text.font = 2)
#dev.off()总结
代码里面很多是参数是可以换的，其实还是数据得质量怎么样，如果数据质量好一些，图片结果应该还是很漂亮得，上面得图得代码还是基本得，需要结合研究得内容，来进行具体得调整，最近明白一个道理，生信得每一步都是走细节，以前只是觉得代码能通就可以了，现在觉得哪里都是问题，都要一点一点得理清楚，所以希望我们不要做调包得那个，要自己去理解其中得含义。
