自己简单计算的变化倍数能跟引用近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"),
value.name = "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:3])
然后我们很容易去拿到金标准算法（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['MKNK1',] ~group_list)
DEG_deseq2['MKNK1',]
df['MKNK1',]
plot(abs(DEG_deseq2$log2FoldChange),log(DEG_deseq2$baseMean+1))
就可以看到，两者其实很容易冲突，因为金标准算法（DESeq2，edgeR，limma-voom）考虑到的事情更多，毕竟是引用近10万的算法啊！
引用近10万的算法
如果加上另外的两个包，轻轻松松破十万啊 ：
- 作者：MD Robinson · 2010 · 被引用次数：39135
- 作者：ME Ritchie · 2015 · 被引用次数：31689
如果大家简简单单的随随便便的就能计算变化倍数， 那要这些引用近10万的算法干啥子？
其实大部分基因并不会冲突
人类的gtf文件里面记录的五六万基因里面，只有不到两万是蛋白质编码基因，去除5000低表达量的，剩下的一万多基因其实使用金标准算法（DESeq2，edgeR，limma-voom）拿到的变化倍数跟自己简单单的随随便便的就能计算变化倍数，并不会冲突很多。
但是有一些时候确实是会出现大面积冲突，因为可能会某个组里面的某一些样品有离群点基因。正常情况下每个转录组测序的样品会两千万个reads，平均到两万个基因应该是每个基因1000左右的reads，但是等于核糖体和线粒体它们每个基因可以有几百万个reads，非常的浪费你的测序。这个时候，如果有一个样品的两千万个reads里面的一千万都是核糖体和线粒体，它就影响了你手动计算变化倍数，但是并不会影响金标准算法（DESeq2，edgeR，limma-voom）。因为金标准算法（DESeq2，edgeR，limma-voom）会考虑到这样的离群点基因。
腾讯云开发者
