文章主要讨论了在处理转录组测序的count矩阵时,手动计算基因表达量的变化倍数与金标准算法(如DESeq2,edgeR,limma-voom)之间的差异。文章强调了金标准算法在处理这种情况时的优势,特别是在处理含有离群点基因的样品时。
文章指出在处理转录组测序的count矩阵时,手动计算基因表达量的变化倍数与金标准算法存在差异,这是正常现象。
文章描述了一种手动计算变化倍数的流程,包括加载count矩阵、数据预处理、差异表达基因的识别等步骤。
文章强调了金标准算法(如DESeq2,edgeR,limma-voom)在处理转录组测序数据时的优势,特别是在处理含有离群点基因的样品时。
文章提到了在某些情况下手动计算和算法计算的变化倍数会出现大面积冲突,主要是由于样品中可能存在离群点基因的影响。
一个小伙伴表示他在处理两分组的转录组测序后的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 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万的算法
如果加上另外的两个包,轻轻松松破十万啊 :