R语言打基础 — boxplot / heatmap

2019-05-05 本文已影响165人黄晶_id

这是我听B站—鲮鱼不会飞视频（R语言相关内容补充）里的笔记哦~

先读入我们的数据

cuffdiff_result = read.table(file = "hela_gene_exp.diff",header = T,sep = "\t")

参数解读：

使用read.table函数读入一个文件；
header = T：第一行设为行名；
sep = "\t"：以Tab（'\t'）为分隔符。

我们读入的文件这是一个用cuffdidd算出来的真实数据，介绍一下各列都代表什么意思：
第1—3列都是基因相关的信息，都是些基因名；
locus：那列是range，就是基因在染色体上的起始位置；
sample_1/sample_2：算的时候给它的一个名字，sample_1 是control组sample_2是treatment组；
value_1和value_1_2 ：是算出来的FPKM，value_1是基因在control组里算出的FPKM，value_2是基因在treatment组里算出的FPKM。
比如，图中第5行A2LD1基因在control组的FPKM是0.7208970(value_1)，在treatment组里的FPKM是1.1026200(value_2)，那该基因在实验组和对照组的fold_change= value_2/value_1 (用处理组的表达量除以对照组就会得到一个倍数，这个倍数就叫做Fold Change )，在作图的时候我们是要对fold_change取一个log2的；

> log2(1.1026200/0.7208970)
[1] 0.6130706

status：当计算过程中有些地方不是很准确的时候，软件不进行统计学检验，在这列会写一个NOTEST ；当是ok时，说明软件进行了统计学检验；
p_value/q_value ：软件的统计结果会给出一个p_value和一个修正之后的p_value也就是q_value ；
significant：默认情况p_value<0.05都是显著的，在实际画图过程中我们找显著性差异表达基因时筛选条件会更严格一点，否则只按照这个标准来选，选出来的差异基因太多了。

找显著的差异基因

筛选出significant那一列是yes的行

> cuffdiff_result$significant == "yes"
   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE 
  [16] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE 
  [31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE

将显著的信息赋值给一个新的向量filter_vector：

filter_vector = cuffdiff_result$significant == "yes"

将所有significant == "yes"的行放入一个新的表中:

cuffdiff_result.sign = cuffdiff_result[filter_vector,]

选出yes列

我们可以看到筛选出了5189个显著的基因，在实际操作过程中，我们不能仅仅以这一个条件为基准，一般有这几个条件：

差异基因筛选条件

p_value<0.05； #统计显著性；

log2(fold_change) > 1 或者 log2(fold_change) < -1 换成人话就是：比起对照组的FPKM必须升高或降低两倍以上； #要存在剧烈变化；

one of tow case fpkm >= 1 #control和treatment里面必须有一个FPKM大于等于1，用这个条件筛选出具有生物统计显著性的差异基因。

我们下面用代码筛选出符合上述三个条件的数据

filter_vector_p_value = cuffdiff_result$p_value < 0.05
filter_vector_fc = abs(cuffdiff_result$log2.fold_change.) > 1
filter_vector_fpkm = cuffdiff_result$value_1 > 1 | cuffdiff_result$value_2 > 1

当同时满足这三个条件时才是我们要找的差异基因，所以我们做一个交集：

filter_vector_all = filter_vector_p_value & filter_vector_fc & filter_vector_fpkm

这时我们再将这些真正显著的放入一个新表里cuffdiff_result.sign.real

cuffdiff_result.sign.real = cuffdiff_result[filter_vector_all,]

生物学意义上的差异基因

通过这种方法就能找到真正生物学上的差异基因，这种方法筛出来的才是比较靠谱的。

用筛出来的差异基因矩阵画图

plot boxplot

看下我们筛选出来的这些显著的基因分布情况（control组和treatment组都要看）

boxplot(cuffdiff_result.sign.real$value_1,cuffdiff_result.sign.real$value_2)

从图中我们可以看到，绝大部分的基因都是低表达的（符合生物体的特征），而且表达量高的特别高，低的又特别低 “贫富差距” 太大，为了让图能更直观的反应问题，我们需要进行一些处理，一般会对表达量取一个log2，但是有些样本在某个基因的表达量是0，所以我们一般会进行加1处理，这就是我们常用的画图方式——对Y轴加1再取个log2：

boxplot(log2(cuffdiff_result.sign.real$value_1+1),log2(cuffdiff_result.sign.real$value_2+1))

这个时候，我们就能比较直观和清楚的看到表达量的分布情况啦，横轴1和2分别是value_1和value_2；纵轴就是FPKM取过log2的，从图中我们可以看到，这些显著的基因处理之后，大部分都是下调的。

然后我们给它加个颜色让它更美观吧—左边是"orange"，右边是"pink"
添加col参数给图形加颜色

boxplot(log2(cuffdiff_result.sign.real$value_1+1),log2(cuffdiff_result.sign.real$value_2+1),col = c("orange","pink"))

heatmap—画热图

先安装加载画热图用的包—pheatmap

install.packages("pheatmap")
library(pheatmap)

pheatmap这个包要求你输入的是一个矩阵，其实value_1/value_2两列就够了,使用矩阵取子集的方法取出第8列和第9列：cuffdiff_result.sign.real[,c(8,9)]，然后用这两列作图：

pheatmap(cuffdiff_result.sign.real[,c(8,9)])

这张图还是那个问题，最大值和最小值还有平均值差距太大了，而且大部分都是低表达，把整张图的色调整体拉下来了，所以还是需要做一步处理—加1再取log2

pheatmap(log2(cuffdiff_result.sign.real[,c(8,9)]+1))

处理后再看图，这就是我们在论文里常看到的heatmap图了，每一行是一个基因，颜色就代表这个基因的表达量，我们可以看到这些基因都是我们挑出来的差异基因。value_1是control组，value_2是treatment组，比较两组的颜色变化我们可以看到处理后大部分的基因表达量下降了。

1556982776970.png
左边的聚类是对行进行聚类，也就是说，行与行之间表达量最接近的放在一起，所以聚类结果就会出现高表达的聚在一起，低表达的聚在一起，由低表达变成高表达的聚在一起。

注：有时候会在文章里看到heatmap有一个 z score是怎么回事？当有多个repeat的时候对于每一行进行了normaliz，normaliz过后的数值叫做 z score，因为我们这里只有两个样本做比较，所以没有办法算 z score。