Shape of PC loadings
I'm analyzing RNAseq data and was examining distribution of my PC1 loadings. It looks like half of the genes are ~0, another half at ~0.06, and then expected tail at both ends. Can anyone explain what might cause the step up in the middle of the plot?
Additional information: ~80 samples, analyzed with DESeq vst method below
# Create DESeq object
dds <- DESeqDataSetFromMatrix(countData = count_data, colData = col_data, design = ~ Condition)
# Filter lowly expressed genes
dds <- estimateSizeFactors(dds)
idx <- rowSums( counts(dds, normalized = TRUE) >= 10 ) >= 11
dds <- dds[idx,]
# Run DESeq
dds <- DESeq(dds)
# VST
vsd <- vst(dds, blind = TRUE)
# PCA plot
RNAseqQC::plot_pca(vsd, n_feats = 500)
# PC loadings
top500genes <- order(apply(assay(vsd), 1, var), decreasing = TRUE)[1:500]
top500mat <- assay(vsd)[top500genes, ]
pca <- prcomp(t(top500mat))
loadings_pc1 <- pca$rotation[ ,1]
loadings_pc1 <- sort(loadings_pc1, decreasing = F)
plot(loadings_pc1, xlab = "Rank", ylab = "Loading", main = "PC1")
Thanks!
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Batch and Sample correction for downstream analysis using DESeq2
written by NorbertK •
-
In DESeq2, using plotPCA how do I go about extracting genes making up PC1 and PC2?
written by mgranada3 •
-
Help for PCA analysis
written by axe880 •
-
Constrained rda for a PCA made with DESeq2?
written by luzglongoria •
-
Help me do a gene enrichment analysis
written by anasjamshed
-
Plot biplot (pca loadings) onto PCA plot from DESeq2?
written by NBR •
-
DESeq2 Batch correction enquiry?
written by prithvi.mastermind •
-
Batch effect problem DEG, DESseq2
written by maria2019
-
Correcting for batch effect in RNA-seq data
written by Rimma •
-
Individual genes plots after batch correction with Deseq2
written by Eugene A
The big step from ~275 to 300 might be expected if, within the 500 genes you used for PCA, there’s a subset that’s much more aligned with PC1 than the rest (i.e., a group of genes that "move together" along that axis) against a background of genes with weaker loadings. In my experience, when rank-sorting and plotting signed loadings, this kind of "knee" is not uncommon: it could signify a clump of genes with very similar PC1 weights, and then there's a transition into a different regime of genes.
If you want to check what PC1 is "about," take the subsets of genes (for example, the step, the top bit, the bottom bit, etc.), do a quick GO, enrichment, and/or pathway check (e.g., using ToppGene if these are human data), and then correlate the sample PC1 scores with obvious covariates: batch, size factor/library depth, source study, mitochondrial fraction (if relevant), ribosomal fraction (if relevant), and your actual biological labels (condition, treatment, time point, etc. [whatever is your design for your experiment with 80 samples]). (If you need more details on this, see the
eigencorplotsection in the PCATools docs.) If PC1 tracks batch, size factor, mito fraction, ribo fraction, etc., then it’s probably a technical axis; if it tracks biology and those genes enrich for a coherent program (immune, cell cycle, etc. [that is, if it’s consistent with the biology domain knowledge for your study]), then PC1 is likely capturing biological signal.Also, keep in mind that with this kind of top-variable selection, where you’re looking at a biased subset (the top 500 variable genes), the "knee" can look more dramatic than when you use a different subset of genes or all genes.
Thank you! I'll try it.