DEseq2差异表达分析

转录组数据绕不过差异分析。
为什么选择这两个包呢？
DEseq2针对有生物学重复的样本。（一般情况下应该是都需要生物学重复的）
edgeR对于单个样本是比较好的。（但是细胞材料真的没办法，细胞是又贵又难养。下一篇介绍。）
一、DEseq2差异表达分析
1、安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
library(DESeq2)

2、准备数据
featureCounts定量后的数据，或者FPKM数据（下一遍讲如何获取FPKM）
定量后的数据的数据格式如下：

image.png

colData，其实就是表型数据。
格式如下

image.png
这里可以看到，control有两个重复，A有三个重复，B有三个重复。
如果想看不同重复之间的是否相似（PCA，热图都行），你可以不计算平均值。
重复之间如果很相似，那就求平均，如果有些特别离谱，就不要了。

3、计算差异基因

library(tidyverse)
library(DESeq2)
#载入定量后的数据
setwd("E:/4/")
mycounts <- read.table("07counts_matrix.txt", header=TRUE)

#如果重复之间需要就平均值，那就计算下均值。
mycounts$Control <- round((mycounts$control2.sorted.bam + mycounts$control3.sorted.bam)/2)
mycounts$A <- round((mycounts$A1.sorted.bam + mycounts$A2.sorted.bam + mycounts$A3.sorted.bam)/3)
mycounts$B <- round((mycounts$B1.sorted.bam + mycounts$B2.sorted.bam + mycounts$B3.sorted.bam)/3)
mycounts <- subset(mycounts, select=c(Geneid,Control,A,B))

#前面我输入的时候没有把ensemble名字修改成symble，所以这里就不需要。
#如果处理成名字之后的数据会存在重复，所以需要把重复的删掉。
rows <- rownames(unique(mycounts['Geneid']))
mycounts <- mycounts[rows,]

#这里有个Geneid，需要去除，先把第一列当作行名来处理
rownames(mycounts)<-mycounts[,1]
#把带Geneid的列删除
mycounts <- mycounts[,-1]
#head(mycounts)

#载入colData文件
colData <-read.csv("08counts_matrix_condition.csv")
#colData <- colData[1:5,]
#colData <- colData[c(1,2,6,7,8),]
#构建数据矩阵
dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
dds <- DESeq(dds)

#接下来，我们要查看treatmentt VS control的总体结果，并根据p-value进行重新排序。
#利用summary命令统计显示一共多少个genes上调和下调（FDR0.1）
res <- results(dds, contrast=c("condition", "control", "treatment"))##或者res= results(dds)
res <- res[order(res$pvalue),]
summary(res)
#所有结果先进行输出
write.csv(res,file="DESeq2_results.csv")

#获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1.5或者小于-1.5的差异表达基因。
#DS和PD都是合理的，所以选择1.5。
diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5) 
##或者diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1.5 | log2FoldChange < -1.5))
dim(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file= "DESeq2_diffExpression.csv")

####提取出你所需要的log2FoldChange和pvalue列。
diff_gene_deseq2 <- na.omit(diff_gene_deseq2)
###这里把需要的两行提取出来先，是log2FoldChange和padj
nrDEG <- diff_gene_deseq2[,c(2,6)]
View(nrDEG)
###然后换个名字即可。
colnames(nrDEG) <- c('log2FoldChange','pvalue')
View( nrDEG)
write.csv(nrDEG,file= "DESeq2_diffExpression_logFC.csv")

4、数据可视化
4.1、火山图
加载包

library(ggplot2)
library(ggrepel)
library(export)

上步骤得到了差异基因，赋值给一个新的参数。
之所以这么干，是因为我这边是两个分开写的脚本。

DEG <- res
dim(DEG)
#这里用的是DESeq2的DEG，删掉NA。
DEG <- na.omit(DEG)
dim(DEG)
# 使用基础函数plot绘图
plot(DEG$log2FoldChange,-log2(DEG$padj))
# 确定差异表达倍数,abs表示绝对值
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)))
# 取前两位小数
logFC_cutoff <- round(logFC_cutoff, 2)
logFC_cutoff
查看一下logFC_cutoff的值。

# 确定上下调表达基因。
#方法1：按照logFC_cutoff绘图
#方法2：按照差异基因的log2FoldChange大于1.5或者其他值绘图。
DEG$change = as.factor(ifelse(DEG$padj < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
                        ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'STABLE'))
                        
DEG$change = as.factor(ifelse(DEG$padj < 0.05 & abs(DEG$log2FoldChange) > 1.5,
                        ifelse(DEG$log2FoldChange > 1.5 ,'UP','DOWN'),'STABLE'))
 
##确定需要在火山图上展示的字。
#this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,2),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
                    
#this_tile <- paste0('Cutoff for logFC is ',round(2.0,1),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))                    

###绘图的时候必须是dataframe，所以需要转换一下。
DEG <- data.frame(DEG)
## 颜色与分组一一对应
g <- ggplot(data=DEG, aes(x=log2FoldChange, y=-log10(padj),color=change)) + 
  geom_point(shape = 16, size=2) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2 fold change") + 
  ylab("-log10 p-value") +
  ggtitle( this_tile ) + ##这个可以不要
  theme(plot.title = element_text(size=15,hjust = 0.5)) + 
  theme_classic()+
  scale_colour_manual(values = c('blue','black','red'))
graph2ppt(file="Volcano_Plot.ppt", width=10, aspectr=1)

4.2、热图
载入包

library(pheatmap)

载入数据

dat <- mycounts
###载入表型数据
View(colData$condition)
# 提取差异倍数
# 根据需要修改DEG的值
FC <- DEG$log2FoldChange
###根据FC筛选。
#FC <- subset(FC, abs(FC)>1)
#View(FC)
names(FC) <- rownames(DEG)
#View(FC)

# 排序差异倍数，提取前100和后100的基因名
#或者全部的基因名字
DEG_200 <- c(names(head(sort(FC),100)),names(tail(sort(FC),100)))
DEG_all <- c(names(sort(FC)))
###将筛选得到的基因，提取并进行标准化。
dat <- t(scale(t(dat[DEG_200,])))
dat <- t(scale(t(dat[DEG_all,])))
###针对不同数据集，求平均值
dat.fram <- data.frame(dat)
dat.fram$Control <- (dat.fram$control2.sorted.bam + dat.fram$control3.sorted.bam)/2
dat.fram$A <- (dat.fram$A1.sorted.bam + dat.fram$A2.sorted.bam + dat.fram$A3.sorted.bam)/3
dat.fram$B <- (dat.fram$B1.sorted.bam + dat.fram$B2.sorted.bam + dat.fram$B3.sorted.bam)/3
dat2 <- subset(dat.fram, select=c(Control,A,B))
#这里需要重新导入表型数据的名字。
colData <-read.csv("09counts_matrix_condition_mean.csv")
colData <- colData[1:2,]
colData <- colData[c(1,3),]
group <- data.frame(colData$condition)
rownames(group)=colnames(dat2)
#绘图
#show_colnames =T，show_rownames = F，cluster_cols = T这三个参数根据需要设置
#颜色根据需要设置colorRampPalette(c("green", "balck", "red"))
g2 <- pheatmap(dat2,show_colnames =T,show_rownames = F, cluster_cols = T,
         annotation_col=group,border=FALSE,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50))     
graph2ppt(file="pheatmap_cluster_mean.ppt", width=10, aspectr=1)