復(fù)習(xí)資料：

表觀調(diào)控13張圖之一證明基因干擾有效性

表觀調(diào)控13張圖之二相關(guān)性熱圖看不同樣本相關(guān)性

表觀調(diào)控13張圖之三卵慰。。榴芳。

在10月份的時(shí)候跟著學(xué)習(xí)了表觀調(diào)控相關(guān)的聯(lián)合組學(xué)分析嗡靡，當(dāng)時(shí)做了點(diǎn)記錄，但是很多東西還是不熟悉窟感，這次跟著公眾號(hào)開(kāi)始復(fù)習(xí)叽躯。

此次表觀調(diào)控課程是整合RNA-seq數(shù)據(jù)和表觀數(shù)據(jù)，比如Chip-seq和ATAC-seq肌括，主要是依托于文章 Global changes of H3K27me3 domains and Polycomb group protein distribution in the absence of recruiters Spps or Pho 文獻(xiàn)解讀在 在果蠅探索 PRC 復(fù)合物（逆向收費(fèi)讀文獻(xiàn) 2019-18） 這一推文点骑。
原文鏈接：https://www.pnas.org/content/115/8/E1839

先回顧一下背景知識(shí)：

• H3K27me3 是發(fā)生在組蛋白 H3 的第 27 位賴氨酸上的三甲基化修飾，從整體染色質(zhì)水平來(lái)講, 大量的證據(jù)支持 H3K27me3 與異染色質(zhì)相關(guān)谍夭；從個(gè)體的基因水平來(lái)講, H3K27me3 常常與轉(zhuǎn)錄抑制相關(guān)黑滴。H3K27me2/me3 甲基化的加載是由 PcG 復(fù)合體催化的。

• 如何精確地募集 Polycomb group（PcG）蛋白到其靶基因仍知之甚少紧索。在果蠅中袁辈，PcG蛋白被募集到由多個(gè) DNA 結(jié)合蛋白的結(jié)合位點(diǎn)組成的 Polycomb 反應(yīng)元件（ PRE ）。為了了解如何將 PcG 蛋白募集到 PRE 并在其上維持珠漂，作者系統(tǒng)地研究了野生型的 PcG 結(jié)合晚缩，相關(guān)的 H3K27me3 和轉(zhuǎn)錄組以及三種 PRE 結(jié)合蛋白的突變體（ Pho、Spps媳危、cg ）荞彼。

• 僅 Pho 結(jié)合位點(diǎn)不足以招募 PcG。PRE 由復(fù)雜的 DNA 結(jié)合蛋白的結(jié)合位點(diǎn)組成待笑，包括Pho / Phol鸣皂，Spps，Cg暮蹂，GAF / Psq寞缝，Adf-1，Grh仰泻，Dsp1和Zeste / Fsh-S荆陆。然而這些蛋白質(zhì)在 PcG 蛋白質(zhì)募集中的作用尚不清楚。Spps 是與 PREs 結(jié)合的 Sp1 / KLF 鋅指蛋白家族的成員

第一張圖：說(shuō)明基因干擾的有效性

當(dāng) PcG 招募成員被中斷后集侯，基因表達(dá)是如何改變的被啼。進(jìn)一步探索差異表達(dá)與 H3K27me3 修飾的改變有何關(guān)系帜消？
作者對(duì)野生型、Spps 和 Pho 突變體 3 齡幼蟲(chóng)進(jìn)行了 RNA-seq 測(cè)序趟据。而本張圖用來(lái)證實(shí)相應(yīng)突變體中 Spps 和 Pho 表達(dá)確實(shí)降低了券犁。

注： 當(dāng)在檢測(cè)時(shí)發(fā)現(xiàn)基因突變體的 RNA-seq 結(jié)果發(fā)生改變的時(shí)候术健，首先在獲得表達(dá)量時(shí)候汹碱，第一件事情是去查看這些數(shù)據(jù)提示的基因在突變體中是否真的被抑制或者不表達(dá)了。

加載本次分析所涉及的包

rm(list = ls())
library(ggpubr)
library(stringr)
library(ggplot2)
library(cowplot)

自定義 barplot 主題

定義過(guò)的主題以后還可以重復(fù)使用

my_bar_theme <- theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 16), 
                      # 因?yàn)?x 軸標(biāo)簽要旋轉(zhuǎn) 90°荞估，所以這里用來(lái)旋轉(zhuǎn)
                      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)) # 基因名要居中咳促，這里用來(lái)居中。

數(shù)據(jù)讀取

