差异表达分析——DEseq2

2022-04-14 本文已影响0人 Wei_Sun

差异基因筛选的常见方法有edgeR, DEseq2, EBseq等，这篇文章主要介绍DEseq2。

官方网站：
https://bioconductor.org/packages/release/bioc/html/DESeq2.html

官方说明书：
https://bioconductor.org/packages/release/bioc/manuals/DESeq2/man/DESeq2.pdf

安装环境：

R（4.1）
Bioconductor version: Release (3.14)

1.安装

DEseq2不在CRAN上，如果没有BiocManager的话需要提前安装。报错缺什么包再装什么包就可以，安装应该不困难。

> if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
> BiocManager::install("DESeq2")
> library(DESeq2)

2.数据准备

需要两个table：

counts.csv：转录组的基因表达值矩阵，第一列为基因名，随后每一列的列名为样本名。

count.csv

> mycounts <- read.csv("counts.csv")
> rownames(mycounts)<-mycounts[,1]
> mycounts<-mycounts[,-1]

colData.csv：样本的分组信息。第一列为样本名，第二列为分组信息。这里要注意，样本名要和counts.csv中的额样本名顺序一致。

colData.csv

> colData<-read.csv("colData.csv", row.names = 1)
#检查count文件和colData文件中的样本顺序是否一致
> all(rownames(colData) == colnames(mycounts))
[1] TRUE
#指定样本分组信息
> condition <- factor(colData$condition)
> condition = relevel( condition, "ref")

3.构建DEseqDataSet对象

这里需要注意，count数据中不能有缺失，如有缺失值，需要替换为0。

> any(is.na(mycounts))
[1] TRUE
> mycounts[is.na(mycounts)] <- 0

> dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)

4.计算差异倍数

> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

查看基因的差异倍数，和显著性p值。alt在前，ref 在后，意为 alt 相较于 ref中哪些基因上调/下调：

> res <- results(dds, contrast = c('condition', 'alt', 'ref'))

对p值重新进行排序：

> res = res[order(res$pvalue),]

res中，包含了基因id、标准化后的基因表达值、log2转化后的差异倍数（Fold Change）值、显著性p值以及校正后p值（默认FDR校正）等主要信息。

> head(res)
log2 fold change (MLE): condition alt vs ref 
Wald test p-value: condition alt vs ref 
DataFrame with 6 rows and 6 columns
                     baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                    <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
Parent=MC4G091700.1   5477.03        2.34803 0.0596426   39.3684  0.00000e+00  0.00000e+00
Parent=MP4G079000.1   2590.57       -2.93073 0.0772923  -37.9175  0.00000e+00  0.00000e+00
Parent=MP4G208400.1   9088.23        2.60702 0.0655940   39.7448  0.00000e+00  0.00000e+00
Parent=MP2G355500.1   3542.03       -2.51194 0.0712561  -35.2523 3.16856e-272 2.46830e-268
Parent=MC7G199500.1   5775.65       -1.85031 0.0585838  -31.5840 6.12845e-219 3.81925e-215
Parent=MC3G324200.1  16542.50       -1.74946 0.0588895  -29.7075 6.14876e-194 3.19325e-190

输出差异分析结果：

> write.csv(res,file="DEseq2_allresults.csv",quote = FALSE)

5.筛选差异表达基因

差异表达基因的界定很不统一，但log2FC是用的最广泛同时也是最不精确的方式，FC=2相对比较可靠；
采取log2FC和padj挑选差异基因；
获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。

标记上调和下调以及无差异基因

显著上调：log2FC ≥ 1 & padj < 0.05
显著下调：log2FC ≤ -1 & padj < 0.05
无差异基因：其余基因都为无差异基因

> res[which(res$log2FoldChange >= 1 & res$padj < 0.05),'sig'] <- 'up'
> res[which(res$log2FoldChange <= -1 & res$padj < 0.05),'sig'] <- 'down'
> res [which(abs(res$log2FoldChange) <= 1 | res$padj >= 0.05),'sig'] <- 'none'

输出差异基因总表

> diff_gene_deseq2 <- subset(res, sig %in% c('up', 'down'))
> write.csv(diff_gene_deseq2,file="diff_gene_deseq2.csv")

分别输出上调和下调基因

> res_up <- subset(res, sig == 'up')
> res_down <- subset(res, sig == 'down')
> write.csv(res_up, file = 'diff_gene_up.csv')
> write.csv(res_down, file = 'diff_gene_down.csv')

6.绘制火山图

> library(ggplot2)
> library(ggrepel)
> dat<-as.data.frame(res)
> pdf("volcano_plot.pdf",height=12,width=11)
> ggplot(dat,aes(x=log2FoldChange,y=-log10(padj),color=sig))+
  geom_point()+
  scale_color_manual(values=c("#CC0000","#BBBBBB","#2f5688"))+  #确定点的颜色
  theme_bw()+  #修改图片背景
  theme(
    legend.title = element_blank()  #不显示图例标题
  )+
  theme(axis.title.x =element_text(size=14,face = "bold"), axis.title.y=element_text(size=14,face = "bold"),axis.text = element_text(size = 14,face = "bold")) +  #调整坐标轴字号
  ylab('-log10 (p-adj)')+  #修改y轴名称
  xlab('log2 (FoldChange)')+  #修改x轴名称
  geom_vline(xintercept=c(-1,1),lty=3,col="black",lwd=0.5) +  #添加垂直阈值|FoldChange|>2
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.5)  #添加水平阈值padj<0.05
> dev.off()

火山图

引用转载请注明出处，如有错误敬请指出。