手动计算logFC（wilcoxon差异分析）

医学和生信笔记

发布于 2023-09-12 18:50:59

7740

发布于 2023-09-12 18:50:59

文章被收录于专栏：医学和生信笔记医学和生信笔记

前段时间有小伙伴问怎么手动计算logFC，今天说一下。

logFC是log fold change的缩写，也就是log之后的差异倍数。这个差异倍数意思是某个基因在A组表达量的平均值是B组表达量平均值的几倍。

这个东西的计算其实很简单的，就是常规的对数计算而已。

一般来说，我们用tpm或者fpkm时，通常都会先进行log2处理，在log2处理后的表达矩阵里，如果某个基因在A组表达量是x，在B组表达量是y，那么这个基因的logFC = x - y。

这不是巧合，只是一个很简单的数学公式log(x/y)=log(x)-log(y)

准备数据

用gse87466这个GEO的数据做演示，下载整理的过程这次就不演示了。数据在粉丝qq群文件，需要的加群下载即可。

library(easyTCGA)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)

load(file = "G:/easyTCGA_test/gse87466.Rdata")

这是一个炎症性肠病的数据集，一共108个样本，21个normal，87个uc（ulcerative colitis）。

探针注释我已经提前做好了，但是有一些探针对应多个symbol，为了方便我这里直接删掉了：

exprSet[1:4,1:4]
##                                           GSM2332098 GSM2332099 GSM2332100
## IGK@ /// IGKC                               13.86197   13.76880   13.95740
##                                             13.95740   13.92619   13.79664
## IGL@                                        13.73797   13.61266   13.86197
## IGH@ /// IGHA1 /// IGHA2 /// LOC100126583   13.79664   13.16844   13.76880
##                                           GSM2332101
## IGK@ /// IGKC                               13.95740
##                                             13.86197
## IGL@                                        13.76880
## IGH@ /// IGHA1 /// IGHA2 /// LOC100126583   13.73797

exprSet <- exprSet[!grepl("/",rownames(exprSet)),]

group <- factor(group_list,levels = c("normal","UC"))
table(group)
## group
## normal     UC 
##     21     87

使用limma做差异分析

首先对这个数据做下差异分析，也是用easyTCGA包，1行代码即可，基因芯片数据也是支持的。

如果是count矩阵，会自动使用DESeq2、limma、edgeR进行差异分析；
如果是tpm、fpkm、基因表达芯片数据，它会自动检测需不需要进行log2转换，然后进行wilcoxon和limma的差异分析：

library(easyTCGA)

diff_res <- diff_analysis(exprset = exprSet
                          , group = group
                          , is_count = F
                          )
## => log2 transform not needed
## => Running limma
## => Running wilcoxon test
## => Analysis done.

# limma的结果
diff_limma <- diff_res$deg_limma

diff_limma <- diff_limma %>% 
  arrange(desc(logFC))

head(diff_limma)
##            logFC   AveExpr         t      P.Value    adj.P.Val        B
## MMP3    5.125542  9.319122 12.754123 2.387453e-23 1.132335e-20 42.56916
## SLC6A14 5.025462  9.414419 21.582496 3.792392e-41 7.554445e-37 82.65877
## DUOX2   4.912367  9.917941 17.377496 3.456356e-33 1.721265e-29 64.81419
## DEFB4A  4.609705  8.710986  8.923018 1.217313e-14 5.214811e-13 22.81705
## CHI3L1  4.554672  9.777758 12.195694 4.281634e-22 1.332658e-19 39.72614
## S100A8  4.259343 10.016995 12.112299 6.602449e-22 1.906098e-19 39.29950
##         genesymbol
## MMP3          MMP3
## SLC6A14    SLC6A14
## DUOX2        DUOX2
## DEFB4A      DEFB4A
## CHI3L1      CHI3L1
## S100A8      S100A8

现在有很多文章中直接使用的wilcoxon检验，但是它并不能计算logFC：

diff_wilc <- diff_res$deg_wilcoxon
head(diff_wilc)
##              pvalue genesymbol
##        9.202462e-03           
## IGL@   1.182949e-06       IGL@
## RPL41  1.708300e-02      RPL41
## RPL24  2.428482e-01      RPL24
## RPL37A 1.903898e-01     RPL37A
## RPS11  8.255636e-03      RPS11

自己计算logFC

根据前面的理论，我们可以自己计算logFC，思路就是分别计算某个基因在两组中的平均表达量，然后直接相减即可。

下面我们用dplyr中的rowwise操作实现这一过程，当然还有其他方法，选择自己喜欢的即可。

library(dplyr)
library(tidyr)

uc <- colnames(exprSet)[group_list == "UC"]
normal <- colnames(exprSet)[group_list == "normal"]

logfc_df <- exprSet %>% 
  rowwise() %>% 
  mutate(mean_uc=rowMeans(across(all_of(uc))),
         mean_normal=rowMeans(across(all_of(normal))),
         logfc = mean_uc - mean_normal, # 这个就是logFC了
         .keep = "none"
            ) %>% 
  bind_cols(genesymbol = rownames(exprSet)) %>% 
  arrange(desc(logfc))
  
head(logfc_df)
## # A tibble: 6 × 4
## # Rowwise: 
##   mean_uc mean_normal logfc genesymbol
##     <dbl>       <dbl> <dbl> <chr>     
## 1   10.3         5.20  5.11 MMP3      
## 2   10.4         5.37  5.01 SLC6A14   
## 3   10.9         5.95  4.92 DUOX2     
## 4    9.60        5.00  4.60 DEFB4A    
## 5   10.7         6.10  4.56 CHI3L1    
## 6   10.8         6.55  4.29 S100A8

可以看到，我们手动计算的这个logfc和上面limma包计算的logFC基本上是一样的（有误差，可以忽略）哦。