美文网首页生信学习
R入门与基础绘图系统

R入门与基础绘图系统

作者: 黄晶_id | 来源:发表于2019-04-27 18:19 被阅读128次

这是我听B站鲮鱼不会飞视频(R入门与基础绘图系统)里的笔记哦~
加个目录吧,因为实在太长了,怕找不到重点:
一、前面是R的基础介绍,特别基础的那种
二、用基础绘图函数做了火山图,有每一步代码的解释
三、用基础函数做了热图——Hi-C可视化,有每一步代码的解释

问:为什么有那么多现成的包可以做热图和火山图,我们还是要用基础函数自己写?
答:就当在练习和熟悉这些基础作图函数吧~其实当遇到特别奇葩的审稿人让你微调图的时候,基础函数特别有优势。

做生物信息学最需要接触3类包

1.CRN 中的R包——CRN是R的官方网站

安装命令:install.packages("package name")

2.bioconductor上的R包——Bioconductor是一个专门做生信R包的平台,里面发布了各种生信分析的R包

安装命令:biocLite("package name")

特别注意:Bioconductor上的course,特别有用。就是Bioconductor上有自己的课程,这些课程不但给代码而且给PPT有的时候还会录视频。这里面单细胞测序的、甲基化数据分析的、CNV、SNP、RNA-Seq分析的....全都有,非常全,它会告诉你用什么包画什么图怎么解决问题,代码非常详细,最最最最重要的是这里面所有的course完全免费。

3.大神写的包在GitHub上
根据自己的情况选择性的使用


今天我们学习R里基础绘图系统

问:有了高级绘图系统为什么还是要学会基础绘图系统?
答:高级绘图系统更多的是快速探索数据的性质,这是他的优势;对于一些精修要画可发表的图,对格式的要求非常多,这时使用基础绘图包非常方便。

R的绘图系统

  • R的基础绘图系统:
  1. 基础绘图系统中的低级绘图函数:创建画布,点,线,多边形
  2. 基础绘图系统中的高级绘图函数:plot(), boxplot(), hist(), density() 等等
  • R的高级绘图系统:
  1. grid绘图系统:基于grid绘图系统开发的ggplot2
  2. lattic绘图系统

R的基础知识介绍

> a <- 1e5  #1乘以10的5次方
> a
[1] 1e+05

> a <- 1e-10 #1乘以10的-10次方
> a
[1] 1e-10

> a <- exp(1) #自然对数的底数(e的几次方的意思)= e的1次方
> a
[1] 2.718282

理解下这两个的不同以及length.out参数

> seq(from=1,to=100,by=5)
 [1]  1  6 11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 96
> seq(from=1,to=100,length.out = 20)
 [1]   1.000000   6.210526  11.421053  16.631579  21.842105
 [6]  27.052632  32.263158  37.473684  42.684211  47.894737
[11]  53.105263  58.315789  63.526316  68.736842  73.947368
[16]  79.157895  84.368421  89.578947  94.789474 100.000000

理解下each参数

> rep(c(1:10),5)
 [1]  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10
[21]  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10
[41]  1  2  3  4  5  6  7  8  9 10
> rep(c(1:10),each=5)
 [1]  1  1  1  1  1  2  2  2  2  2  3  3  3  3  3  4  4  4  4  4
[21]  5  5  5  5  5  6  6  6  6  6  7  7  7  7  7  8  8  8  8  8
[41]  9  9  9  9  9 10 10 10 10 10

d[-1] 删掉d中的第一个元素
d=[c(-1,-3,-5)] 删掉d中的第1/3/5个元素
a = c(1:100) 我想把a里面的奇数位全都删掉
a[-1 * seq(1,100,2)] 把这个代码拆开看下

> a = c(1:100)
> -1 * seq(1,100,2)
 [1]  -1  -3  -5  -7  -9 -11 -13 -15 -17 -19 -21 -23 -25 -27 -29
[16] -31 -33 -35 -37 -39 -41 -43 -45 -47 -49 -51 -53 -55 -57 -59
[31] -61 -63 -65 -67 -69 -71 -73 -75 -77 -79 -81 -83 -85 -87 -89
[46] -91 -93 -95 -97 -99
> a[-1 * seq(1,100,2)]
 [1]   2   4   6   8  10  12  14  16  18  20  22  24  26  28  30
[16]  32  34  36  38  40  42  44  46  48  50  52  54  56  58  60
[31]  62  64  66  68  70  72  74  76  78  80  82  84  86  88  90
[46]  92  94  96  98 100

R中的 if 判断和 for 循环

d <- seq(from=1,to=100,by=5)
for(i in d){
  if(i >= 50){
    print(i)
  }
}

d <-seq(from=1,to=100,by=5) :1-100每隔5个选一个,选出来的所有数赋值给d;for循环加if循环的意思是:d里面的数,如果大于等于50就打印出来

自己创建函数

我们自己创建一个名为my_hj的函数,功能是:给定一个数x,当这个数小于等于0时,返回-x,否则就返回x