options(stringsAsFactors = F)
a <- read.table('~/Downloads/fly/RNA-seq/counts/all.counts.id.txt', header = T)
# 查看數(shù)據(jù)行列信息勘伺，大致看下多少基因
dim(a)

> dim(a)
[1] 17714    16

因?yàn)槲覀兊臄?shù)據(jù)來(lái)源是一樣的跪腹，所有結(jié)果也應(yīng)該一樣，由于存放的位置不同飞醉，路徑做修改即可冲茸。

數(shù)據(jù)清洗、可視化

這一步非常重要缅帘，而且有幾個(gè)難點(diǎn)轴术，首先是獲取數(shù)據(jù)，然后是排序钦无，最后是取色逗栽，這里前輩用了Takecolor 軟件進(jìn)行對(duì)準(zhǔn)樣圖取色，我是不會(huì)的失暂，下去還需要學(xué)習(xí)一下彼宠，這里以及后面都可以可以參考曾老師的御用視頻剪輯師的方法來(lái)完成。
鏈接：生信 | IGV導(dǎo)出svg文件后期美化

# 提取 `pho` 基因?qū)?yīng)的行以及表達(dá)量信息
cg <- a[a[,1] == 'pho',7:16]

# 構(gòu)建新的數(shù)據(jù)框以便進(jìn)行可視化
dat <- data.frame(gene = as.numeric(cg),
                  sample = str_split(names(cg),'\\.',simplify = T)[,1],
                  group = str_split(names(cg),'_',simplify = T)[,1]
)

# 重定義 x 軸標(biāo)簽順序
labels = c(paste0("WT", "_",1:3), paste0("PhoKO", "_", 1:3), paste0("SppsKO", "_", 1:4))

# 然后使用 `ggbarplot()` 函數(shù)進(jìn)行繪圖弟塞，其實(shí)我個(gè)人更喜歡用 `ggplot2` 來(lái)繪制凭峡，這個(gè)包也是基于 `ggplot2` 來(lái)的，所以對(duì)于我來(lái)說(shuō)沒(méi)啥區(qū)別决记，反而還變麻煩了想罕。
ggbarplot(dat,x = 'sample', y = 'gene', 
          color = 'black', fill = "group", 
          size = 1) +
  # 使用 Takecolor 軟件進(jìn)行對(duì)準(zhǔn)樣圖取色
  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 = "")) 

得到結(jié)果

封裝畫(huà)圖函數(shù)

為了方便以后調(diào)用，因?yàn)橐?huà)多個(gè)基因的圖來(lái)展示（也是偷懶的一種方式霉涨，哈哈）

my_barplot <- function(gene){
  cg <- a[a[,1] == gene, 7:16] # 提取候選的表達(dá)量對(duì)應(yīng)的行和列

  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 軸標(biāo)簽的順序按价，這里是按照原圖的順序來(lái)的
  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 
}

組合兩個(gè)基因的圖

Pho_plot <- my_barplot("pho")
Spps_plot <- my_barplot("Spps")
gene_exp_plot <- plot_grid(Pho_plot, Spps_plot, 
                           labels = c("A", "B"))

# 保存到本地，由于我的電腦沒(méi)有 `Arial` 字體笙瑟，所以這里是報(bào)錯(cuò)的楼镐，故把 `family = "Arial"` 去掉
ggsave("gene_exp_plot.pdf", gene_exp_plot, width = 10, height = 5) 

結(jié)果2

批量多個(gè)基因組圖

gene_list <- a[, 1][1:5]
test <- lapply(gene_list, my_barplot)
# 每一行 最多五個(gè)圖
five_gene <- plot_grid(plotlist = test, ncol = 5)
# 當(dāng)然字體需自己調(diào)整一下了。

ggsave("five_gene_exp_plot.pdf", five_gene, width = 20, height = 5* length(test) %/% 5)

結(jié)果3

第二張圖：相關(guān)性熱圖看不同樣本相關(guān)性

當(dāng)我們拿到數(shù)據(jù)時(shí)候凄杯，除了前面的質(zhì)控等分析外，我們一般需要查看樣品內(nèi)的重復(fù)性怎么樣秉宿，一般目前市面上的 RNA-seq戒突、ChIP-seq 測(cè)序樣品內(nèi)的相關(guān)性都能高達(dá) 0.9 以上。

