【Rna-seq 分析流程】01.差异基因分析

1. 简介

差异分析是两组数据集间的比较，可以分为对照组样本(Control)，和目标组样本(可以是处理样本 - Treat，或疾病样本 - Illness)。
这里研究的是患有抗中性粒细胞胞浆抗体相关性血管炎（AAV ）与正常人之间存在的基因差异，这些差异基因可能是造成疾病的原因。

2. 数据信息

要使用上步处理后的数据：

1. exp_train.csv 筛选后的基因表达量
2. train_group.csv 样本分组信息

3. 思路

1. 利用 R 包 “limma” 获取差异基因，设定显著差异的条件是：|log2FC| > 1，adj.P < 0.05。
2. 根据 log2FC 值的大小判定基因是上调还是下调，并标记在刚刚生成的差异基因表格中。
3. 绘制火山图将上述结果可视化
4. 绘制热图将上述结果可视化

4. 代码

##---- 1.获取差异基因 -----
suppressMessages(library(limma))
suppressMessages(library(dplyr))
suppressMessages(library(ggplot2))
library(pheatmap)
library(tidyverse)
library(ggplotify)


###---- group_list --------
p2 = identical(colnames(exp_train),train_group$Sample);p2
train_group_list = train_group$label

###---- 差异基因 ----
design = model.matrix(~train_group_list)
fit = lmFit(exp_train,design)
fit = eBayes(fit)
DEG_o = topTable(fit,coef=2,number = Inf) 


##---- 2. 为deg数据框添加tags ----
DEG <- mutate(DEG_o,Symbol=rownames(DEG_o))
logFC_t=1
adj.P.Val_t = 0.05
k1 = (DEG$adj.P.Val < adj.P.Val_t)&(DEG$logFC < -logFC_t)
k2 = (DEG$adj.P.Val < adj.P.Val_t)&(DEG$logFC > logFC_t)
change = ifelse(k1,
                "down",
                ifelse(k2,"up","stable"))
DEG <- mutate(DEG,change)


###----- 3.火山图 -----
dat  = DEG
for_label_up <- dat%>% 
  filter(logFC > logFC_t & adj.P.Val < adj.P.Val_t) %>% 
  head(10)
for_label_down <- dat%>% 
  filter(logFC < -logFC_t  & adj.P.Val < adj.P.Val_t) %>% 
  head(10)
for_label = rbind(for_label_up,for_label_down)

p1 <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(adj.P.Val))) +
  geom_point(alpha=0.4, size=3, 
             aes(color=change)) +
  ylab("-log10(adj.Pvalue)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(adj.P.Val_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p1
p1 + geom_point( size=3, shape = 1, data = for_label) + 
  ggrepel::geom_label_repel(aes(label = Symbol), data = for_label, color="black", max.overlaps = Inf) 

ggsave("volcano.png", units="in", dpi=500, width=8, height=5, device="png")
ggsave("volcano.pdf", width = 8, height = 6)

###---- 4.热图 ----
x = DEG$logFC 
names(x)= rownames(DEG)
cg=c(names(head(sort(x),10)),names(tail(sort(x),10)))
n=exp_train[cg,]

annotation_col= column_to_rownames(train_group,"Sample") 
p3 = identical(colnames(n),rownames(annotation_col));p3  
colnames(annotation_col)[1] = "Group"
annotation_col[annotation_col == "AAV"] = 'AAV'
annotation_col[annotation_col == "C"] = 'Control'

heatmap_plot <- as.ggplot(pheatmap(n,
                                   name = "expression", 
                                   show_colnames = F,
                                   show_rownames = F,
                                   scale = "none", 
                                   annotation_col=annotation_col))
heatmap_plot
ggsave(heatmap_plot,filename = "heatmap.png", dpi=500, width=8, height=5, device="png")
ggsave(heatmap_plot,filename = "heatmap.pdf", width = 8, height = 6)

##---- 5. save ----
write.csv(DEG, file = "DEG.csv")

5. 结果展示

volcano.png

heatmap.png