我可以使用DESeq2的结果子集来计算Bray-Curtis相异性，然后绘制PCoA吗

Question

在比较疾病和对照之间细菌基因的差异表达，然后计算Bray-Curtis相异性并随后绘制PCoA时，我很难使用DESeq2的结果。

我已将DESeq2的输出保存为数据帧。它由6000行基因名称和两列组成，一列表示p值(均<0.05(，另一列表示log2FOldChange>1。数据帧称为siggenes1。在运行Bray Curtis和PCoA之前，我需要规范化我的数据吗？我以为这已经通过DESEq2完成了，但看看我可以提供的代码，在执行DESEq2时，我还没有包括归一化=t。

或者，在使用DESeq2之前，我是否需要使用扫描函数对初始数据进行归一化？

Bray-Curtis差异的我的代码

vegDistOut=vegdist(t(siggenes1),"bray")

以上得到1个值，即0.995。现在我有点不知道如何设计出用它绘制PCoA的代码，因为我的下一段代码是错误的。

pcoaOut=pcoa(vegDistOut)

数组错误(STATS，dims[perm](："dims"的长度不能为0

由于以上步骤，我无法继续。

如果有人能帮忙，我将不胜感激。谢谢

Answer 1

欢迎使用StackOverflow。

Bray-Curtis相似性通常用于确定两个样本的物种组成有多相似。典型的输入包括来自kraken、clark等软件的每个样本的物种计数，用于高通量数据，或者在16S的情况下，Qiime2：或dada2。

|     Genus     |  Sample 1   |  Sample 2  |
|---------------|-------------|------------|
|  Pseudomonas  |     200     |    100     |
| Streptococcus |      50     |     20     |

当然，你可以计算基因表达数据的这个指标，但这不是常见的方法，我需要更多关于为什么你想这样做的信息。

据我所知，您对在主成分分析图中可视化样本表达之间的距离很感兴趣。使用DESeq2，您可以：

library(DESeq2)
# Get a DESeqDataSet from somewhere
dds <- DESeqDataSetFrom...(...)
# You don't need to run `DESeq()` on the dds for a PCA, just transform your data 
# into a homoscedastic dataset with either VST or rlog
vsd <- varianceStabilizingTransformation(dds, blind=TRUE)
rld <- rlogTransformation(dds, blind=TRUE)
# 'xxx' here takes the place of your condition of interest from your 
# design data frame
plotPCA(vsd, intgroup=c('xxx'))

好吧，假设你想要的是PCA中的基因，而不是样本。在这种情况下，您可以从VST或rlog对象中获取转换后的表达式值，并自己运行PCA代码：

library(DESeq2)
library(ggplot2)
# Get gene expression post VST
vst_expr <- assay(vsd)
# Or - if you want to select some genes
vst_expr <- assay(vsd)[c(...), ]
# Perform PCA
pca <- prcomp(vst_expr)
# Calculate explained % variation
pvar_expl <- round(((pca$sdev ^ 2) / sum(pca$sdev ^ 2)) * 100, 2)
ggplot(as.data.frame(pca$x), aes(x = PC1, y = PC2)) + 
geom_point() +
xlab(paste("PC1: ", pvar_expl[1], "%")) +
ylab(paste("PC2: ", pvar_expl[2], "%"))

最后一点，在进行探索性数据分析之前，通常不建议只选择一些基因，尤其是以你正在思考的方式。你已经测试了这些基因在DESeq2中的差异表达，所以你知道它们是不同的。使用PCA或热图执行盲可视化要好得多。按照此步骤了解有关DESeq2的所有信息，并查看https://support.bioconductor.org/