可以通過(guò)兩種策略來(lái)計(jì)算樣品內(nèi)的相關(guān)性
1描睦、根據(jù)基因的表達(dá)量信息來(lái)計(jì)算樣品之間的相關(guān)性膊存，比如 RNA-seq 。
2忱叭、將全基因組等分 bin 的方法隔崎，然后計(jì)算每個(gè) bin 里面的 reads 數(shù), 然后通過(guò)均一化等過(guò)程，再對(duì)數(shù)據(jù)進(jìn)行計(jì)算相關(guān)性, 比如 ChIP-seq 等 DNA 類(lèi)型測(cè)序數(shù)據(jù)韵丑。

熱圖一 通過(guò)基因的表達(dá)量來(lái)計(jì)算樣品相關(guān)性

rm(list = ls())
options(stringsAsFactors = F)
a = read.table('~/Downloads/fly/RNA-seq/counts/all.counts.id.txt', header = T)
dim(a)
dat = a[,7:16]
library(stringr)
ac = data.frame(group=str_split(colnames(dat),'_',simplify = T)[,1])
rownames(ac) = colnames(dat)
M=cor(log(dat+1))
pheatmap::pheatmap(M,
                   annotation_col = ac,
                   # breaks = seq(0, 100, length.out = 100)
) 

熱圖1

熱圖二 分析deeptools軟件的multiBigwigSummary和plotCorrelation得到的相關(guān)性結(jié)果

• linux/終端中運(yùn)行：

multiBigwigSummary  bins -b  *.bw -o results.npz -p 4
plotCorrelation -in results.npz \
  --corMethod spearman --skipZeros \
  --plotTitle "Spearman Correlation of Read Counts" \
  --whatToPlot heatmap --colorMap RdYlBu --plotNumbers \
  --plotFileFormat pdf \
  -o heatmap_SpearmanCorr_readCounts.pdf   \
  --outFileCorMatrix SpearmanCorr_readCounts.tab

• 將結(jié)果導(dǎo)入 R

b = read.table('SpearmanCorr_readCounts.tab', header = T)
colnames(b) <- str_split(colnames(b),'\\.',simplify = T)[, 1]
rownames(b) <- str_split(colnames(b),'\\.',simplify = T)[, 1]
bc = data.frame(group = str_split(colnames(b),'_',simplify = T)[, 1])
rownames(bc) <- colnames(b)
pheatmap::pheatmap(b,
                   annotation_col = bc) 

熱圖2

可以很清楚的看到，樣品內(nèi)的重復(fù)性都是極高的撵彻。同一樣品都聚類(lèi)在一起钓株。樣品內(nèi)的相關(guān)性顯著高于樣品間的相關(guān)性。說(shuō)明數(shù)據(jù)重復(fù)性很好陌僵，可以進(jìn)行進(jìn)下一步轴合。（與前輩的有不一樣，他展示了分析后的拾弃，我的就是原來(lái)是數(shù)據(jù)值桩，沒(méi)有進(jìn)一分析。主要是關(guān)于這個(gè)）

第三張圖 DESeq2尋找差異基因并可視化

DESeq2尋找差異基因

if(T){
  # 賦值表達(dá)矩陣和分組信息為一個(gè)新的變量豪椿，方便以后直接調(diào)用
  exprSet = dat
  suppressMessages(library(DESeq2)) # 加載包

  (colData <- data.frame(row.names = colnames(exprSet), 
                         group_list = group_list) )

  # 構(gòu)建一個(gè)表達(dá)矩陣
  dds <- DESeqDataSetFromMatrix(countData = exprSet,
                                colData = colData,
                                design = ~ group_list)
  dds <- DESeq(dds)

# 下面我們得到 Spps 突變后的差異基因    
  res <- results(dds, 
                 contrast=c("group_list", "SppsKO", "WT")) 
# 注意這里是前面比后面奔坟，別把順序搞反了，到時(shí)候一不注意結(jié)果就是反的搭盾。所以拿差異倍數(shù)對(duì)著原始 reads 比較一下咳秉。

  resOrdered <- res[order(res$padj),] # 按照 padj 值進(jìn)行排序
  head(resOrdered)
  DEG_SppsKO = as.data.frame(resOrdered) # 將數(shù)據(jù)轉(zhuǎn)變?yōu)?data.frame 數(shù)據(jù)框
  DEG_SppsKO = na.omit(DEG_SppsKO) # 去除包含 NA 值的行

# 下面我們得到 Pho 突變后的差異基因 
  table(group_list)
  res <- results(dds, 
                 contrast = c("group_list","PhoKO","WT"))
  resOrdered <- res[order(res$padj),]
  head(resOrdered)
  DEG_PhoKO = as.data.frame(resOrdered)
  DEG_PhoKO = na.omit(DEG_PhoKO)

# 將關(guān)鍵結(jié)果以 Rdata 形式保存到本地，下次如有需要就不需要重新用 DESeq2 計(jì)算了    
  save(DEG_PhoKO, DEG_SppsKO,file = 'deg_output.Rdata')
}