my_hj <- function(x){
  if(x <= 0){
    return(-x)
  }else{
    return(x)
  }
}

创建矩阵:matrix

g <- matrix(c(1:100), nrow = 10)  #将向量c(1:100),生成一个10行10列的矩阵
h <- matrix(c(1:15),5,5)  #将向量c(1:15)生成一个5行5列的矩阵,按说一共15个数不能生成5行5列的矩阵但R里并不会报错。数不够了就使用循环补齐,不报错给我们带来很大的困扰

#取矩阵的上三角和下三角
g_up <- upper.tri(g,diag = T)  #上三角 
#参数diag = T就是保留对角线的意思
g_low <- lower.tri(g)  # 下三角
t(g) # 转置
det(h) # 求方阵行列式
eigen(g) #特征值

创建data.frame

gene_id = c(1:100)
gene_fpkm = rnorm(100,10,5)
gene_fpkm = abs(gene_fpkm)  # fpkm基因表达量
sample_id = c(rep(1,50),rep(2,50))  # 50个1,,50个2
sample_id2 = round(runif(100,1,10)) # 从1-10中随机取100个呈均匀分布的数,再用函数round()把这些数变成整数。
gene_table = data.frame(gene_id,
                        gene_fpkm,
                        sample_id,
                        sample_id2)
colnames(gene_table)   # 查看数据框的列名
rownames(gene_table)  # 查看数据框的行名
dim(gene_table)  #查看数据框的维度,也就是数据框有多少行多少列

对数据框进行画图操作

这里有一个小重点,我想统计数据框gene_tablesample_id sample_id2中各个样本出现的次数,怎么办?用table()函数

> table(gene_table$sample_id)

 1  2 
50 50 
# sample_id这一列1出现了50次,2出现了50次

> table(gene_table$sample_id2)

 1  2  3  4  5  6  7  8  9 10 
 8 15 15 14  7 11  8 10  6  6 
#同理也给出了1-10个出现了多少次

只看这些数字没意思,我想把每个样本出现的次数画成条形统计图:用barplot()函数

barplot(table(gene_table$sample_id))
barplot(table(gene_table$sample_id2))

我想看下基因表达量(gene fpkm)的分布(distribution)情况:用hist()函数画一个直方图

hist(gene_table$gene_fpkm) ##提取数据框中基因表达量(gene_fpkm)那一列
基因表达量的分布
到现在我们都是在用这些包的基础用法,下面我们加一些参数画出稍微复杂一点的图
第一版,加上平均值线和中位数线
mean(gene_table$gene_fpkm) #求基因fpkm的平均数
median(gene_table$gene_fpkm) #求基因fpkm的中位数
abline(v = mean(gene_table$gene_fpkm)) #把平均值这条线竖直的加到图上
abline(v = median(gene_table$gene_fpkm)) #把中位数也加到图上
第一版图
第二版,对图进行精修:图加上颜色,线也加上颜色并加粗变成虚线
网址:https://coolors.co/browser/latest/1
hist(gene_table$gene_fpkm,col = "#C09BD8",border = F)  #柱体加上颜色去掉边框
abline(v = mean(gene_table$gene_fpkm),col="#001DFF",lwd=3,lty=3)  #平均值线加上颜色,加粗(lwd),变成虚线(lty)

各参数解释:

  • border = F #去掉柱体上的边框
  • col = "#C09BD8" #颜色参数,给一个色号"#C09BD8"
  • lwd 加粗
  • lty 变成虚线

推荐一个选颜色的网站:https://coolors.co/browser/latest/1
去这个网站选一个自己喜欢的配色,点 view 查看该配色里面的几个颜色,复制你喜欢的颜色的编号,比如我这里柱体选的是我喜欢的浅紫色#C09BD8,平均值线选的是显眼的蓝色#001DFF

第二版精修图
第三版,加一条拟合线
hist(gene_table$gene_fpkm,col = "#C09BD8",border = F) #柱体加上颜色去掉边框
abline(v = mean(gene_table$gene_fpkm),col="#001DFF",lwd=3,lty=3) #平均值线加上颜色,加粗(lwd),变成虚线(lty)
par(new=T) #下张图在上一张图的基础上继续画
plot(density(gene_table$gene_fpkm),xaxt="n",yaxt="n",bty="n") #xaxt="n"/yaxt="n" xy坐标去掉,bty="n"边框去掉

参数解释:

  • abline()函数是在加线,v参数意思是说,加一条竖直线(vertical curve)
  • xaxt="n"/yaxt="n"xy坐标去掉,bty="n"边框去掉

画拟合线:

plot(density(gene_table$gene_fpkm))
拟合线 第三版,加拟合线

用基础绘图系统作图

绘制RNA-Seq基因表达量的散点图

volcano plot 火山图

