R语言 DESeq2 基因差异分析 简单备注版 火山图

RT

#加载数据 linba代表淋巴
lnc <- read.csv("K0.csv")
#保证ID的唯一性，删除重复
table(duplicated(lnc$id))
index<-duplicated(lnc$id)
us_count<-lnc[!index,]  #反向筛选
#转化成矩阵，rowname转化成第一列
A<-us_count
rownames(A)<-A[,1] 
A<-A[,-1] 

#Count with fewer than 10 reads were removed.
A[A < 10] <- 0
#将输入数据取整
A<-round(A,digits=0) 
A
#准备，将数据转换为矩阵格式
A<-as.matrix(A) 
## 设置分组信息，建立环境(个样本，组处理)
colnames(A)
dim(A)
condition <- factor(c(rep("control",179),rep("case",80)))
coldata<-data.frame(row.names=colnames(A),condition)  
coldata 
#展示coldata内容
condition  
#展示环境

# 然后判断现在表达矩阵与样本信息是否一致，特别是针对很多样本时。
all(rownames(coldata) == colnames(A))

library("DESeq2") 
#使用library函数加载DEseq2包
##构建dds矩阵 
dds<-DESeqDataSetFromMatrix(A,coldata,design=~condition)
head(dds) 
#去除表达量非常少的行，比如设置阈值为每行表达量为10。
#基因至少在某一些文库的count超过10 ~ 15 才被认为是表达。
keep <- rowSums(counts(dds)) >= 2000
dds <- dds[keep,]
dim(dds)  #看过滤了
#之前不设置时，默认control在case后面。 设置后，
#control组(untreated)排在了前面
#dds$condition <- factor(dds$condition, levels = c("control","case"))
dds$condition <- relevel(dds$condition, ref = "control")
levels(dds$condition)


##进行差异分析
#对原始的dds进行标准化
dds<-DESeq(dds) 
#查看结果名称
resultsNames(dds)  
#用results函数提取结果，并赋值给res变量
res<-results(dds) 
#查看结果
plotMA(res,ylim=c(-2,2)) 
mcols(res,use.names=TRUE)
#简单绘制火山图
plot(res$log2FoldChange,res$pvalue) 

#提取差异基因
res <- res[order(res$padj),]
resdata<-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),
               by="row.names",sort=FALSE)
deseq_res<-data.frame(resdata)

#提取上调差异表达基因
up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 2)) 
#提取下调差异表达基因
down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -2)) 
#输出上调基因
write.csv(up_diff_result,"linbaUP.csv") 
#输出下调基因
write.csv(down_diff_result,"linbaDOWN.csv") 
#输出总的差异
write.csv(deseq_res,"linbaUpDown.csv")


#简单验证下方向问题，看前面设置的对照组问题，
#和上面导出的数据癌和对照的方向有没有错
library(tidyverse)
train <- up_diff_result[,-c(1:7)]
train <- t(train)
dim(train)
train <- data.frame(train)
train$group <-'case'  
dim(train)
train$group[1:179] <- 'control'
train$group <- factor(train$group)
#通过平均值验证方向。
#aggregate(x, by, FUN)
aggregate(train[,1],list(train[,19]),mean)
#

#火山图
#提取差异的数据
resdata <- mutate(resdata,select_change=if_else(abs(log2FoldChange) >= 1,'y','n'),
                  select_pvalue=if_else( padj < 0.05,'y','n')  )

resdata <- mutate(resdata,select=paste(resdata[,'select_change'], 
                                       resdata[,'select_pvalue'], sep = ''))
table(resdata$select)

library("ggplot2")

##ggplot2 差异火山图
resdata$select <- factor(resdata$select, levels = c('nn', 'ny', 'yn', 'yy'), 
                         labels = c('p >= 0.05, FC < 2', 'p < 0.05, FC < 2', 
                                    'p >= 0.05, FC >= 2', 'p < 0.05, FC >= 2'))

#纵轴为显著性 p 值
volcano_plot_pvalue <- ggplot(resdata, aes(log2FoldChange, -log(padj, 10))) +
  geom_point(aes(color = select), alpha = 0.6) +
  scale_color_manual(values = c('gray30', 'green4', 'blue2','red2')) +
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', 
                                        fill = 'transparent')) +
  theme(legend.title = element_blank(), 
        legend.key = element_rect(fill = 'transparent'), 
        legend.background = element_rect(fill = 'transparent')) +
  geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.5) + 
  geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.5) +
  labs(x = 'log2 Fold Change', y = '-log10 p-value')
volcano_plot_pvalue
#保存为PDF
#ggsave('volcano_plot_pvalue.pdf', volcano_plot_pvalue, width = 6, height = 5)