可視化

MAplot: 圖 a 和 圖 b
繪制 Spps

# apeglm 方法需要安裝 apeglm 包
BiocManager::install("apeglm")

# ashr 方法同樣需要安裝額外依賴的包
BiocManager::install("ashr")

Spps_resLFC <- lfcShrink(dds, coef  = "group_list_SppsKO_vs_PhoKO", type = "apeglm")
Spps_resNorm <- lfcShrink(dds, coef = "group_list_SppsKO_vs_PhoKO", type = "normal")
Spps_resAsh <- lfcShrink(dds, coef  = "group_list_SppsKO_vs_PhoKO", type = "ashr")

par(mfrow=c(1,3), mar=c(4,4,2,1))
xlim <- c(1,1e5); ylim <- c(-3,3)

plotMA(Spps_resLFC, xlim  = xlim, ylim = ylim, main = "apeglm")
plotMA(Spps_resNorm, xlim = xlim, ylim = ylim, main = "normal")
plotMA(Spps_resAsh, xlim  = xlim, ylim = ylim, main = "ashr")

dev.off()

plotMA繪制3種算法下的Spps結(jié)果

MAplot 在DESeq2 包中表示的是變化倍數(shù) log2(Fold change) 與所有樣本標(biāo)準(zhǔn)化后的 counts 數(shù)的平均值之間的關(guān)系鸯隅。
怎么得到這個(gè)結(jié)果呢：
lfcShrink 收縮 FC 三種方法如下：

• normal is the the original DESeq2 shrinkage estimator, an adaptive Normal distribution as prior. This is currently the default, although the default will likely change to apeglm in the October 2018 release given apeglm’s superior performance.

• apeglm is the adaptive t prior shrinkage estimator from the apeglm package (Zhu, Ibrahim, and Love 2018).

• ashr is the adaptive shrinkage estimator from the ashr package (Stephens 2016). Here DESeq2 uses the ashr option to fit a mixture of Normal distributions to form the prior, with method="shrinkage".

這三種方法澜建，后兩種都是可以隨便選擇的。
同樣的方法繪制另一個(gè)基因的圖

Pho_resLFC <- lfcShrink(dds, coef  = "group_list_WT_vs_PhoKO", type = "apeglm")
Pho_resNorm <- lfcShrink(dds, coef = "group_list_WT_vs_PhoKO", type = "normal")
Pho_resAsh <- lfcShrink(dds, coef  = "group_list_WT_vs_PhoKO", type = "ashr")

par(mfrow=c(1,3), mar=c(4,4,2,1))
xlim <- c(1,1e5); ylim <- c(-3,3)

plotMA(Pho_resLFC, xlim  = xlim, ylim = ylim, main = "apeglm")
plotMA(Pho_resNorm, xlim = xlim, ylim = ylim, main = "normal")
plotMA(Pho_resAsh, xlim  = xlim, ylim = ylim, main = "ashr")
dev.off()

plotMA繪制3種算法下的Pho結(jié)果

par(mfrow=c(1,2), mar=c(4,4,2,1))
plotMA(Spps_resAsh, xlim = xlim, ylim = ylim, main = "Spps_Ash")
plotMA(Pho_resAsh, xlim = xlim, ylim = ylim, main = "Pho_Ash")
dev.off()

同時(shí)用第三種方法展示2個(gè)基因的結(jié)果

差異基因韋恩圖:UpSetR/VennDiagram

第一種：我們使用 R 包 UpSetR 來(lái)繪制差異基因之間的韋恩圖（ 多組時(shí)候蝌以，這種更加一目了然 ）

library(UpSetR)
load(file = 'deg_output.Rdata')

colnames(DEG_PhoKO)
DEG_PhoKO$change = ifelse(DEG_PhoKO$padj > 0.05,'stable',
                          ifelse(DEG_PhoKO$log2FoldChange > 0, 'up', 'down'))

DEG_SppsKO$change=ifelse(DEG_SppsKO$padj>0.05, 'stable',
                         ifelse(DEG_SppsKO$log2FoldChange > 0, 'up', 'down'))

SppsKO_up   = rownames(DEG_SppsKO[DEG_SppsKO$change == 'up',])
SppsKO_down = rownames(DEG_SppsKO[DEG_SppsKO$change == 'down',])
PhoKO_up    = rownames(DEG_PhoKO[DEG_PhoKO$change == 'up',])
PhoKO_down  = rownames(DEG_PhoKO[DEG_PhoKO$change == 'down',])

