使用 prcomp 分组的 PCA 置信区间答案

【问题标题】：Confidence intervals in PCA by group with prcomp使用 prcomp 分组的 PCA 置信区间
【发布时间】：2020-08-05 13:23:51
【问题描述】：

我正在进行 RNA-seq 分析，我对驱动基因表达的组织特异性变异的基因感兴趣。 PCAs 在 RNA-seq 分析中很常见，但大多数软件包（例如 DESeq2）仅将其用于 2D 图。因此，我使用了 prcomp 和 fviz_pca_ind 来生成我的 2D 图，这样我就可以进一步分析。但是，由于输入数据结构，我似乎无法将我的分组信息（即组织）包含在我的 prcomp 对象中，因此无法对我的样本进行颜色编码或在我的组周围产生置信区间。

包：

library(ggplot2)
library(factoextra)
library(Deseq2)

我从一个 DESeq2 对象的方差稳定转换开始（很抱歉这是多长时间......）：

即使是子集，数据集也太大 - 这是一个链接（希望这是可以接受的） https://drive.google.com/file/d/1Gtw5GUCAyBVr3MI6CpgsF4n81KZuq8WI/view?usp=sharing

还有我的 PCA 代码（这实际上只是 DESeq2 中的函数）

####################################################################################################
# Principle component analysis - PRComp
####################################################################################################

# set number of genes to include in PCA
ntop = 500

# calculate the variance for each gene
rv <- rowVars(assay(vst))

# select the ntop genes by variance
select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]

# perform a PCA on the data in assay(x) for the selected genes
tissue_pca <- prcomp(t(assay(vst)[select,]))

还有我的可视化代码：

# Visualize samples on PC1 and PC2
fviz_pca_ind(tissue_pca,
            # col.var = ,
             palette = cbPalette,
             ellipse.type = "confidence",
             repel = TRUE,
             mean = FALSE)

产生这个：

正如您将在 prcomp 对象中看到的，没有分组信息。关于如何 1) 将此信息包含在 prcomp 对象中或 2) 将此信息合并到图形代码中的任何想法？fa

【问题讨论】：

标签： r pca prcomp

【解决方案1】：

您的 dput 似乎有点大，但这里有一个可重现的示例 - 只需添加 habillage 和颜色即可。

suppressPackageStartupMessages(invisible(
  lapply(c("DESeq2", "factoextra"),
         require, character.only = TRUE)))
set.seed(123)
dds <- makeExampleDESeqDataSet(n=5e3, betaSD=.8)
dds$group <- gl(4,3, labels = paste0("Group_", LETTERS[1:4]))
counts(dds)[, 1:3] <- as.integer(counts(dds)[, 1:3] *
  stats::runif(1.5e4, 0.8, 1.2))
vst <- vst(dds)
ntop = 500
rv <- rowVars(assay(vst))
select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]
tissue_pca <- prcomp(t(assay(vst)[select,]))

fviz_pca_ind(tissue_pca, 
  habillage = vst$group,
  palette = c("blue", "red", "cyan", "magenta"),
  ellipse.type = "confidence",
  repel = TRUE,
  mean = FALSE)

^{由reprex package (v0.3.0) 于 2020 年 8 月 5 日创建}

【讨论】：