都是FPKM进行差异分析，为啥差异感觉这么大呢？

生信菜鸟团

发布于 2023-01-05 21:04:23

6.5K1

发布于 2023-01-05 21:04:23

文章被收录于专栏：生信菜鸟团

缘起

上周我们比较了FPKM与count分别进行差异分析的区别，发现两者差异分析的结果基本一致。恰巧曾老师叫我好好看看一篇发表在frontiers in Molecular Neuroscience的文章，他的转录组数据处理使用的就是FPKM矩阵。那这周的推文的话，我们就来看看我们的流程分析结果与作者有何不同。当然，强调一下哈，有count还是尽量用count为好。

探究

今天，我们使用2020年发表在frontiers in Molecular Neuroscience杂志上，文献名称为Differential Gene Expression Patterns Between Apical and Basal Inner Hair Cells Revealed by RNA-Seq的转录组FPKM数据（附表4）进行差异分析，感兴趣的小伙伴也可以自身下载一下下面的上游数据定量，同上一篇推文一样，比较两种数据类型进行差异分析的异同。

转录组数据集介绍

该数据集提交在ENA官网，项目号是PRJNA593359，如有小伙伴想看上游，链接见https://www.ebi.ac.uk/ena/browser/view/PRJNA593359 。该数据集组别上分为两组，一组是Basal组，一组是Apical组，每组均有三个样本。

正式分析

1.利用fpkm值进行差异分析

进行差异分析时，阈值按照文章的选法，FC为2，p<0.001

# 1.下载文中的补充数据集table 4.xlsx
rm(list=ls())
library(openxlsx)
library(stringr)
a=read.xlsx('table 4.xlsx',colNames = T) # 提取表达矩阵
class(a)
## [1] "data.frame"
colnames(a)=a[1,]
a=a[-1,]
table(duplicated(a$`Feature ID`))
## 
## FALSE  TRUE 
## 15307     1
a=a[!duplicated(a$`Feature ID`),]
rownames(a)=a$`Feature ID`
a=a[,-1]
keep <- rowSums(a>0) >= floor(0.5*ncol(a))
expMatrix <- a[keep,] #获得filter_count矩阵
#(防止0值影响FPKM转换为TPM值)

# 2.将FPKM值转换为TPM值
fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

class(expMatrix)
## [1] "data.frame"
expMatrix=na.omit(expMatrix)
numer=function(x){as.numeric(x)}
expMatrix=apply(expMatrix, 1, numer) 
#必须将矩阵中的字符类型由向量型改为字符型
expMatrix=t(expMatrix)
tpms <- apply(expMatrix,2,fpkmToTpm) 
tpms[1:3,]
##                    [,1]       [,2]     [,3]      [,4]     [,5]       [,6]
## 0610007P14Rik  47.67969 145.695029 22.27847   0.00000 117.7912  80.804891
## 0610009B22Rik 156.09341 300.652882 96.91044 104.00990 252.0514 224.045029
## 0610009L18Rik  30.88343   8.785541 39.00099  30.98678  10.0835   5.584516
colSums(tpms)
## [1] 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
colnames(tpms)=c("Apical1","Apical2","Apical3","Basal1","Basal2","Basal3")

# 3.差异分析
group_list=c(rep('Apical',3),rep('Basal',3))
## 强制限定顺序
group_list <- factor(group_list,levels = c("Apical","Basal"),ordered = F)
#表达矩阵数据校正
exprSet <- tpms
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

library(limma) 
head(exprSet)
##                 Apical1    Apical2  Apical3       Basal1    Basal2     Basal3
## 0610007P14Rik  47.67969 145.695029 22.27847   0.00000000 117.79116  80.804891
## 0610009B22Rik 156.09341 300.652882 96.91044 104.00989902 252.05143 224.045029
## 0610009L18Rik  30.88343   8.785541 39.00099  30.98678285  10.08350   5.584516
## 0610009O20Rik  70.01261  83.931085 46.90684   0.00000000  35.58989  11.294728
## 0610010B08Rik  10.48073  14.035653 27.61583   0.06863075  10.89523   1.544271
## 0610010F05Rik  13.90093  11.348727 14.68230   0.00000000  90.48093   0.000000
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

#判断数据是否需要转换
exprSet <- log2(exprSet+1)
range(exprSet)
## [1]  0.00000 14.31486
#差异分析：
dat <- exprSet
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
##          logFC AveExpr       t   P.Value adj.P.Val     B
## Adam8    5.364   2.682  13.180 4.176e-06   0.03458 4.131
## Ocm     -5.603   8.401 -12.813 5.033e-06   0.03458 4.018
## Nwd1     4.124   2.062   9.869 2.756e-05   0.06021 2.856
## Gm4869   3.901   1.951   9.800 2.882e-05   0.06021 2.823
## Rubcnl   3.858   1.929   9.769 2.942e-05   0.06021 2.807
## Igsf10   3.885   1.943   9.293 4.048e-05   0.06021 2.563
## Igfals   3.617   1.811   9.198 4.324e-05   0.06021 2.512
## Gm13547  3.582   1.791   9.123 4.556e-05   0.06021 2.471
## Bbox1    4.626   2.313   9.080 4.693e-05   0.06021 2.448
## Mmp23   -4.935   4.383  -9.050 4.794e-05   0.06021 2.431
bp=function(g){
  library(ggpubr)
  df=data.frame(gene=g,stage=group_list)
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg) 
##         logFC AveExpr       t   P.Value adj.P.Val     B
## Adam8   5.364   2.682  13.180 4.176e-06   0.03458 4.131
## Ocm    -5.603   8.401 -12.813 5.033e-06   0.03458 4.018
## Nwd1    4.124   2.062   9.869 2.756e-05   0.06021 2.856
## Gm4869  3.901   1.951   9.800 2.882e-05   0.06021 2.823
## Rubcnl  3.858   1.929   9.769 2.942e-05   0.06021 2.807
## Igsf10  3.885   1.943   9.293 4.048e-05   0.06021 2.563
#save(deg,file = 'deg.Rdata')

