search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
120 stars 30 forks source link

转录组差异分析P值与FDR值区别有多大 #2681

Closed ixxmu closed 2 years ago

ixxmu commented 2 years ago

https://mp.weixin.qq.com/s/N3EYU8TNANh7nSzxt0MnIw

ixxmu commented 2 years ago

转录组差异分析P值与FDR值区别有多大 by 生信菜鸟团

缘起

上一期转录组专题我们讲了 转录组数据三分组的差异分析方式，有小伙伴针对其中的内容提到为什么edgeR做差异分析用的是P值，而不用校正后的FDR值，不是说FDR经过多重校正，能够更好地降低假阳性吗？这个小伙伴说的非常对，但为了筛选更多的差异基因，以及避免丢失一些基因，我们通常也会考虑过滤条件不那么严格的P值，所以这是一个假阳与假阴的权衡。关于这个权衡，我们将用一个转录组数据集通过edgeR进行差异分析，去展示P值与FDR值差异分析的区别。

试验

今天，我们使用标题为 NAT10 depletion suppresses fatty acid metabolism in MCF7 breast cancer cells的数据集GSE210086进行分析，数据集的介绍页面如下：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE210086，感兴趣的小伙伴可以进入看看作者的总体设计。时间紧的话，可粗略了解下，作者主要通过敲低NAT10基因看其对MCF7乳腺癌细胞系脂肪酸代谢的影响，具体脂肪酸代谢与乳腺癌表型的描述可能需要小伙伴去查查代谢与肿瘤相关的文章了，因为这个数据集的文章目前还没有发表。

转录组数据集介绍

GSE210086数据集的样本分组如下，敲低对照组有三个重复，敲低NAT10组有三个重复

处理数据的话，庆幸的是，作者上传了基因count矩阵，我们就可以直接走上次count数据处理的差异分析流程。值得注意的点是，作者上传的是基因表达矩阵，这意味着我们无需进行ID转换，直接过滤进行差异分析即可

以下为基因count矩阵下载链接：https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE210086&format=file&file=GSE210086%5Fgene%5Fcount%2Exls%2Egz

接下来，我们重点关注P值与多重校正后的FDR值差异分析的区别

正式分析

0.清控环境，加载包

rm(list = ls())
library(edgeR)
library(DESeq2)
library(FactoMineR)
library(factoextra)
library(org.Hs.eg.db)
library(stringr)
library(stringi)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(VennDiagram)
library(RColorBrewer)

1.读取数据，获得基因表达矩阵以及cpm矩阵

rawcount <- read.delim("./GSE210086_gene_count.xls.gz",header = T,sep="\t")
rawcount <- rawcount[,2:8] #获得正式的基因表达矩阵
table(duplicated(rawcount$gene_name))#随机降重
## 
## FALSE  TRUE 
## 57167  1568
rawcount$median <- apply(rawcount[1:6],1,median)
rawcount <- rawcount[!duplicated(rawcount$gene_name),]
rownames(rawcount) <- rawcount$gene_name
rawcount <- rawcount[,1:6] #获得正式的基因表达矩阵
head(rawcount)
##        siNAT10_MCF2 siNAT10_MCF3 siC_MCF1 siNAT10_MCF1 siC_MCF2 siC_MCF3
## EEF1A1        82029        79914   133846        79414    96249    81294
## PKM          113137        87066   103953       142330    53897    45121
## ENO1         106702        81082   111201       119585    41414    34417
## FN1           54799        88855    33452        66222    86689   119331
## RPLP0         96851        63370    90415        93733    38245    33715
## IGFBP3        67465        79051    60625        65019    69043    61142
# 将组别换成正常组在前，对照组在后的习惯去处理额
rawcount=rawcount[,c(3,5,6,1,2,4)]
head(rawcount) #获取正常组在前，疾病组在后的基因表达矩阵
##        siC_MCF1 siC_MCF2 siC_MCF3 siNAT10_MCF2 siNAT10_MCF3 siNAT10_MCF1
## EEF1A1   133846    96249    81294        82029        79914        79414
## PKM      103953    53897    45121       113137        87066       142330
## ENO1     111201    41414    34417       106702        81082       119585
## FN1       33452    86689   119331        54799        88855        66222
## RPLP0     90415    38245    33715        96851        63370        93733
## IGFBP3    60625    69043    61142        67465        79051        65019
keep <- rowSums(rawcount>0) >= floor(0.5*ncol(rawcount))
filter_count <- rawcount[keep,] #获得filter_count矩阵
express_cpm <- log2(cpm(filter_count)+ 1)
express_cpm[1:4,1:4] #获得cpm矩阵
##        siC_MCF1 siC_MCF2 siC_MCF3 siNAT10_MCF2
## EEF1A1 13.46855 13.32252 13.18932     12.72864
## PKM    13.10394 12.48606 12.34009     13.19245
## ENO1   13.20117 12.10605 11.94950     13.10797
## FN1    11.46852 13.17161 13.74302     12.14676
save(express_cpm,filter_count,file="Step03filterGeneCountCpm.Rdata")

