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
