如果拿到的表達譜不是原始的read counts數(shù)據(jù)而是TPM值，就不能用包（比如DESeq2）來進行差異表達基因分析。我們可以手動用wilcox.test 函數(shù)手動進行分析烫映。我的數(shù)據(jù)為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
這樣就可以進行后續(xù)分析了。
