差异基因筛选的常见方法有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()
火山图
引用转载请注明出处，如有错误敬请指出。
网友评论