DESeq - Changing PCA Axis to PC3 and PC4.

[Chequer]

Hi all, I'm new to R and I'm currently going through the DESeq introduction guide by Anders and Huber entitled 'Differential expression of RNA-Seq data at the gene level.'

Simply, I have created a PCA plot using:
print(plotPCA(vsdFull, intgroup=c("condition", "libType")))

My question is how do I produce a graph showing other eigenvectors like PC3 and PC4, ideally using DESeq and gplots.

Thanks in advance.

[fanli]

You can modify/rewrite the existing plotPCA function. Bolded sections are modified to include PC3 and PC4:

> DESeq2:::plotPCA.DESeqTransform
function (object, intgroup = "condition", ntop = 500, returnData = FALSE) 
{
    rv <- rowVars(assay(object))
    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, 
        length(rv)))]
    pca <- prcomp(t(assay(object)[select, ]))
    percentVar <- pca$sdev^2/sum(pca$sdev^2)
    if (!all(intgroup %in% names(colData(object)))) {
        stop("the argument 'intgroup' should specify columns of colData(dds)")
    }
    intgroup.df <- as.data.frame(colData(object)[, intgroup, 
        drop = FALSE])
    group <- if (length(intgroup) > 1) {
        factor(apply(intgroup.df, 1, paste, collapse = " : "))
    }
    else {
        colData(object)[[intgroup]]
    }
    d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], [B]PC3 = pca$x[, 3], PC4 = pca$x[, 4],[/B] group = group, 
        intgroup.df, name = colnames(object))
    if (returnData) {
        attr(d, "percentVar") <- percentVar[1:2]
        return(d)
    }
    ggplot(data = d, aes_string(x = "[B]PC3[/B]", y = "[B]PC4[/B]", color = "group")) + 
        geom_point(size = 3) + xlab(paste0("[B]PC3[/B]: ", round(percentVar[[B]3[/B]] * 
        100), "% variance")) + ylab(paste0("[B]PC4[/B]: ", round(percentVar[[B]4[/B]] * 
        100), "% variance")) + coord_fixed()
}

[Chequer]

You can modify/rewrite the existing plotPCA function. Bolded sections are modified to include PC3 and PC4:

> DESeq2:::plotPCA.DESeqTransform
function (object, intgroup = "condition", ntop = 500, returnData = FALSE) 
{
    rv <- rowVars(assay(object))
    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, 
        length(rv)))]
    pca <- prcomp(t(assay(object)[select, ]))
    percentVar <- pca$sdev^2/sum(pca$sdev^2)
    if (!all(intgroup %in% names(colData(object)))) {
        stop("the argument 'intgroup' should specify columns of colData(dds)")
    }
    intgroup.df <- as.data.frame(colData(object)[, intgroup, 
        drop = FALSE])
    group <- if (length(intgroup) > 1) {
        factor(apply(intgroup.df, 1, paste, collapse = " : "))
    }
    else {
        colData(object)[[intgroup]]
    }
    d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], [B]PC3 = pca$x[, 3], PC4 = pca$x[, 4],[/B] group = group, 
        intgroup.df, name = colnames(object))
    if (returnData) {
        attr(d, "percentVar") <- percentVar[1:2]
        return(d)
    }
    ggplot(data = d, aes_string(x = "[B]PC3[/B]", y = "[B]PC4[/B]", color = "group")) + 
        geom_point(size = 3) + xlab(paste0("[B]PC3[/B]: ", round(percentVar[[B]3[/B]] * 
        100), "% variance")) + ylab(paste0("[B]PC4[/B]: ", round(percentVar[[B]4[/B]] * 
        100), "% variance")) + coord_fixed()
}

Unfortunately on using that code, I get this error:

> print(plotPCA(vsdFull, intgroup=c("condition", "libType")))
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘plotPCA’ for signature ‘"ExpressionSet"’

I understand that I need to use an argument but simply what code do I need to make it work. I do apologise for being so new to R.

[fanli]

Hmm, I should have been more specific. I think if you copy and paste the following code into your R session it should work:

plotPCA <- function (object, intgroup = "condition", ntop = 500, returnData = FALSE) 
{
    rv <- rowVars(assay(object))
    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, 
        length(rv)))]
    pca <- prcomp(t(assay(object)[select, ]))
    percentVar <- pca$sdev^2/sum(pca$sdev^2)
    if (!all(intgroup %in% names(colData(object)))) {
        stop("the argument 'intgroup' should specify columns of colData(dds)")
    }
    intgroup.df <- as.data.frame(colData(object)[, intgroup, 
        drop = FALSE])
    group <- if (length(intgroup) > 1) {
        factor(apply(intgroup.df, 1, paste, collapse = " : "))
    }
    else {
        colData(object)[[intgroup]]
    }
    d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], PC3 = pca$x[, 3], PC4 = pca$x[, 4], group = group, 
        intgroup.df, name = colnames(object))
    if (returnData) {
        attr(d, "percentVar") <- percentVar[1:2]
        return(d)
    }
    ggplot(data = d, aes_string(x = "PC3", y = "PC4", color = "group")) + 
        geom_point(size = 3) + xlab(paste0("PC3: ", round(percentVar[1] * 
        100), "% variance")) + ylab(paste0("PC4: ", round(percentVar[2] * 
        100), "% variance")) + coord_fixed()
}