# 4.FPKM转换为TPM后，再利用limma进行差异分析

## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
if(T){
  logFC_t=1
  deg$g=ifelse(deg$P.Value>0.01,'stable',
               ifelse( deg$logFC > logFC_t,'UP',
                       ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
  )
  table(deg$g)
  head(deg)
  
}
##         logFC AveExpr       t   P.Value adj.P.Val     B    g
## Adam8   5.364   2.682  13.180 4.176e-06   0.03458 4.131   UP
## Ocm    -5.603   8.401 -12.813 5.033e-06   0.03458 4.018 DOWN
## Nwd1    4.124   2.062   9.869 2.756e-05   0.06021 2.856   UP
## Gm4869  3.901   1.951   9.800 2.882e-05   0.06021 2.823   UP
## Rubcnl  3.858   1.929   9.769 2.942e-05   0.06021 2.807   UP
## Igsf10  3.885   1.943   9.293 4.048e-05   0.06021 2.563   UP
fpkm_deg=deg

2.绘制文章的两张图

#  01复现第一张图，fpkm>1的Veen图，可以看看组别独有的基因
# 绘制思路：分别获得针对两个组别fpkm>1的基因，取个交集即可

apical=exprSet[,1:3]
apical=as.data.frame(apical)
apical$mean=apply(apical,1,mean)

basal=exprSet[,4:6]
basal=as.data.frame(basal)
basal$mean=apply(basal,1,mean)

# 韦恩图取个交集
Apical <- rownames(apical[apical$mean>1,])
Basal <- rownames(basal[basal$mean>1,])
library(ggvenn)
deg<-list(`Apical`=Apical,
          `Basal`=Basal)
p1 <- ggvenn(deg,
             show_percentage = F,
             stroke_color = "white",
             fill_color = c("#b2d4ec","#d3c0e2"),
             set_name_color = c("#ff0000","#4a9b83"))
p1 #30个FDR分析的差异基因全在利用P值分析的差异基因里


## 02差异基因箱线图
# 绘制箱线图
table(fpkm_deg$g=="UP")
## 
## FALSE  TRUE 
## 13454   287
table(fpkm_deg$g=="DOWN")
## 
## FALSE  TRUE 
## 13560   181
data=data.frame(c(287,181),row.names = c("up","down"))
colnames(data)=c("deg_number")
p2 <- ggplot(data,aes(x=rownames(data),y=deg_number,fill=rownames(data)))+
  geom_boxplot(width=0.6,alpha=0.8)
p3 <- p2+scale_y_continuous(breaks = seq(0,800,200),
                            limits=c(0,800))+theme_classic()
p3

3.用作者提出来的上调与下调基因在我的结果中验证

图就不修了哈，火山图的绘制很简单前面已经画了很多遍了，由于作者没有提供差异基因列表，我们就简单地针对其在文中提到的几个基因进行验证吧

cg=c("Prkd1","Sncg","Nme2","Fbxo32")
deg_fc=fpkm_deg[cg,]
head(deg_fc)
##          logFC AveExpr      t  P.Value adj.P.Val      B      g
## Prkd1  -1.1511   5.959 -2.853 0.025267    0.4069 -3.381 stable
## Sncg   -2.0205   5.932 -4.412 0.003308    0.1865 -1.377   DOWN
## Nme2    0.8174   7.355  1.333 0.225381    0.7532 -5.430 stable
## Fbxo32  1.3724   6.863  2.908 0.023378    0.4002 -3.305 stable

# 此处重点看logFC,再与文章

以上我们对该篇文章的FPKM数据集进行了处理，还绘制了文章的两幅结果图以及运用了文中的几个基因进行验证。韦恩图中可能是由于我设置了过滤参数导致我们的结果有一定区别。但是与原文相比，我们的上调差异基因有287个，下调差异基因有181个。与文章的644个上调，45个下调在数量上存在显著不同，这点就很值得咀嚼。

检查的几个差异基因源于文中的描述，如下。与3中我们的差异分析结果相比较，logFC大致是对应了，但是p值两者存在很大不同。为啥同一个数据，一样的阈值两者差异基因的数量相差如此显著，难道是统计方法的不同导致的差异基因在作者的分析结果与我的分析结果中存在不同吗？感兴趣的小伙伴们可以尝试尝试，欢迎点评哈。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-12-21，如有侵权请联系 cloudcommunity@tencent.com 删除

腾讯云开发者社区

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

腾讯云开发者社区

登录后参与评论

0 条评论

热度