火山图(Volcano Plot)是大家在文章里面经常看到的图,特别是在展示差异表达的基因时,常常出现在芯片、测序等组学检测技术的结果中,与热图等等常一起出现。
我们看看火山图长什么样子

百度图片搜volcano plot随便找了一张

下面我们就说明下如何看懂和绘制火山图。
标准的火山图常用于展示显著差异表达的基因,这里有两个关键词:显著是指P<0.05,差异表达一般我们按照Fold Change(倍数变化)>=2.0作为标准。
当我们拿到基因表达的P值和倍数后,为了用火山图展示结果,一般需要把倍数进行Log2的转化,比如某基因在实验组表达水平是对照组的4倍,log2(4)=2,同样的如果是1/4,也就是0.25,转换后的结果就是-2。
同样的道理,对P值进行-log10的转化,-log10(0.05)约等于1.30103,由于P值越小表示越显著,所以我们进行-log10(P value)转化后,转化值越大表示差异约显著,比如-log10(0.001)=3 > -log10(0.01)=2 > -log10(0.05)=1.30。
比如:
RNA-Seq测序时会有一个对照组的样本,一个处理组的样本,处理组和对照组相比较,会发现某个基因表达量或变高了或变低了。用处理组的表达量除以对照组就会得到一个倍数,这个倍数就叫做Fold Change ,差异表达一般我们按照Fold Change(倍数变化)>=2.0作为标准。为了用火山图更直观的显示结果我们要将横轴的Fold Change进行Log2处理,对纵轴的P值进行-log10的处理
log2(fold change)=0也就是图中间的位置,就相当于:处理组fpkm/对照组fpkm=1也就是说对照组和处理组之间的基因表达量是没有变化的。所以说出现在中间位置的点是表达量没有发生变化的基因,越往右边走就说明这些基因的表达量在处理组比对照组大,越往左走说明相对于对照组来说,处理组基因表达量变低。
在该图中,不显著的被设置成了粉色,显著的设置成了绿色。
Y轴是cuffdiff算的统计检验的P-value取了log10之后又取了一个负数,所以说-log10(P-value)越高代表着越显著,也就是说越往左上角和越往右上角的点是越显著的点。
火山图(Volcano Plot)是在做RNA-Seq分析的时候特别常用的一张图,因为它可以非常清晰的展示出哪些点是我们真正想要的显著性的点。
我们用来自于cuffdiff的输出文件gene_exp.diff画图,读入R:

volcano_plot = read.table(file="gene_exp.diff文件路径",header = TRUE)

注:读表用read.table函数
读入后,看看这个表张什么样子:

gene_exp.diff
可以看到:
  • 第一列、第二列都是基因名;
  • 第三列gene没写因为已经给了gene_id了;
  • 第四列locus是该基因在基因组上位置;
  • sample_1 sample_2就是取了个名字不用注意;
  • status那列显示的是,是否通过检验,显示OK的就是进行检验的,显示NOTEST就说明表达量特别特别低不进行检验;
  • p-value是做统计检验看看是否显著;
  • q-value是对p-value进行修正,最终以修正过的q-value为准,q-value<0.05才认为是显著的。

开始做火山图

先准备X轴及log2(fold change)

前面我们交代过了fold change=处理组fpkm/对照组fpkm
所以:

control_FPKM = volcano_plot$value_2
treat_FPKM = volcano_plot$value_1 
log2_foldchange = log2(treat_FPKM / control_FPKM)

我们可以看到,我们这里重新赋值求的log2_foldchange,原来的表格里有log2.fold_change这一列,为什么还要重新求呢?因为原来表格里求反了求的是:对照组fpkm/处理组fpkm,我们这里改正下。
我们看下计算好的log2_foldchange里面有很多的InfNaN,这是怎么回事?

> log2_foldchange
   [1]  0.3382477530 -0.4605051707 -0.0964487964 -0.0249217524  0.7042396503           NaN           Inf
   [8] -0.9036291456           NaN           NaN  0.4595512248 -0.2834890803 -0.1511719401          -Inf
  [15]           Inf           NaN           NaN           NaN  0.4606267258  0.1979788743  0.1460252515
  [22] -0.3713316852 -0.9538916822 -0.2045978915 -0.2938850235  0.2349539904  0.3777615202  0.4461331227
  [29] -0.1815706646 -1.0014515800           Inf -0.3531645487  0.1261928713  2.8901005944          -Inf
  [36] -0.0146452448 -0.8561317308 -0.1874764797  0.5134263089 -0.0978467676 -0.3708083209  0.8298739826

这是因为有些control组基因的FPKM是0,除数是零,是不符合数学规则的,所以会出现这种情况。
我们需要处理下:
control_FPKM == 0 判断一下control_FPKM的值是否等于0
control_FPKM == 0即除数是0的全都筛选出来

> log2_foldchange[control_FPKM == 0 ]
   [1] NaN Inf NaN NaN Inf NaN NaN NaN Inf Inf NaN NaN NaN Inf NaN NaN NaN NaN Inf Inf NaN NaN NaN NaN NaN NaN
  [27] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN Inf Inf NaN NaN NaN NaN
  [53] NaN NaN Inf NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN Inf NaN NaN NaN Inf NaN

