gzh原文点我，欢迎同时关注

---

本系列为 RNA-Seq 分析具体实现过程，将具体描述流程及可执行代码～
本次复现文章来自 Nature communication，
题目为：Temporal and Spatial Heterogeneity of Host Response to SARS-CoV-2 Pulmonary Infection

本篇将根据经过Deseq2分析过的数据和基因，以火山图和热图的形式进行展示
文章共提供了88个sra文件，共计61533个基因，在进行差异分析之前，先去除了表达量小于50的基因，剩余24126个，最终得到5785个高表达基因和261个低表达基因，在今天的内容中，我们将以火山图和热图的方式，展示关键基因

Volcano 火山图

---

1、输入数据

#构建dds矩阵,标准化,以及进行差异分析
dds <- DESeqDataSetFromMatrix(mycounts_data, coldata, design= ~ condition)
dds$condition <- relevel(dds$condition, ref = "high")
dds <- DESeq(dds)
res= results(dds, alpha = 0.01)
res_data<-data.frame(res)

#基因名转换
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
mms_symbols_res<- getBM(attributes=c('ensembl_gene_id','external_gene_name','description'),
                    filters = 'ensembl_gene_id', values =row.names(res_data), mart = mart)
res_data<-rownames(res_data) %>% cbind(res_data)
colnames(res_data)[1]<-c("ensembl_gene_id")
res_data<-merge(res_data,mms_symbols_res,by="ensembl_gene_id")

2、对数据进行处理

#选中关键基因
selected_genes <- c('IFIF3','IFI6','IFI44L','OAS3','ADAR','SFTPB','MUC2','KRT4','SFTPB',
                    'OAS2','IFIH1','STAT1','MLPH','RELN','ACOXL','PTPRN','IFIT3','RSAD2','ISG15','RSAD2','SFTPB')

#去除padj为na的行，并对处于不同位置的基因进行分类设置
my_de_result <- res_data %>%
  filter(!is.na(padj)) %>%
  mutate(padj = replace_na(padj, 1)) %>%
  mutate(direction = if_else(padj < 0.01 & abs(log2FoldChange) >= 1,
                             'star',
                             if_else(abs(log2FoldChange) >= 1,
                                     'green',
                                     if_else(padj >= 0.01 ,'blue','grey'
                                     ))))%>%
  mutate(selected = if_else(external_gene_name %in% selected_genes,
                            'Yes', direction))

3、画图：

ggplot(data = my_de_result, aes(x = log2FoldChange, 
                                y = -log10(padj))) +
  geom_point(size = 2,
             aes(color = direction),
             show.legend = F,
             alpha=0.6) +
  geom_point(data = filter(my_de_result, 
                           selected == 'Yes'), #filter()过滤数据；为选择的点加圆框
             shape = 21, 
             size = 2, 
             color = 'black', 
             stroke = 0.8) + # 边框粗细
  geom_text_repel(data = filter(my_de_result, 
                                selected == 'Yes'),
                  size = 5, box.padding = 0.5,
                  aes(label = external_gene_name)) + #名称，使用ggrepel避免重合
  geom_hline(yintercept = -log10(0.01), 
             linetype = 'dotdash', color = 'grey30') + #画虚横线
  geom_vline(xintercept = c(-1, 1), 
             linetype = 'dotdash', color = 'grey30') + #画两条垂直竖线
  scale_color_manual(
    breaks = c("blue", "star","green", "grey"),
    values = c('#4a67d9','#db3b22','#689f5d','#4d4d4d')) +
  ylim(0, 15) +
  labs(x = 'Log2(fold change)', 
       y = '-log10(p-value)') +
  theme_bw()

4、火山图展示：

Volcano

ComplexHeatmap

---

在 pheatmap 的基础上， 我们使用更为高阶的 ComplexHeatmap 自定义地展示数据：

1、数据处理部分：

#从res 中选出 491个最显著的基因 padj < 0.0001 & abs(log2FoldChange) > 1，基因名存储在diff_name 中
#使用经过Deseq2归一化后的数据进行绘制
DeseqNorm_data<-counts(dds,normalized=TRUE)
diff<-data.frame(diff_name)
heatmap_data<-data.frame(DeseqNorm_data)
colnames(heatmap_data)<-coldata$address
heatmap_data<- heatmap_data %>%
  mutate(ensembl_gene_id = row.names(heatmap_data))
heatmap_data<-merge(heatmap_data,diff,by="ensembl_gene_id") %>%
  dplyr::select(-c(48:53)) %>%
  filter(external_gene_name!="") %>%
  dplyr::select(-c(49,1)) %>%
  column_to_rownames(var="external_gene_name")

2、自定义部分：

#定义颜色
col_fun = colorRamp2(c(-0.5, 1, 2.5), c('#2878B5','white','#c82423'),
                     space = "RGB")
#定义注释：文字(样本名称注释)和类别注释
ha = HeatmapAnnotation(bar=coldata$condition,
                       col = list(bar=c("high"='#c82423',"low"='#2878B5')),
                       show_annotation_name =FALSE,
                       show_legend = FALSE,
                       which='column',
                       text = anno_text(colnames(heatmap_data),
                                        just = "left", offset = unit(1, "mm"),
                                        gp = gpar(fontsize = 9)))
#自定义类别图例
lgd3 = Legend(labels = c("Virus high", "Virus low"),
              legend_gp = gpar(fill=c('#2878B5','#c82423')),
              direction="vertical",
              row_gap = unit(3, "mm"),
              grid_width= unit(8, "mm"）)

3、设置主图样式并绘制：

#主图绘制
heatmap=Heatmap(log10(heatmap_data+1),
        col = col_fun,#配色
        top_annotation = ha,
        column_names_side = "top",
        column_dend_side = "top",
        show_row_names = F,
        show_column_names = F,
        row_dend_width = unit(2, "cm"),#改变树宽
        column_dend_height = unit(2, "cm"),#改变树高，
        heatmap_legend_param = list( #热图图例自定义
          at = c(-2,0,2),
          title="log10(DeseqNorm+1)",
          legend_direction="horizontal",
          legend_width = unit(4, "cm"),
          title_position="topcenter",
          title_gp = gpar(fontsize = 10,fontface = "bold")))
#绘制并定义热图图例位置
draw(heatmap, heatmap_legend_side = "bottom")
#绘制类别图例
draw(lgd3,just = "top",x = unit(1.25, "cm"), y = unit(18, "cm"))

4、最终结果：

ComplexHeatmap

往期推荐

5 类图形为你展示—变量间相互关系—提供示例与完整代码
来势汹汹的Omicron-BA.2谱系
简易代码！一网打尽ggplot2箱线图，小提琴图，蜂窝图和云雨图
读好文｜零门槛了解RNA-seq