This is a test version of Biostars. For the public version, visit https://www.biostars.org.
all samples are aligned vertically on the PCA

Hi,

I'm analysing a dataset that does not have great experimental design, and on the PCA all samples align in one vertical line. Does anyone know what could be going on? Thanks!

Edit: I see on the pca that sample_5 seems like an outlier, and this is the only sample sequenced in sequencing run3. However, when removing this sample and rerunning the PCA, the others are still all aligned just like below.

Edit2: I also know that there is perfect confounding between the batch and the genotype

Edit3: first the raw counts coming from Dragen are changed into vst counts which are used for the PCA

model_design <- ~ genotype

dds <- DESeqDataSetFromMatrix(countData = round(raw_data), SamplePlan, design = model_design)
vst <- vst(dds)
avst <- assay(vst)

# plotPCA() is a wrapper we use for FactoMineR::PCA() where scale.unit = TRUE
# tbh I think scale.unit should be FALSE as it's already vst here, see example below

plotPCA(avst, SamplePlan, dims = c(1, 2), VarColor = "genotype",  VarShape = "genotype")

samples:

enter image description here

PCA 1 (vst count for all genes, with PCA(scale.unit = TRUE)):

enter image description here

Edit3:

PCA 2 (vst counts for filtered genes, with PCA(scale.unit = TRUE)):
added filtering of the raw data based on a cpm of at least 1 in at least 2 samples -> PCA looks 'normal' again

raw_data <- raw_data[rowSums(cpm_data >= 1) >= 2,]

enter image description here

PCA 3 (vst counts for all genes, with PCA(scale.unit = FALSE)):
removing the scaling in the PCA function also makes the PCA look normal again -> remaining question: is the weird PCA above due to double-scaling (vst + in PCA function)? should we filter the count data to generate a PCA (intuitively I would say no)?

enter image description here

Edit4: as suggested, I've generated the PCAs with DESeq2::plotPCA to confirm (this functions does not have a scaling option).

PCA4 (vst count for all genes, with DESeq2::plotPCA()) -> very similar to PCA plot above

DESeq2::plotPCA(vst, intgroup = "genotype", ntop = nrow(assay(vst)))

enter image description here

PCA5 (vst count for 500 most variable genes, with DESeq2::plotPCA()) -> more variance explained in PC1, which makes sense as the most variant genes were selected (see also this post and this post).

DESeq2::plotPCA(vst, intgroup = "genotype", ntop = 500)

enter image description here

pca rna-seq
• 161 views • link
updated by ATpoint ATpoint • written by marieke • marieke

Show code please and tell what the inout data are (normalized, log-transformed, scaled?). Almost certainly a code or parsing or input error.

• 0 views • link
written by ATpoint ATpoint

yes of course, my bad for not adding the code. I've added it above, and by looking closer, realised that it might be due to the combination of using vst counts and PCA(scale.unit = TRUE). Does it make sense that 'double scaling' produces this linear PCA?

• 0 views • link
written by marieke • marieke

1 answer

I do not know the FactoMineR package or its functions. To be save, if you want to use an established function, then simply use DESeq2::plotPCA(). That ensures the vst output is handled correctly. Sidenote: Avoid overwriting functions with variables, e.g. vst <- vst(x) -- rather do vsd < vst(x). DESeq2::plotPCA()can directly use vsd. It by default uses the top-500 most variable genes. Maybe in your linear plot it's an unfortunate combination of poor gene selection and scaling, idk.

• 0 views • link
written by ATpoint ATpoint

thanks for all the advise, this helped me a lot! I've added the DESeq2::plotPCA() plots. Based on those, I would think that the issue was mostly the scaling of the vst by the FactoMineR::PCA() function.

• 0 views • link
written by marieke • marieke

Log in to answer this question.

More posts like this