我们可以看到经过log2_foldchange[control_FPKM == 0 ] 操作,把log2_foldchange里的所有 NaN / Inf 选出来了,之后我们把筛选出的 NaN / Inf 强行赋值为0

log2_foldchange[control_FPKM == 0 ] = 0 

查看新生成的log2_foldchange发现还有一些-Inf

> log2_foldchange
   [1]  0.3382477530 -0.4605051707 -0.0964487964 -0.0249217524  0.7042396503  0.0000000000  0.0000000000
   [8] -0.9036291456  0.0000000000  0.0000000000  0.4595512248 -0.2834890803 -0.1511719401          -Inf
  [15]  0.0000000000  0.0000000000  0.0000000000  0.0000000000  0.4606267258  0.1979788743  0.1460252515
  [22] -0.3713316852 -0.9538916822 -0.2045978915 -0.2938850235  0.2349539904  0.3777615202  0.4461331227
  [29] -0.1815706646 -1.0014515800  0.0000000000 -0.3531645487  0.1261928713  2.8901005944          -Inf
  [36] -0.0146452448 -0.8561317308 -0.1874764797  0.5134263089 -0.0978467676 -0.3708083209  0.8298739826
  [43]  1.8493688172  0.0889640074  0.3113824129  0.0000000000  1.1102702040 -0.2531163954  0.0957645721
  [50] -0.6889365003  0.0000000000  0.0000000000  0.4828096335          -Inf -0.0058634211 -0.1436717399
  [57] -0.2029754567  1.9504284898 -0.2426462728 -0.3974959807  0.7211822730 -0.2428425868  0.5005116514
  [64]  0.3829284571 -0.0027983348  0.2696816982 -0.1575258190 -0.5746841946 -0.2751886483  0.0894799450
  [71] -0.1773208535 -1.6980255555  0.0000000000  0.5061142478  0.3639842228          -Inf -0.1930460419

这是因为当treat_FPKM=0即分母是0时,求log2就等于负无穷-Inf
所以,同理,我们把treat_FPKM也处理一下,把log2_foldchange里所有的-Inf赋值为0,这时我们就把画图过程中的异常点过滤掉了。

log2_foldchange[treat_FPKM == 0 ] = 0

此时log2_foldchange里所有的数都是有效数字了,原来NaN / Inf的地方都变成0啦~

现在我们准备Y轴。

前面我们也交代Y轴是p-value取log10,这时得到的是一个负值,我们再把它变成整数所以乘以-1即我们的Y轴。

log10_p_value = log10(cuffdiff_result$p_value) * -1

此时,X,Y轴都准备好了,我们开始作图

plot(x=log2_foldchange,y=log10_p_value)
volcano plot 初稿

火山图的精修

这张草图很明显不符合paper的要求,有许多需要改进的点:

  • 坐标轴需要修改,一般X轴我们只要[-4,4]之间的点;
  • 在0的地方有很多基点,也就是p-value=0的点,这些点一般来说是不画的;
  • p-value=0.05这根辅助线添加上;
  • 图上的颜色太难看了,我们给定它一个颜色,显著的用什么颜色,不显著的用什么颜色。
    下面我们一个一个修改

过滤坐标

下面我们对-log10(p-value)进行筛选,把等于0的点过虑掉(即不要最下面那一横排y=0的点)。

log10_p_value.filter = log10_p_value[log10_p_value >= 0.001]
log2_foldchange.filter = log2_foldchange[log10_p_value >= 0.001]

用过滤后的X/Y轴再次画图,使用xlim/ylim参数限定X/Y坐标轴,X轴取[-4,4];Y轴取[0,4]

plot(x=log2_foldchange.filter,y=log10_p_value.filter,xlim=c(-4,4),ylim=c(0,4))

此时从图中可以看出来,-log10(p-value)=0的点去掉了

显著/不显著加不同的颜色

接下来,给图加上漂亮的颜色

思路:所有点先改成蓝色 -> 找出显著性的点 -> 将显著性的点变成红色

所有点变成蓝色
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
     )
修改版

参数解释:

  • col=rgb(0,0,1,0.1)这是一个颜色参数,r是红色;g是绿色;b是蓝色。所以(0,0,1,0.1)代表,没有红色、没有绿色、全是蓝色,透明度是0.1
  • pch=16图上的点是空心的,这个参数就是让它变成实心的点。
找出显著性的点

首先要知道被我们过滤完之后,一共还有多少个点:(12769个)

> length(log2_foldchange.filter)
[1] 12769

创建一个储存颜色的向量,向量里存放着重复了12769次的颜色参数——rgb(0,0,1,0.1)

> col_vector = rep(rgb(0,0,1,0.1),length(log2_foldchange.filter))
> col_vector
   [1] "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A"
   [8] "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A"
  [15] "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A" "#0000FF1A"