print(plotPCA(vsdFull, intgroup=c("condition", "libType")))

By the way, just to verify, are you using DESeq2 or DESeq?

[Chequer]

Hmm, I should have been more specific. I think if you copy and paste the following code into your R session it should work:

plotPCA <- function (object, intgroup = "condition", ntop = 500, returnData = FALSE) 
{
    rv <- rowVars(assay(object))
    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, 
        length(rv)))]
    pca <- prcomp(t(assay(object)[select, ]))
    percentVar <- pca$sdev^2/sum(pca$sdev^2)
    if (!all(intgroup %in% names(colData(object)))) {
        stop("the argument 'intgroup' should specify columns of colData(dds)")
    }
    intgroup.df <- as.data.frame(colData(object)[, intgroup, 
        drop = FALSE])
    group <- if (length(intgroup) > 1) {
        factor(apply(intgroup.df, 1, paste, collapse = " : "))
    }
    else {
        colData(object)[[intgroup]]
    }
    d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], PC3 = pca$x[, 3], PC4 = pca$x[, 4], group = group, 
        intgroup.df, name = colnames(object))
    if (returnData) {
        attr(d, "percentVar") <- percentVar[1:2]
        return(d)
    }
    ggplot(data = d, aes_string(x = "PC3", y = "PC4", color = "group")) + 
        geom_point(size = 3) + xlab(paste0("PC3: ", round(percentVar[1] * 
        100), "% variance")) + ylab(paste0("PC4: ", round(percentVar[2] * 
        100), "% variance")) + coord_fixed()
}

print(plotPCA(vsdFull, intgroup=c("condition", "libType")))

By the way, just to verify, are you using DESeq2 or DESeq?

I'm using DESeq but I also have DESeq2 installed. Unfortunately, the same error is returned. I'm working through the DESeq guide and once I use this new code, my original PCA fails to work at all with the error. I massively appreciate your help. I just feel as if I need to learn the code a bit more to understand why its not working.

[fanli]

Oh, sorry. I was basing everything on DESeq2. Try this:

plotPCA <- function (x, intgroup = "condition", ntop = 500, ax1="PC1", ax2="PC2") 
{
    rv = rowVars(exprs(x))
    select = order(rv, decreasing = TRUE)[seq_len(ntop)]
    pca = prcomp(t(exprs(x)[select, ]))
    fac = factor(apply(pData(x)[, intgroup, drop = FALSE], 1, 
        paste, collapse = " : "))
    colours = if (nlevels(fac) >= 3) 
        brewer.pal(nlevels(fac), "Paired")
    else c("green", "blue")
    xyplot(as.formula(sprintf("%s ~ %s", ax2, ax1)), groups = fac, data = as.data.frame(pca$x), 
        pch = 16, cex = 2, aspect = "iso", col = colours, main = draw.key(key = list(rect = list(col = colours), 
            text = list(levels(fac)), rep = FALSE)))
}

plotPCA(vsdFull, intgroup=c("condition", "libType"), ax1="PC3", ax2="PC4")

[Chequer]

Oh, sorry. I was basing everything on DESeq2. Try this:

plotPCA <- function (x, intgroup = "condition", ntop = 500, ax1="PC1", ax2="PC2") 
{
    rv = rowVars(exprs(x))
    select = order(rv, decreasing = TRUE)[seq_len(ntop)]
    pca = prcomp(t(exprs(x)[select, ]))
    fac = factor(apply(pData(x)[, intgroup, drop = FALSE], 1, 
        paste, collapse = " : "))
    colours = if (nlevels(fac) >= 3) 
        brewer.pal(nlevels(fac), "Paired")
    else c("green", "blue")
    xyplot(as.formula(sprintf("%s ~ %s", ax2, ax1)), groups = fac, data = as.data.frame(pca$x), 
        pch = 16, cex = 2, aspect = "iso", col = colours, main = draw.key(key = list(rect = list(col = colours), 
            text = list(levels(fac)), rep = FALSE)))
}

plotPCA(vsdFull, intgroup=c("condition", "libType"), ax1="PC3", ax2="PC4")

Worked perfectly! Thank you ever so much!