2.获取差异分析的分组信息

#根据列的信息，提取分组信息
colnames(rawcount) 
## [1] "siC_MCF1"     "siC_MCF2"     "siC_MCF3"     "siNAT10_MCF2" "siNAT10_MCF3"
## [6] "siNAT10_MCF1"
group=rep(c("shCON","siNAT10"),each=3)
group_list=factor(group,levels = c("shCON","siNAT10"))
table(group_list)#检查一下组别数量
## group_list
##   shCON siNAT10 
##       3       3

3.初看两组分析的数据与样本分组分布（箱线图与PCA图）

## 01绘制整体表达的箱线图
exprSet <- express_cpm
data <- data.frame(expression=c(exprSet),
                   sample=rep(colnames(exprSet),each=nrow(exprSet)))
p1 <- ggplot(data = data,aes(x=sample,y=expression,fill=sample))+ geom_boxplot() + theme(axis.text.x = element_text(angle = 90)) + 
  xlab(NULL) + ylab("log2(CPM+1)")+theme_bw()

## 02绘制PCA图
dat <- express_cpm
dat <- as.data.frame(t(dat))
dat <- na.omit(dat)
dat$group_list <- group_list
dat_pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#画图仅需数值型数据，去掉最后一列的分组信息
p2 <- fviz_pca_ind(dat_pca,
                   geom.ind = "point", # 只显示点，不显示文字
                   col.ind = dat$group_list, # 用不同颜色表示分组
                   palette = c("#00AFBB", "#E7B800"),
                   addEllipses = T, # 是否圈起来，少于4个样圈不起来
                   legend.title = "Groups") + theme_bw()
p1+p2

4.分别利用P值与FDR值进行差异分析

# 01先获取整体的差异分析表达结果
exprSet <- filter_count
design <- model.matrix(~0+rev(factor(group_list))) 
rownames(design) <- colnames(exprSet)
colnames(design) <- levels(factor(group_list))
DEG <- DGEList(counts=exprSet,  #构建edgeR的DGEList对象
               group=factor(group_list))
DEG$samples$lib.size 
## [1] 11808834  9396381  8704124 12087295 10388406 12481825
DEG <- calcNormFactors(DEG)
DEG$samples$norm.factors
## [1] 0.9260817 1.1034025 1.1501644 0.9244119 1.0255480 0.8975017
# 计算线性模型的参数
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)
# 拟合线性模型
fit <- glmFit(DEG, design)
lrt <- glmLRT(fit, contrast=c(1,-1)) 
# 提取过滤差异分析结果
DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG)))
head(DEG_edgeR)
##           logFC   logCPM       LR       PValue          FDR
## TAGLN  2.279278 7.298149 56.93375 4.507442e-14 6.611515e-10
## CRYAB  1.513615 5.833886 32.39699 1.256810e-08 9.217444e-05
## SAMD11 2.424716 2.454854 24.80752 6.334993e-07 3.097389e-03
## LMCD1  1.313811 6.322715 23.30588 1.381766e-06 5.066936e-03
## TSSC4  1.637568 4.057177 22.79470 1.802621e-06 5.288170e-03
## BOK    1.411509 6.456652 21.09287 4.375522e-06 1.069669e-02
# 设定阈值，筛选显著上下调差异基因
fc <- 2.0
p <- 0.05
## 02利用P值进行差异分析
DEG_edgeR$regulated <- ifelse(DEG_edgeR$logFC>log2(fc)&DEG_edgeR$PValue<p,
                              "up",ifelse(DEG_edgeR$logFC<(-log2(fc))&DEG_edgeR$PValue<p,"down","normal"))
table(DEG_edgeR$regulated)
## 
##   down normal     up 
##    252  14055    361
## 03利用FDR值进行差异分析
fc <- 2.0
p <- 0.05
DEG_edgeR$regulated1 <- ifelse(DEG_edgeR$logFC>log2(fc)&DEG_edgeR$FDR<p,
                               "up",ifelse(DEG_edgeR$logFC<(-log2(fc))&DEG_edgeR$FDR<p,"down","normal"))
table(DEG_edgeR$regulated1)
## 
##   down normal     up 
##      5  14638     25

5.展示P值与FDR值分析的差异分析结果（火山图与热图）

