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的基础绘图系统：

基础绘图系统中的低级绘图函数：创建画布，点，线，多边形
基础绘图系统中的高级绘图函数：plot(), boxplot(), hist(), density() 等等

R的高级绘图系统：

grid绘图系统：基于grid绘图系统开发的ggplot2
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_table中sample_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里面有很多的Inf和NaN，这是怎么回事？

> 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_1或value_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入门与基础绘图系统

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

R的绘图系统

R的基础知识介绍

R中的 if 判断和 for 循环

自己创建函数

对数据框进行画图操作

用基础绘图系统作图

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

volcano plot 火山图

开始做火山图

先准备X轴及`log2(fold change)`

现在我们准备Y轴。

火山图的精修

过滤坐标

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

所有点变成蓝色

找出显著性的点

Heatmap plot 热图—Hi-C可视化

用`rect()`函数画小方块儿

使用`rect()`函数在画布上画小方块儿

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

确定所有小方块儿的坐标

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

获取颜色浓度矩阵

去除外边框和内边框

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

真实的Hi-C数据做热图

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

生信学习

R入门与基础绘图系统

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

R的绘图系统

R的基础知识介绍

R中的 if 判断和 for 循环

自己创建函数

对数据框进行画图操作

用基础绘图系统作图

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

volcano plot 火山图

开始做火山图

先准备X轴及log2(fold change)

现在我们准备Y轴。

火山图的精修

过滤坐标

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

所有点变成蓝色

找出显著性的点

Heatmap plot 热图—Hi-C可视化

用rect()函数画小方块儿

使用rect()函数在画布上画小方块儿

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

确定所有小方块儿的坐标

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

获取颜色浓度矩阵

去除外边框和内边框

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

真实的Hi-C数据做热图

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

生信学习

先准备X轴及`log2(fold change)`

用`rect()`函数画小方块儿

使用`rect()`函数在画布上画小方块儿