我们可以看到rgb(0,0,1,0.1)在运行后会变成#0000FF1A,前两位数字00就代表r(红色);再往后两位00代表g(绿色);再往后FF代表b()蓝色;1A就表示透明度。这就是把rgb(0,0,1,0.1)参数变成了一个颜色编号,无所谓了,对我们作图没有影响。

满足什么样条件的点通常被认为是显著性的点呢?

  • p_value <= 0.05 ;即Y轴:log10_p_value.filter >= 1.30103:-log10(0.05)=1.30103
  • fold change 大于2,或小于二分之一 ,即X轴的绝对值大于等于1:abs(log2_foldchange.filter) >= 1
  • 对照组和实验组的FPKM > 0;
    通知满足这三个条件才能判断为显著,现在我们把同时满足这三个条件的基因筛选出来并写进一个向量select_sign_vector
select_sign_vector = (volcano_plot$value_1 > 0 ) & (volcano_plot$value_2 >0) & (abs(log2_foldchange) >= 1) & (volcano_plot$value_1 >= 1 | volcano_plot$value_2 >= 1) & (volcano_plot$p_value <= 0.05)

abs(log2_foldchange) >= 1意思是log2_foldchange的绝对值大于等于1,那也就是说:我们只取,相对于对照组,处理组要么升高两倍,要么缩小为原来的二分之一的这些点;
volcano_plot$value_1 >= 1 | volcano_plot$value_2 >= 1的意思是value_1value_2有一个大于等于1就可以了。
把具有显著差异的点筛选出来后,我们使用table()函数看下显著的基因有多少个

> table(select_sign_vector)
select_sign_vector
FALSE  TRUE 
22774    83

我们的目的是将这些显著差异的点赋值为红色rgb(1,0,0)

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] #上面对X/Y轴进行过滤了,所以这里也需要按照童谣标准过滤一下。

col_vector[select_sign_vector.filter] = rgb(1,0,0)  #将显著差异的点所对应的颜色参数换成红色 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
)
将显著差异的点标上红色
加一条p-value的辅助线
abline(h=-1*log10(0.05),lwd=3,lty=3,col="#4C5B61")
abline(v=log2(2) ,lwd=3,lty=3,col="#4C5B61")
abline(v=log2(1/2) ,lwd=3,lty=3,col="#4C5B61")

解释各参数:

  • h是加一条水平线(horizontal line);
  • v是加一条竖直线(vertical curve);
  • lwd线条加粗,且选3号样式;
  • lty线条做成虚线;
  • col给线一个色号"#4C5B61"
    精修图
    基础绘图系统做火山图——以上代码整合:
rm(list=ls())  #好习惯要养成,先清空下环境变量
volcano_plot = read.table(file="文件路径",header = TRUE)
log2_foldchange = log2(volcano_plot$value_1 / volcano_plot$value_2) #X轴

log2_foldchange[volcano_plot$value_2 == 0 ] = 0  #筛选X轴
log2_foldchange[volcano_plot$value_1 == 0 ] = 0 #筛选X轴

log10_p_value = log10(volcano_plot$p_value) * -1 #准备Y轴
#X/Y轴都准备好了,画草图
plot(x=log2_foldchange,y=log10_p_value)

##修正X/Y轴,去掉-log10(p-value)=0的点
log10_p_value_qc = log10_p_value[log10_p_value >= 0.001]
log2_foldchange_qc = log2_foldchange[log10_p_value >= 0.001]

plot(x=log2_foldchange_qc, y=log10_p_value_qc, xlim=c(-4,4), ylim=c(0,4))

#思路:所有点先改成灰色  -> 找出显著性的点 -> 将显著性的点变成蓝色

#变颜色-灰色
plot(x=log2_foldchange_qc, y=log10_p_value_qc,
     xlim=c(-4,4),ylim=c(0,4),
     col="#BCBABE",pch=16
     )

#把显著性的点写进一个向量`sign_point`
sign_point = (abs(log2_foldchange_qc) >= 1) & (log10_p_value_qc >= 1.30103)

#将所有点的颜色编号`#BCBABE`放入到向量`col_point`里
col_point = rep("#BCBABE", length(log2_foldchange_qc))

#把筛选出来的显著性的点`sign_point`变成蓝色
col_point[sign_point] = "#1B2CC1"

plot(x=log2_foldchange_qc, y=log10_p_value_qc,
     xlim=c(-4,4),ylim=c(0,4),
     col=col_point,pch=16
     )

#加辅助线

abline(h=-1*log10(0.05),lwd=3,lty=3,col="#4C5B61")
abline(v=log2(2) ,lwd=3,lty=3,col="#4C5B61")
abline(v=log2(1/2) ,lwd=3,lty=3,col="#4C5B61")

Heatmap plot 热图—Hi-C可视化