## 01绘制P值与FDR值差异分析的火山图
# 绘制第一组用P值进行差异分析的火山图
data1 <- DEG_edgeR
p3.1 <- ggplot(data=data1, aes(x=logFC, y=-log10(PValue),color=regulated)) + 
  geom_point(alpha=0.5, size=1.8) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2FC") + ylab("-log10(Pvalue)") +
  scale_colour_manual(values = c('blue','black','red'))+theme_bw()

# 绘制第二组用FDR值进行差异分析的火山图
p3.2 <- ggplot(data=data1, aes(x=logFC, y=-log10(FDR),color=regulated1)) + 
  geom_point(alpha=0.5, size=1.8) + 
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2FC") + ylab("-log10(FDRvalue)") +
  scale_colour_manual(values = c('blue','black','red'))+theme_bw()

# 拼图
p3.1+p3.2
## 02绘制P值与FDR值差异分析的热图
# 绘制P值分析后比较的差异基因热图
edgeR1_sigGene1 <- DEG_edgeR[DEG_edgeR$regulated!="normal",]
edgeR1_sigGene1 <-rownames(edgeR1_sigGene1)
dat1 <- express_cpm[match(edgeR1_sigGene1,rownames(express_cpm)),]
dat1[1:4,1:4]
##        siC_MCF1 siC_MCF2 siC_MCF3 siNAT10_MCF2
## TAGLN  6.011491 5.856408 5.584638     8.187205
## CRYAB  4.968841 5.055129 4.974774     6.439499
## SAMD11 1.678857 1.045610 1.504981     3.555214
## LMCD1  5.702199 5.617758 5.506579     6.830937
group1 <- data.frame(group=group_list)
rownames(group1)=colnames(dat1)
# 绘制热图
library(pheatmap)
p4.1 <- pheatmap(dat1,scale = "row",show_colnames =T,show_rownames = F, 
                 cluster_cols = T, 
                 annotation_col=group1,
                 main = "edgeR's DEG")
## 绘制FDR值分析后比较的差异基因热图
edgeR1_sigGene2 <- DEG_edgeR[DEG_edgeR$regulated1!="normal",]
edgeR1_sigGene2 <-rownames(edgeR1_sigGene2)
dat2 <- express_cpm[match(edgeR1_sigGene2,rownames(express_cpm)),]
dat2[1:4,1:4]
##        siC_MCF1 siC_MCF2 siC_MCF3 siNAT10_MCF2
## TAGLN  6.011491 5.856408 5.584638     8.187205
## CRYAB  4.968841 5.055129 4.974774     6.439499
## SAMD11 1.678857 1.045610 1.504981     3.555214
## LMCD1  5.702199 5.617758 5.506579     6.830937
group1 <- data.frame(group=group_list)
rownames(group1)=colnames(dat2)
# 绘制热图
library(pheatmap)
p4.2 <- pheatmap(dat2,scale = "row",show_colnames =T,show_rownames = F, 
                 cluster_cols = T, 
                 annotation_col=group1,
                 main = "edgeR's DEG")
# 热图拼图
library(ggplotify)
p4.1 <- as.ggplot(p4.1)
p4.2 <- as.ggplot(p4.2)
p4.1+p4.2

6.整体看一下P值与adjustP值分析差异基因的交集情况吧

degP <- rownames(DEG_edgeR[DEG_edgeR$regulated!="normal",])
degFDR <- rownames(DEG_edgeR[DEG_edgeR$regulated1!="normal",])
table(degP%in%degFDR)
## 
## FALSE  TRUE 
##   583    30
library(ggvenn)
deg<-list(`P`=degP,
          `FDR`=degFDR)
p5 <- ggvenn(deg,
             show_percentage = F,
             stroke_color = "white",
             fill_color = c("#b2d4ec","#d3c0e2"),
             set_name_color = c("#ff0000","#4a9b83"))
p5 #30个FDR分析的差异基因全在利用P值分析的差异基因里

以上就是P值与FDR值进行差异分析的区别，咱们可以看见当组别之间样本没有很明显的区分时（聚类效果不好），用P值与用FDR值进行差异分析具有巨大区别，差异基因数量显著减少。

接着，回到咱们关于P值与FDR值进行差异分析的假阳与假阴的权衡上来，如果咱们以FDR值来过滤，会得到很少的显著基因，这些基因可信度可能很高，但是有些显著基因可能因为这条件太严被过滤了。

然而，如果用P值进行差异分析的话，就可能会面临许多显著基因中存在假阳性，且可能会被某些杂志的reviewer argue或者拒收，所以如果数据分组效果可以，想得到较少的假阳基因，推荐用FDR值。

但是如果想进行筛选，数据质量不行，可以考虑用P值，当然这是小编的一家之言，见仁见智，大家也可以谈谈用FDR值以及P值进行差异分析的见解，以及试试用两种值进行差异分析，欢迎大家在留言区阐述自身的看法以及结果哈。