1.加载所需数据

1.1 加载表达矩阵

gene_exp <- read.table("D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/genes.TMM.EXPR.matrix",
                       header = T, row.names = 1)

1.2 加载基因信息表

library(readr)
gene_info <- read_delim("D:/share/R_data/data/rnaseq-apple/query_seqs.fa.emapper.annotations", 
                                    delim = "\t", escape_double = FALSE, 
                                    col_names = FALSE, comment = "#", trim_ws = TRUE)

2.Tidyverse处理差异分析结果

2.1 导入差异分析结果

de <- read.table(file = 'D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/genes.counts.matrix.KID_S1_vs_BLO_S1.DESeq2.DE_results',
                 header = T)

这里不能使用rownanme = 1,因为tidyverse不能识别有行名的数据

2.2 选择保留列

colnames(de)
deg <- select(de, id, log2FoldChange, pvalue, padj) %>% 
  filter(abs(log2FoldChange) > 1 & padj < 0.05) %>% 
  mutate(FC = 2 ** log2FoldChange, direction = if_else(log2FoldChange > 1, 'up', 'down'))

使用colnames()函数加载差异分析结果中的列名，使用tidyverse包中的select()函数挑选需要保存的列，按照log2FoldChange > 1和padj < 0.05这两个条件进行筛选，使用mutate()函数在原表格基础上添加列，在这里添加了两列，一列是将log2FoldChange还原为常数，一列是对上下调基因做一个统计（通过if_else()条件进行判断）

2.3 合并表格

colnames(gene_info)
gene_info_se <- select(gene_info, Gene_Id, Gene_Symbol, Gene_Name)

对gene_info这个表格进行处理，使用select()函数挑选目标行并保存到gene_info_se这个向量中

de_result <- left_join(deg, gene_info_se, by = c('id' = 'Gene_Id'))

使用left_join()函数将deg和gene_info_se两个表组合，关联列是deg中的'id'和gene_info_se中的'Gene_Id'

3.代码简化

gene_exp <- read.table("D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/genes.TMM.EXPR.matrix", header = T, row.names = 1)
#加载表达矩阵
library(readr)
gene_info <- read_delim("D:/share/R_data/data/rnaseq-apple/query_seqs.fa.emapper.annotations", 
                                    delim = "\t", escape_double = FALSE, 
                                    col_names = FALSE, comment = "#", trim_ws = TRUE)
#加载基因信息表
de_result <- read.table(file = 'D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/genes.counts.matrix.KID_S1_vs_BLO_S1.DESeq2.DE_results',
                 header = T)
#加载差异分析结果
library(tidyverse)
de_result <- mutate(de_result, direction = if_else(padj > 0.05, 'ns', 
            if_else(abs(log2FoldChange) < 1, 'ns', if_else(log2FoldChange >= 1, 'up', 'down')))) %>%
#添加上下调信息
left_join(gene_info, by = c('id' = 'Gene_Id')) %>%
#关联基因信息
left_join(rownames_to_column(gene_exp, var = 'id'), by = 'id') %>%
#关联表达量信息
dplyr::select(-c(2:6, 8:9)) %>%
#去除无用的列
arrange(desc(abs(log2FoldChange)))
#按log2FoldChange绝对值降序排列
分组统计
group_by(de_result, direction) %>%
summarise(count = n())

保存与加载数据

save(gene_exp, gene_info, de_result, file = 'path/to/your/fileloader/rsave.rdata')
load(file = 'path/to/your/fileloader/rsave.rdata')

4.绘制火山图

library(EnhancedVolcano)
#加载EnhancedVolcano包
p1 <- EnhancedVolcano(de_result, 
                lab = de_result$id,
                x = 'log2FoldChange', 
                y = 'padj',
                legendPosition = "bottom",
                FCcutoff = 1, 
                pCutoff = 0.01,
                subtitle = NULL,
                title = NULL)

lab表示点的标签（以de_result向量中的id列为标签），x表示X轴上的数据（log2FoldChange），y表示Y轴上的数据（padj）,legendPosition表示图列所处的位置，FCcutoff，pCutoff表示对限制值进行限制， subtitle表示副标题为NULL， title表示标题为NULL

5.拼图

install.packages("cowplot")
library(cowplot)
plot_grid(p1, p2, labels = c('A', 'B', 'C', 'D', 'E'))

如果有好几张图片想把他们拼接到一张图上可以使用cowplot包实现，，p1,p2是要拼接的图片，labels表示图片标识