这是一张很典型的Hi-C可视化的图——横坐标和纵坐标都是染色体位置坐标;蓝色虚线三角形内为TAD结构。文献里报道这张图是用ggplot2做的,我们今天的任务就是,用基础绘图系统,做出这样的图。

7号染色体不同分辨率的染色质相互作用热图

首先明确何为 heatmap plot 热图 ?heatmap: marix -> rect + color
heatmap拆开看就是一个矩阵变成很多小方块儿加颜色
heatmap说白了就是一个矩阵,矩阵里有很多数,数值大对应一个颜色,数值小的对应一个颜色。

rect()函数画小方块儿

先创建画布

dev.off()  #关闭之前画的图
plot(x=c(1:10),y=c(1:10)) #创建一个画布
创建的画布
去掉画布上的空心点,加参数type="n"即可:
plot(x=c(1:10),y=c(1:10),type="n")
去掉画布上的点
使用rect()函数在画布上画小方块儿

我们运行?rect()查看下它的help文档里面的Usage部分:

rect(xleft, ybottom, xright, ytop, density = NULL, angle = 45,
     col = NA, border = NULL, lty = par("lty"), lwd = par("lwd"),
     ...)

发现有好多参数,其中参数 xleft, ybottom, xright, ytop 决定了小方块的位置和大小,坐标(xleft,ybottom)是方块左下角那个点的坐标;坐标(xright,ytop)是方块右上角那个点的坐标。两个对角线的坐标不但决定了小方框在画布的位置,还决定了小方框的大小。
举例子说明:
给定对角线坐标画小方格,并添加颜色

rect(xleft = 0,ybottom = 0,xright = 1,ytop = 1,col=rgb(1,0,0,0.5))
rect(xleft = 5,ybottom = 5,xright = 6,ytop = 6,col="red")

得到图形


理解xleft, ybottom, xright, ytop参数

创建画多个小方块儿的函数

现在我们创建一个函数,该函数的功能是:给它一个矩阵它就可以输出一张图。
要完成以上要求,首先需要确定所有的 xleft, ybottom, xright, ytop ,也就是我们要确定每个小方块的对角线坐标;其次我们需要确定color即每个格子的颜色是什么。
先准备一个矩阵:

input_matrix = matrix(c(1:36),6,6)
6行6列的矩阵
确定所有小方块儿的坐标

第一步,用dim()函数确定图画多大的:6*6大小的~

# set image size 
> x_size = dim(input_matrix)[1]
> y_size = dim(input_matrix)[2]
> x_size
[1] 6
> y_size
[1] 6

我们从0,0 向 6,6方向找规律
先画,第1行,再画第2行,就用这两行的对角线坐标找到规律


帮助理解取点

第一行小方块左下角的坐标即(xleft, ybottom)分别是:
(0,0);(1,0);(2,0);(3,0);(4,0);(5,0)
所以:

  • xleft:0,1,2,3,4,5
  • ybottom:0,0,0,0,0,0
    第一行小方块右上角的坐标即(xright, ytop)分别是:
    (1,1);(2,1);(3,1);(4,1);(5,1);(6,1)
    所以:
  • xright:1,2,3,4,5,6
  • ytop:1,1,1,1,1,1
    综上所述,第一行的xleft, ybottom, xright, ytop分别为:
  • 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
    同理,我们可以得到第二行的xleft, ybottom, xright, ytop分别为:
  • 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
    根据这两行,我们找到xleft, ybottom, xright, ytop规律:
my_xleft = rep(c(0:(x_size-1)),each = x_size)  #获取每一行的xleft
my_xright = my_xleft + 1
my_ybottom = rep(c((y_size-1):0),y_size) #获取每一行的ybottom
my_ytop = my_ybottom + 1

我们确定了所有格子的xleft, ybottom, xright, ytop可以根据格子坐标画格子了

plot(x=c(0:x_size),y=c(0:y_size),type="n")  #准备一个6行6列的画布
rect(xleft = my_xleft,ybottom = my_ybottom,xright = my_xright,ytop = my_ytop) #画布上画上6*6个方格
画上方格子

精修小方格——加颜色、去边框

画最简单的颜色,最小值是白色,最大值是红色,中间呈线性变化。
需要得到颜色矩阵,所谓颜色矩阵就是矩阵里放着每个格子的色号,为了得到颜色矩阵我们首先需要得到颜色浓度矩阵input_matrix.rate

获取颜色浓度矩阵

我们设定最大值是红色,那用input_matrix中每一个格子除以矩阵里的最大值,得到颜色浓度矩阵input_matrix.rate,这样操作后颜色浓度矩阵里数越大红色越深。

mat.max = max(input_matrix)
input_matrix.rate = input_matrix / mat.max #这步得到的就是颜色的浓度
col.mat = rgb(1,0,0,input_matrix.rate) #rgb(1,0,0)代表红色,一般再后面的参数放的就是透明度,这里就是颜色的浓度。

此时就得到了颜色矩阵col.mat

