表观调控13张图之一证明基因干扰有效性

rm(list = ls()) # 日常运行代码前清理一下环境
library(ggpubr)
library(stringr)
library(ggplot2)
library(cowplot)

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)) # 基因名要居中，这里用来居中。

options(stringsAsFactors = F)
a <- read.table('all.counts.id.txt', header = T)
# 查看数据行列信息，大致看下多少基因
dim(a)
[1] 17714    16

# 提取 `pho` 基因对应的行以及表达量信息
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]
               )

# 重定义 x 轴标签顺序
labels = c(paste0("WT", "_",1:3), paste0("PhoKO", "_", 1:3), paste0("SppsKO", "_", 1:4))

# 然后使用 `ggbarplot()` 函数进行绘图，其实我个人更喜欢用 `ggplot2` 来绘制，这个包也是基于 `ggplot2` 来的，所以对于我来说没啥区别，反而还变麻烦了。
ggbarplot(dat,x = 'sample', y = 'gene', 
          color = 'black', fill = "group", 
          size = 1) +
  # 使用 Takecolor 软件进行对准样图取色
  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) 

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)