表观调控13张图之一证明基因干扰有效性
不知不觉进入了学徒培养第10周，前面该分享的各种NGS组学数据分析都已经发布在生信技能树公众号了，现在正式进入多组学整合阶段，其中本次系列分享的主角：表观调控，就是整合RNA-seq数据和表观数据，比如Chip-seq和ATAC-seq，主要是依托于文章 Global changes of H3K27me3 domains and Polycomb group protein distribution in the absence of recruiters Spps or Pho 文献解读在 在果蝇探索 PRC 复合物（逆向收费读文献 2019-18） 这一推文。
我把表观调控数据分析，拆分成为了13张图，分别录制为13个视频，即将免费发布在B站，这个期间我们的视频编辑师还在兢兢业业的奋斗，希望这13张图能带领大家学会表观调控数据分析的一般流程， 然后应用到自己的课题哈!
同时我也招募了4位视频审查员和细节知识补充员。接下来我们就连载第一位视频审查员的13个笔记：
关于视频审查员
第一位视频审查员大家也许并不陌生了，早在2018-08-29我发布给学徒的ATAC-seq数据实战（附上收费视频） 他就学完了全部课程内容，还写了笔记在简书，大家搜索二货潜，就可以看到。
背景
现代表观遗传学的定义为：不依赖于 DNA 序列改变的可遗传的变异 (Bonasio et al 2010)。这句定义揭示经典遗传和表观遗传的共性和个性，所谓共性：即可遗传的变异，这是遗传学的基础，只有将基因表达被改变的信息传递给下一代才具有遗传学和生物学意义；然而这两种遗传也有个性：即表观遗传学基因表达的改变不依赖 DNA 序列，经典遗传学建立在 DNA 序列的改变。
表观调控包括 DNA 甲基化、组蛋白修饰和转录因子的调节等 (Maleszewska and Kaminska 2013)。
组蛋白N末端的赖氨酸上除了可以进行多种乙酰化的修饰外也可以进行多种甲基化的修饰(Bhaumik et al 2007)。每个位点的甲基化均存在三种甲基化修饰形式，分别是单、二和三甲基化。不同位点的甲基化以及甲基化的程度直接影响了生物体内染色质的状态及基因的表达。
H3K27me3是发生在组蛋白 H3 的第 27 位赖氨酸上的三甲基化修饰，从整体染色质水平来讲, 大量的证据支持H3K27me3与异染色质相关；从个体的基因水平来讲,H3K27me3常常与转录抑制相关。H3K27me2/me3甲基化的加载是由PcG复合体催化的。如何精确地募集
Polycomb group（PcG）蛋白到其靶基因仍知之甚少。在果蝇中，PcG蛋白被募集到由多个DNA结合蛋白的结合位点组成的Polycomb反应元件（PRE）。为了了解如何将PcG蛋白募集到PRE并在其上维持，作者系统地研究了野生型的PcG结合，相关的H3K27me3和转录组以及三种PRE结合蛋白的突变体（ Pho、Spps、cg ）。仅
Pho结合位点不足以招募PcG。PRE由复杂的DNA结合蛋白的结合位点组成，包括Pho / Phol，Spps，Cg，GAF / Psq，Adf-1，Grh，Dsp1和Zeste / Fsh-S。然而这些蛋白质在PcG蛋白质募集中的作用尚不清楚。Spps是与PREs结合的Sp1 / KLF锌指蛋白家族的成员
第一张图：说明基因干扰的有效性
要强调的一点我们所画的图都是为了说明生物学问题，让人一目了然就能看出你想展示什么，切勿忘本。
当
PcG招募成员被中断后，基因表达是如何改变的。特别是，差异表达与H3K27me3修饰的改变有何关系？作者对野生型、Spps 和 Pho 突变体 3 龄幼虫进行了RNA-seq测序。
而本张图用来证实相应突变体中 Spps 和 Pho 表达确实降低了。
注：当我测基因突变体的RNA-seq时候，我们自然在获得表达量时候第一件事情是去查看，这些想对于的基因在突变体中是否真的被抑制或者不表达了。
当我们拿到转录组数据分析
featurecounts或者htseq等工具计算得到的基因的reads数目时候（这里我们先不考虑TPM或者reads的均一化之类哈），我们可以类似qPCR一样来展示每个样品中基因的表达量信息。那么怎么做呢？就看下面吧。本图仅仅是关注下面的RNA-seq数据。
rm(list = ls()) # 日常运行代码前清理一下环境
library(ggpubr)
library(stringr)
library(ggplot2)
library(cowplot)
自定义我的 `barplot` 主题
一般绘制某一类型图时候，当我们确定他大概所涉及的主题时，都会首先定义一个属于自己的主题，下次就可以直接用于此类图形。
my_bar_theme <- theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 16),
# 因为 x 轴标签要旋转 90°，所以这里用来旋转
axis.text.y = element_text(size = 16),
axis.title.y = element_text(size = 18,
face = "bold",),
legend.position = "none") +
theme(plot.title = element_text(size = 20, face = "bold", hjust = 0.5)) # 基因名要居中，这里用来居中。
数据读取
主要是RNA-seq数据上游分析后，featurecounts得到的表达矩阵：all.counts.id.txt
options(stringsAsFactors = F)
a <- read.table('all.counts.id.txt', header = T)
# 查看数据行列信息，大致看下多少基因
dim(a)
[1] 17714 16
数据清洗、可视化
cg <- a[a[,1] == 'pho',7:16]
dat <- data.frame(gene = as.numeric(cg),
sample = str_split(names(cg),'\\.',simplify = T)[,1],
group = str_split(names(cg),'_',simplify = T)[,1]
)
labels = c(paste0("WT", "_",1:3), paste0("PhoKO", "_", 1:3), paste0("SppsKO", "_", 1:4))
ggbarplot(dat,x = 'sample', y = 'gene',
color = 'black', fill = "group",
size = 1) +
scale_fill_manual(values = c(WT = "#9C4B25",
PhoKO = "#DE2C1C",
SppsKO = "#43A542")) +
scale_color_manual(values = "black") +
scale_x_discrete(limits = labels) +
labs(y = "Normalized read count",
x = "",
title = "Pho") +
my_bar_theme +
scale_y_continuous(expand = c(0, 0)) +
geom_text(aes(y = gene * 1.1, label = ""))
封装画图函数
对于会点点编程的人来说，肯定不会同意每画一个基因就重新运行一遍撒，所以我们将上面封装成函数
my_barplot <- function(gene){
cg <- a[a[,1] == gene, 7:16] # 提取候选的表达量对应的行和列
dat <- data.frame(gene = as.numeric(cg),
sample = str_split(names(cg),'\\.',simplify = T)[,1], # 这里将样品后面的 bam文件后缀去掉
group = str_split(names(cg),'_',simplify = T)[,1]
)
# x 轴标签的顺序，这里是按照原图的顺序来的
labels = c(paste0("WT", "_",1:3), paste0("PhoKO", "_", 1:3), paste0("SppsKO", "_", 1:4))
ggbarplot(dat, x = 'sample', y = 'gene',
color = 'black', fill = "group",
size = 1) +
scale_fill_manual(values = c(WT = "#9C4B25",
PhoKO = "#DE2C1C",
SppsKO = "#43A542")) +
scale_color_manual(values = "black") +
scale_x_discrete(limits = labels) +
labs(y = "Normalized read count",
x = "",
title = gene) +
scale_y_continuous(expand = c(0, 0)) +
geom_text(aes(y = gene * 1.1, label = "")) +
my_bar_theme
}
组合两个基因的图
Pho_plot <- my_barplot("pho")
Spps_plot <- my_barplot("Spps")
gene_exp_plot <- plot_grid(Pho_plot, Spps_plot,
labels = c("A", "B"))
# 保存到本地，由于我的电脑没有 `Arial` 字体，所以这里是报错的，故把 `family = "Arial"` 去掉
ggsave("gene_exp_plot.pdf", gene_exp_plot, width = 10, height = 5)
批量多个基因组图
一般我们需要挑大概十个左右基因来验证转录组数据的结果，这时候就可以这样做。
如果我们要绘制多个基因呢？就留点作业大家自己去参考 [Y 叔的王八拳教程](同一数据多变量分组的 boxplot?) 吧
给点提示：不能再多了
gene_list <- a[, 1][1:10]
test <- lapply(gene_list, my_barplot)
# 每一行 最多五个图
ten_gene <- plot_grid(plotlist = test, ncol = 5)
# 哈哈，当然字体就得自己调整一下了哈。
ggsave("ten_gene_exp_plot.pdf", ten_gene, width = 20, height = 5* length(test) %/% 5)
此文是二稿，第一稿就是只画了个图，才有了下面这句话
最后的最后作图一定要为了生物学问题而作图，切勿忘本！
后记
如果RNA-seq分析流程是自己走的，就会有bam文件或者bw文件，直接载入IGV也可以查看：
这里有一个番外，会使用IGV载入bam文件或者bw文件，但是你不一定有作者出图这样漂亮，我们后期会同步举办AI和PS学习交流群，一起来攻克发表级图像处理！
关于ATAC-seq数据处理系列连载目录：
关于RAN-seq和ChIP-seq数据处理
看B站免费视频就好！
号外：生信技能树全国巡讲11月在福州和上海，点击了解报名哈：（福州、上海见！）全国巡讲第19-20站（生信入门课加量不加价）