> col.mat
 [1] "#FF000007" "#FF00000E" "#FF000015" "#FF00001C" "#FF000023" "#FF00002B" "#FF000032" "#FF000039"
 [9] "#FF000040" "#FF000047" "#FF00004E" "#FF000055" "#FF00005C" "#FF000063" "#FF00006A" "#FF000071"
[17] "#FF000078" "#FF000080" "#FF000087" "#FF00008E" "#FF000095" "#FF00009C" "#FF0000A3" "#FF0000AA"
[25] "#FF0000B1" "#FF0000B8" "#FF0000BF" "#FF0000C6" "#FF0000CD" "#FF0000D5" "#FF0000DC" "#FF0000E3"
[33] "#FF0000EA" "#FF0000F1" "#FF0000F8" "#FF0000FF"

有了颜色矩阵就可以把每个小方格加上颜色了。

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)
加上颜色了

图虽然画出来了,但是这张热图问题很大的:

  • 所有格子的边框不要
  • 留白不要

去除外边框和内边框

rect()函数中加上参数border = F可以去掉每个小方格的黑框;
plot()函数中加上参数frame.plot = F,xaxt="n",yaxt="n",xlab="",ylab=""可以去掉一些我们不想要的留白和坐标,每个代表什么含义呢,看图

解释各参数
plot(x=c(0:x_size),y=c(0:y_size),type="n",frame.plot = F,xaxt="n",yaxt="n",xlab="",ylab="")
rect(xleft = my_xleft,ybottom = my_ybottom,xright = my_xright,ytop = my_ytop,col=col.mat,border = F)

此时,不想要的东西都删掉了

合格图
这时热图就绘制完成了,之后用真实数据画热图的时候只需要把数据文件读进去就可以啦~上面我们是一条一条代码拆开讲的,下面我们把它整合到一起,之后用真实数据画图的时候方便改动。

基础绘图系统热图代码整合:

# 把该清的清一清
rm(list=ls())
dev.off()

#读入数据文件
input_matrix = matrix(abs(rnorm(10000,1000,500)),100,100) #这里只是随机生成了一个矩阵,换成自己的数据文件即可。

# set image size 
x_size = dim(input_matrix)[1]
y_size = dim(input_matrix)[2]

# 确定坐标
my_xleft = rep(c(0:(x_size-1)),each = x_size)  #获取每一行的xleft
my_xright = my_xleft + 1
my_ybottom = rep(c((y_size-1):0),y_size) #获取每一行的ybottom
my_ytop = my_ybottom + 1

# 精修
## 得到浓度矩阵
mat.max = max(input_matrix)
input_matrix.rate = input_matrix / mat.max #这步得到的就是颜色的浓度
col.mat = rgb(1,0,0,input_matrix.rate) #rgb(1,0,0)代表红色,一般再后面的参数放的就是透明度,这里就是颜色的浓度。


## 画图——加颜色、去外边框、内边框、留白
plot(x=c(0:x_size),y=c(0:y_size),type="n",frame.plot = F,xaxt="n",yaxt="n",xlab="",ylab="")
rect(xleft = my_xleft,ybottom = my_ybottom,xright = my_xright,ytop = my_ytop,col=col.mat,border = F)


#圈出某个特定的区域
segments(x0 = 20,x1 = 20,y0 = 50,y1 = 80,lwd=3)
segments(x0 = 20,x1 = 50,y0 = 80,y1 = 80,lwd=3)
segments(x0 = 20,x1 = 50,y0 = 50,y1 = 50,lwd=3)
segments(x0 = 50,x1 = 50,y0 = 50,y1 = 80,lwd=3)
segments(x0 = 20,x1 = 50,y0 = 50,y1 = 80,lwd=3)

使用真实数据在画图的时候,除了把真实数据读入,还在得到颜色矩阵的地方进行了调整~
将之前的代码:

## 得到浓度矩阵
mat.max = max(input_matrix)  
input_matrix.rate = input_matrix / mat.max #这步得到的就是颜色的浓度
col.mat = rgb(1,0,0,input_matrix.rate) #rgb(1,0,0)代表红色,一般再后面的参数放的就是透明度,这里就是颜色的浓度。

微调为:

## 得到浓度矩阵
mat.max = quantile(input_matrix,prob=0.9) #下面解释
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))) #rgb(1,0,0)代表红色,一般再后面的参数放的就是透明度,这里就是颜色的浓度。

这样就可以成功出图啦~~
参数解释:

  • quantile(x,probs)求分位数。其中 x 为待求分位数的数值型向量, probs 为一个由[0,1]之间的概率值组成的数值向量。例子,# 求 x 的 30%和 84%分位点:y <- quantile(x, c(0.3,0.84))

这里为什么进行这样的修改?

