这是我听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_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)