这个文章跟着之前的文章完整转录组RNAseq分析流程（tophat2+cufflink+cuffdiff）用了之前分析的 cuffdiff 的结果,参考视频 https://www.bilibili.com/video/BV1gW411Y7Qf
文中用到的hic的数据由于是别人的东西，就不方便放出来，可以看一下孟叔视频，加群后可以在群文件下载
###############################################
#rnaseq 可视化 孟b站 R绘图
#############################################
#content
#1. volcano plot
#2. heatmap plot
#3. cytogram plot
#4. plot with layout
#############################################
rm(list = ls())
#############################################
#part one
#1.input table cuffdiff
#############################################
#volcano plot
#############################################
cuffdiff_result = read.table(file = "../Desktop/test_data/rnaseq_test_date/diff_out1/gene_exp.diff",header = T)
ctrl_fpkm = cuffdiff_result$value_1
treat_fpkm = cuffdiff_result$value_2
log2_foldchange = log2(treat_fpkm / ctrl_fpkm)
#此时会出现NaN（无穷）,除数为0 -inf(负无穷)，被除数为0，log2无法计算
log2_foldchange[ctrl_fpkm==0]=0 #考虑ctrl组有等于0的情况剔除，将其赋值为0
log2_foldchange[treat_fpkm==0]=0 #同理
log10_p_value = log10(cuffdiff_result$p_value)*-1
######################################## 1st edition
plot(x=log2_foldchange,y=log10_p_value)
# 1.p_value 等于0的值一般不画
# 2.显著和不显著的颜色要标出
plot(x=log2_foldchange,y=log10_p_value,xlim = c(-4,4),ylim = c(0.001,4))
#限制坐标轴，x轴-4到4 ，y轴 0.001到4
######################################## 2nd edition
log10_p_value.filter = log10_p_value[log10_p_value >= 0.001]
log2_foldchange.filter = log2_foldchange[log10_p_value >= 0.001]
plot(x=log2_foldchange.filter,y=log10_p_value.filter,xlim = c(-4,4),ylim = c(0,4))
######################################## 3rd edition
#上色,实心
plot(x=log2_foldchange.filter,y=log10_p_value.filter,
xlim = c(-4,4),ylim = c(0,4),
col = rgb(0,0,1,0.1),pch=16)
####################################### 4th edition
#将显著性的点标出来,需要写一个判断
length(log2_foldchange.filter)
col_vector = rep(rgb(0,0,1,0.1),length(log2_foldchange.filter))
#选出大于等于-1*log10（0.001）的点将其标为红色
col_vector[log10_p_value.filter >= -1*log10(0.001)] = rgb(1,0,0)
plot(x=log2_foldchange.filter,y=log10_p_value.filter,
xlim = c(-4,4),ylim = c(0,4),
col = col_vector,pch=16)
######################################## 5th edition
#真实情况下考虑的情况很多
#1. p-value <=0.05
#2. ctrl,treat FPKM >0
#3. foldchange > 2 | <0.5
select_sign_vector = (cuffdiff_result$value_1 > 0) & (cuffdiff_result$value_2 > 0) & (abs(log2_foldchange) >= 1) & (cuffdiff_result$value_1 >=1 | cuffdiff_result$value_2 >=1) & (cuffdiff_result$p_value <=0.05)
log10_p_value.filter = log10_p_value[log10_p_value >= 0.001]
log2_foldchange.filter = log2_foldchange[log10_p_value >= 0.001]
select_sign_vector.filter = select_sign_vector[log10_p_value >= 0.001]
col_vector = rep(rgb(0,0,1,0.1),length(log2_foldchange.filter))
col_vector[select_sign_vector.filter] = rgb(1,0,0)
plot(x=log2_foldchange.filter,y=log10_p_value.filter,
xlim = c(-4,4),ylim = c(0,4),
col = col_vector,pch=16)
#加辅助线
abline(h = -1*log10(0.05),lwd =2 , lty = 3 , col="#4C5B61")
####################################################################
# part two
# heatmap
#######################################################################
#step 1 画方框
plot(x = c(1:10),y=c(1:10),type="n")
# plot (0,0) (1,0) (1,1) (0,1)
rect(xleft = 0 , ybottom = 0,xright = 1 , ytop = 1,col = rgb(1,0,0))
rect(xleft = 5 , ybottom = 5,xright = 6 , ytop = 6,col = rgb(0,1,0))
# input matrix
# xleft, ybottom,xright,ytop
# color
# output image
input_matrix = matrix(c(1:36),6,6)
#set image size
x_size = dim(input_matrix)[1]
y_size = dim(input_matrix)[2]
#如果，我们从0，0 向6，6方向找规律
#先画，第一行，再画第二行
# 1st row
# xleft 0,1,2,3,4,5
# ybottom 0,0,0,0,0,0
# xright 1,2,3,4,5,6
# ytop 1,1,1,1,1,1
#2nd row
# xleft 0,1,2,3,4,5
# ybottom 1,1,1,1,1,1
# xright 1,2,3,4,5,6
# ytop 2,2,2,2,2,2
#其他同理
my_xleft = rep(c(0:(x_size-1)),y_size)
my_xright = my_xleft + 1
my_ybottom = rep(c(0:(y_size - 1)),each = y_size)
my_ytop = my_ybottom + 1
plot(x=c(0:x_size),y=c(0:y_size),type="n")
rect(xleft = my_xleft ,ybottom = my_ybottom ,xright = my_xright,ytop = my_ytop)
# 下一步，确定颜色
# 先画最简单的，最小值是白色，最大值是红色，中间线性变换
mat.max = max(input_matrix)
input_matrix.rate = input_matrix / mat.max
col.mat = rgb(1,0,0,input_matrix.rate)
plot(x=c(0:x_size),y=c(0:y_size),type="n")
rect(xleft = my_xleft ,ybottom = my_ybottom ,xright = my_xright,ytop = my_ytop,col = col.mat)
######此时画的热图与表格对应不上，应该将表格转置
##去掉方格线
plot(x=c(0:x_size),y=c(0:y_size),type="n",xlab = "",ylab = "",axes = F,frame.plot = F)
rect(xleft = my_xleft ,ybottom = my_ybottom ,xright = my_xright,ytop = my_ytop,col = col.mat,border = F)
# 真实的Hi—c数据
hic_mat.raw = read.table(file="../Desktop/test_data/R基础绘图（live5）/20171203-Live-R_partII/data_file/chr_16_100000_MAPQ20.txt",sep ="," ,header=F)
input_matrix = as.matrix(hic_mat.raw)
x_size = dim(input_matrix)[1]
y_size = dim(input_matrix)[2]
my_xleft = rep(c(0:(x_size-1)),each = x_size)
my_xright = my_xleft + 1
my_ybottom = rep(c((y_size - 1):0),y_size)
my_ytop = my_ybottom + 1
mat.max = quantile(input_matrix,prob=0.95)
input_matrix.rate = input_matrix / mat.max
input_matrix.rate[input_matrix.rate>1]=1
col.mat = rgb(1,0,0,as.vector(as.matrix(input_matrix.rate)))
plot(x=c(0:x_size),y=c(0:y_size),type="n",xlab = "",ylab = "",axes = F,frame.plot = F)
rect(xleft = my_xleft ,ybottom = my_ybottom ,xright = my_xright,ytop = my_ytop,col = col.mat,border = F)
# 为了快速画图，我们尝试直接把图片生成
png(file="../Desktop/test_hic.png",width = 2000,height = 2000)
plot(x=c(0:x_size),y=c(0:y_size),type="n",xlab = "",ylab = "",axes = F,frame.plot = F)
rect(xleft = my_xleft ,ybottom = my_ybottom ,xright = my_xright,ytop = my_ytop,col = col.mat,border = F)
#记得关掉绘图函数
dev.off()
网友评论