因为如果不改,画出的图特别特别的淡,大部分都是白色的区域,造成这样的原因主要是,最大值与平均值之间差距太大了,我们得到颜色浓度的时候,直接用矩阵中的每一个数除以最大值,就会造成把颜色浓度拉的太低的情况。
所以我们用quantile(x,probs)函数,去矩阵中90%分位点的数,可以这样理解:我们把矩阵中所有的数从大到小进行排序,去掉最大的前10%,求出来的是剩下那90%里最大的数。这样就会避免出现那样的问题。
input_matrix.rate[input_matrix.rate>1] = 1 上面是取90%分位点的数,所以必须加上这步操作。颜色浓度矩阵是这样算的:矩阵中每个数除以90%分位点的数input_matrix / mat.max得到一个从0-1的分数。所以当遇到最大的前10%的数时,这个比率就会大于1,这样不行,所以我们把所有大于1的强行赋值为1。

真实的Hi-C数据做热图

上面已经把代码微调好了,我们这里直接搬过来用
总结下,就是在上面的基础上做了两个改动:

  • 得到颜色浓度矩阵的时候改了,为什么那样改,上面提到了。
  • 在画图代码前加了一条命令:png(file="./hic_png",width = 1000,height = 1000)这样生成的图就存在本地电脑上了,不在R的画布上,大数据作图的时候这样操作可以节省很多时间。
# 把该清的清一清
rm(list=ls())
dev.off()

#读入数据文件
hic_mat.raw = read.table(file="./20171203-Live-R_partII/data_file/chr_16_100000_MAPQ20.txt",sep = ",",header = F)
input_matrix = as.matrix(hic_mat.raw)

# set image size 
x_size = dim(input_matrix)[1]
y_size = dim(input_matrix)[2]

# 确定坐标
my_xleft = rep(c(0:(x_size-1)),each = x_size)  #获取每一行的xleft
my_xright = my_xleft + 1
my_ybottom = rep(c((y_size-1):0),y_size) #获取每一行的ybottom
my_ytop = my_ybottom + 1

# 精修
## 得到浓度矩阵
mat.max = quantile(input_matrix,prob=0.9)
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))) #rgb(1,0,0)代表红色,一般再后面的参数放的就是透明度,这里就是颜色的浓度。


#大数据画图处理方式
png(file="./hic_png",width = 1000,height = 1000)
plot(x=c(0:x_size),y=c(0:y_size),type="n",frame.plot = F,xaxt="n",yaxt="n",xlab="",ylab="")
rect(xleft = my_xleft,ybottom = my_ybottom,xright = my_xright,ytop = my_ytop,col=col.mat,border = NA)
dev.off()

#圈出某个特定的区域
segments(x0 = 20,x1 = 20,y0 = 50,y1 = 80,lwd=3)
segments(x0 = 20,x1 = 50,y0 = 80,y1 = 80,lwd=3)
segments(x0 = 20,x1 = 50,y0 = 50,y1 = 50,lwd=3)
segments(x0 = 50,x1 = 50,y0 = 50,y1 = 80,lwd=3)
segments(x0 = 20,x1 = 50,y0 = 50,y1 = 80,lwd=3)

相关文章

  • R入门与基础绘图系统

    这是我听B站鲮鱼不会飞视频(R入门与基础绘图系统)里的笔记哦~加个目录吧,因为实在太长了,怕找不到重点:一、前面是...

  • R入门1 R入门与基础绘图系统学习笔记

    我的R学习领路人是北大生信在读博士孟浩巍,学习R是从参加他的知乎live课程《R入门与基础绘图系统》开始的,从第一...

  • R绘图函数

    R语言四大作图系统: 基础绘图系统 lattice包 ggplot2包 grid包 R绘图分类: 高级绘图(搭好框...

  • 2.qplot入门

    qplot入门 qplot是quick plot的缩写,与R系统的基本绘图函数plot类似。 diamonds数据...

  • R语言入门与基础绘图系统 1

    1. 什么是R语言? R语言是一种自由软件编程语言与操作环境,主要用于统计分析、绘图和数据挖掘。 R语言是从S语言...

  • Day—4

    R 编程语言,是进行统计计算和绘图的环境 R语言的基础绘图系统主要由基础包graphics提供 Rstudio 图...

  • R语言入门--第十六节(ggplot2绘图)

    之前学习的绘图方法是基于R的基础绘图系统。在R中一共有四种作图系统,分别为base(之前学的)、grid、latt...

  • R语言可视化及作图1--基础绘图(par函数,散点图,盒形图,条

    R语言绘图系统基础绘图包 ⚠️Lattice (语法复杂)ggplot2家族 ⚠️其他,比如:sjplot; pl...

  • 数据挖掘20210111学习笔记

    R语言作图 低级绘图函数建立在高级绘图函数基础上,不能单独使用 ggplot2语法 1.入门级绘图模板2.映射-颜...

  • R|绘图-学习笔记(二)

    tags: 绘图 R的三大绘图系统 1. 基本绘图系统 (base plotting system) 绘图始于空白...

网友评论

    本文标题:R入门与基础绘图系统

    本文链接:https://www.haomeiwen.com/subject/zlpvgqtx.html