Доверительные интервалы в PCA по группам с prcomp - PullRequest
Новая
       64

Доверительные интервалы в PCA по группам с prcomp

0 голосов
/ 05 августа 2020

Я работаю над анализом последовательности РНК, и меня интересует, какие гены определяют тканеспецифические c вариации в экспрессии генов. PCAs распространены в анализах РНК-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

1 Ответ

1 голос
/ 05 августа 2020

Ваш dput казался немного большим, но вот воспроизводимый пример - просто добавьте habillage и цвета, и вам будет хорошо go. 

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)

Created on 2020-08-05 by the пакет (v0.3.0)

...