自己简单计算的变化倍数能跟引用近10万的算法对比？

自己简单计算的变化倍数能跟引用近10万的算法对比？

一个小伙伴表示他在处理两分组的转录组测序后的count矩阵的时候，发现手动计算的变化倍数跟金标准算法（DESeq2，edgeR，limma-voom）计算的不一样！

这个超级正常的，我们首先让人工智能大模型解释一下手动计算两分组的count矩阵中基因表达量的变化倍数（fold change）的过程：

变化倍数（fold change）

我自己写过一个代码：

load(file = './symbol_matrix.Rdata')
table(group_list)
group_list
symbol_matrix[1:4,1:4] ## 基因名字的样品，矩阵
dat = log2(edgeR::cpm(symbol_matrix)+1)
dat[1:4,1:4]
# 转换为长格式
library(reshape2)
melted_data <- melt(dat, 
                    varnames = c("Gene", "Sample"), 
                     = "ExpressionLevel")
library(ggplot2)
ggplot(melted_data, aes(x = ExpressionLevel, fill = Sample)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot with Genes of Interest", 
       x = "Expression Level", y = "Density") +
  theme_minimal()

mylogFC=rowMeans(dat[,4:6])-rowMeans(dat[,1:])

然后我们很容易去拿到金标准算法（DESeq2，edgeR，limma-voom）的差异分析结果进行对比：

load( file =  'deg-raw/DEG_deseq2.Rdata' )
head(DEG_deseq2)
df=data.frame(
  deseq2=DEG_deseq2$log2FoldChange,
  my=mylogFC[match(rownames(DEG_deseq2),names(mylogFC))]
)
head(df)
plot(df)
boxplot(dat['MKK1',] ~group_list)
DEG_deseq2['MKK1',]
df['MKK1',]
plot(abs(DEG_deseq2$log2FoldChange),log(DEG_deseq2$baseMean+1))

就可以看到，两者其实很容易冲突，因为金标准算法（DESeq2，edgeR，limma-voom）考虑到的事情更多，毕竟是引用近10万的算法啊！

引用近10万的算法

如果加上另外的两个包，轻轻松松破十万啊 ：

• 作者：MD Robinson · 2010 · 被引用次数：915
• 作者：ME Ritchie · 2015 · 被引用次数：1689

如果大家简简单单的随随便便的就能计算变化倍数， 那要这些引用近10万的算法干啥子？

其实大部分基因并不会冲突

人类的gtf文件里面记录的五六万基因里面，只有不到两万是蛋白质编码基因，去除5000低表达量的，剩下的一万多基因其实使用金标准算法（DESeq2，edgeR，limma-voom）拿到的变化倍数跟自己简单单的随随便便的就能计算变化倍数，并不会冲突很多。

但是有一些时候确实是会出现大面积冲突，因为可能会某个组里面的某一些样品有离点基因。正常情况下每个转录组测序的样品会两千万个reads，平均到两万个基因应该是每个基因1000左右的reads，但是等于核糖体和线粒体它们每个基因可以有几百万个reads，非常的浪费你的测序。这个时候，如果有一个样品的两千万个reads里面的一千万都是核糖体和线粒体，它就影响了你手动计算变化倍数，但是并不会影响金标准算法（DESeq2，edgeR，limma-voom）。因为金标准算法（DESeq2，edgeR，limma-voom）会考虑到这样的离点基因。