allG = unique(c(SppsKO_up, SppsKO_down, PhoKO_up, PhoKO_down))

df = data.frame(allG = allG,
              SppsKO_up   = as.numeric(allG %in% SppsKO_up),
              SppsKO_down = as.numeric(allG %in% SppsKO_down),
              PhoKO_up    = as.numeric(allG %in% PhoKO_up),
              PhoKO_down  = as.numeric(allG %in% PhoKO_down))

upset(df)

UpSetR 來(lái)繪制差異基因之間的韋恩圖

第二種：我們使用 VennDiagram來(lái)展示炕舵，也是就是文中那種花瓣圖

# 這里直接 copy 琪音同學(xué)的
# 鏈接: https://mp.weixin.qq.com/s/Pg0mjz7mD73atMnHz7jv1A
# 導(dǎo)入本地字體，不然 `Arial` 字體會(huì)警告
#BiocManager::install("extrafont")
library("extrafont")
font_import()
loadfonts()

load(file = 'deg_output.Rdata')

colnames(DEG_PhoKO)
DEG_PhoKO$change = ifelse(DEG_PhoKO$padj > 0.05,'stable',
                          ifelse(DEG_PhoKO$log2FoldChange > 0, 'up', 'down'))

DEG_SppsKO$change=ifelse(DEG_SppsKO$padj>0.05, 'stable',
                         ifelse(DEG_SppsKO$log2FoldChange > 0, 'up', 'down'))

SppsKO_up   = rownames(DEG_SppsKO[DEG_SppsKO$change == 'up',])
SppsKO_down = rownames(DEG_SppsKO[DEG_SppsKO$change == 'down',])
PhoKO_up    = rownames(DEG_PhoKO[DEG_PhoKO$change == 'up',])
PhoKO_down  = rownames(DEG_PhoKO[DEG_PhoKO$change == 'down',])

library(VennDiagram)
venn.diagram(
  x = list(
    "Up in phoKO"    = PhoKO_up,
    "Down in phoKO"  = PhoKO_down,
    "Up in SppsKO"   = SppsKO_up,
    "Down in SppsKO" = SppsKO_down
    ),
  filename       = "Venn.png", # 保存到當(dāng)前目錄
  # stroke 顏色 形式 粗細(xì)
  col            = "#000000", lwd = 3, lty = 1,
  # 填充 透明度
  # 顏色可以參考前一篇跟畅，使用 takecolor 自己取色
  fill           = c("#D3E7F0", "#9FBEDB", "#A0D191", "#6DAE8A"),
  alpha          = 0.50,
  # 里外字體大小形式
  cex            = 1.5,
  fontfamily     = "Arial",
  fontface       = "bold",
  cat.cex        = 1.5,
  cat.dist       = 0.2,
  cat.fontfamily = "Arial",
  # 圖像整體位置 分辨率
  margin         = 0.2,
  resolution     = 300)

花瓣韋恩圖就保存在當(dāng)前目錄下

與文章趨勢(shì)基本一致咽筋。可以看出 Spps 和 Pho 共同調(diào)控許多基因徊件，說(shuō)明這兩基因有一定的關(guān)系奸攻。

兩樣本 log2FC 相關(guān)性散點(diǎn)圖

load(file = 'deg_output.Rdata')
library(ggpubr)
DEG_PhoKO = DEG_PhoKO[rownames(DEG_SppsKO),]
po = data.frame(SppsKO = DEG_SppsKO$log2FoldChange,
              PhoKO = DEG_PhoKO$log2FoldChange)

sp <- ggscatter(po, 'SppsKO', 'PhoKO',
                add        = "reg.line",  # Add regressin line
                add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
                conf.int   = TRUE # Add confidence interval
)
# Add correlation coefficient
sp + stat_cor(method  = "pearson", 
              label.x = -10, label.y = 10) # 相關(guān)系數(shù)顯示位置

相關(guān)性系數(shù)還是可以的蒜危，很有顯著性

兩者的相關(guān)系數(shù)高達(dá)0.53，這在兩個(gè)基因間是算具有很強(qiáng)的相關(guān)的相關(guān)性了睹耐，更加佐證了上圖的韋恩圖的結(jié)果辐赞。

友情宣傳：

• 廣州專(zhuān)場(chǎng)（全年無(wú)休）GEO數(shù)據(jù)挖掘課，帶你飛（1.11-1.12）
• 生信入門(mén)課全國(guó)巡講2019收官--長(zhǎng)沙站