wilcox.test 进行 RNAseq 差异表达基因分析

如果拿到的表达谱不是原始的read counts数据而是TPM值，就不能用包（比如DESeq2）来进行差异表达基因分析。我们可以手动用wilcox.test 函数手动进行分析。我的数据为log2TPM的表达矩阵

ExM # log2TPM 表达矩阵
s1 # 属于类型1（如 tumor）的所有样本ID
s2 # 属于类型2（如 normal）的所有样本ID
cat("wilcox.test\n")
pvalue = padj = log2FoldChange = matrix(0, nrow(ExM), 1)
for(i in 1:nrow(ExM)){

    pvalue[i, 1] = p.value = wilcox.test(ExM[i, s1], ExM[i, s2])$p.value
    log2FoldChange[i, 1] = mean(ExM[i, s1]) - mean(ExM[i, s2])
}

padj = p.adjust(as.vector(pvalue), "fdr", n = length(pvalue))
rTable = data.frame(log2FoldChange, pvalue, padj, row.names = rownames(ExM))
treatment_Log2TPM <- signif(apply(ExM[rownames(rTable), s1], 1, mean), 4)
control_Log2TPM <- signif(apply(ExM[rownames(rTable), s2], 1, mean), 4)

cat("mark DGE\n") 
DGE <- rep("NC", nrow(ExM))
DGE[((rTable$padj) < 0.05) & (rTable$log2FoldChange > 0)] = "UP"
DGE[((rTable$padj) < 0.05) & (rTable$log2FoldChange < 0)] = "DN"
gene = rownames(ExM)
rTable = data.frame(treatment_Log2TPM, control_Log2TPM, rTable[, c("log2FoldChange", "pvalue", "padj")], DGE)
head(rTable)
                         treatment_Log2TPM control_Log2TPM log2FoldChange
ENSG00000166535.20 A2ML1             1.4870          1.8410  -0.3536611345
ENSG00000175899.15 A2M              14.3500         14.1000   0.2527242657
ENSG00000197953.6 AADACL2            0.1622          0.1491   0.0131439534
ENSG00000204518.2 AADACL4            0.1487          0.1492  -0.0005166819
ENSG00000115977.19 AAK1              9.6250          9.7070  -0.0817361189
ENSG00000127837.9 AAMP              11.2000         11.2000   0.0019474297
                              pvalue      padj DGE
ENSG00000166535.20 A2ML1  0.19797430 0.3930997  NC
ENSG00000175899.15 A2M    0.13120671 0.2997906  NC
ENSG00000197953.6 AADACL2 0.09516201 0.2405597  NC
ENSG00000204518.2 AADACL4 0.52208746 0.7067366  NC
ENSG00000115977.19 AAK1   0.44455824 0.6436331  NC
ENSG00000127837.9 AAMP    0.91096435 0.9563070  NC

这样就可以进行后续分析了。