写在前面的话
分享一下利用DESeq2包寻找差异基因的小趣事，更让我觉得生物学意义和统计学意义结合的重要性......尤其处理大型非模式动物的数据时，要多长个心眼~~~
library(DESeq2)
library(dplyr)
library(tibble)
library(tidyr)
mycounts <- read.csv("7MDEGs.csv", header = T)
exprSet<-aggregate(x=mycounts[,2:(ncol(mycounts))],by=list(mycounts$gene.id),FUN=mean)
exprSet <- column_to_rownames(exprSet, var='Group.1')
exprSet<- exprSet[-1,]
mycounts =exprSet
keep <- rowSums(mycounts>0) >= floor(0.75*ncol(mycounts))
table(keep)
filter_count <- rawcount[keep,]
condition <- factor(c(rep("CM", 3), rep("IM", 3)), levels = c("CM","IM"))
colData <- data.frame(row.names = colnames(mycounts), condition)
dds <- DESeqDataSetFromMatrix(round(keep), colData, design = ~condition)
dds <- DESeq(dds)
res <- results(dds,contrast=c("condition", "IM", "CM"))
resOrdered <- res[order(res$padj),]
DEG =as.data.frame(resOrdered)
DEG
有意思的是会发现很多NA
官方文档的note
NOTE: on p-values set to NA
1.If within a row, all samples have zero counts, the baseMean column will be zero, and the log2 fold change estimates, p-value and adjusted p-value will all be set to NA.
2.If a row contains a sample with an extreme count outlier then the p-value and adjusted p-value will be set to NA. These outlier counts are detected by Cook’s distance.
3.If a row is filtered by automatic independent filtering, for having a low mean normalized count, then only the adjusted p-value will be set to NA.
直接用DESeq2结果进行后续分析
deg = na.omit(DEG)
fc_cutoff <- 2
p <- 0.05
deg$regulated <- "normal"
loc_up <- intersect(which(deg$log2FoldChange>log2(fc_cutoff)),
which(deg$pvalue<p))
loc_down <- intersect(which(deg$log2FoldChange< (-log2(fc_cutoff))),
which(deg$pvalue<p))
deg$regulated[loc_up] <- "up"
deg$regulated[loc_down] <- "down"
write.csv(deg,"deg.csv")
最近因分析单细胞数据关注到某个基因，并经实验验证，却一直以为转录组水平检测不到差异水平。细看了一下经DESeq2分析后的结果将该基因p值判断为NA。
查看该counts值
发现存在明显离群值，故都将p值和adj值判断为NA，符合官方文档说明的第二点。但生物学意义是人为赋予的，允许在该模型中存在这样离群现象，且具有生物学意义。可能有时候完全依赖数据工具分析结果，是远远不够的，也要反过来从源头看数据。
对该现象，所以为什么转录组筛选差异基因时候有多种方式，DESeq2包，limma包，edgeR包一起使用进行筛选。或者单独可视化自身感兴趣的基因集吧~~~~~
