1
$\begingroup$

I'm a biologist, not a programmer so please be gentle.

So I have a dataset that looks like

Genes  Patient1   Patient2   Patient3
A          324      433         343
B          431       342        124
Z          232       234        267

then I have the sample sheet where it contains sample info like:


Patient1 - Healthy
Patient2 - Disease
Patient3 - Healthy 
(110 samples)

I am using:

library(ggfortify)
df <- dataset
pca_res <- prcomp(df, scale. = TRUE)

autoplot(pca_res)

Then I want to do

autoplot(pca_res, data = ?, colour = '?')

I wish to use the info from the sample sheet to colour my PCA based on the state (healthy/disease) using the autoplot function. Is there a way to do this?

$\endgroup$
2
  • 1
    $\begingroup$ I think you need to provide more information for a minimal reproducible example. $\endgroup$
    – zorbax
    Jun 9, 2021 at 7:58
  • $\begingroup$ Given that you probably do not want to invest lots of time in such a simple task I'd just read the PCAtools vignette over at Bioconductor. That package can probably do all you want including making pretty plots, and the heavy work is coded under the hood so for the end user it comes down to running a couple of commands in R. $\endgroup$
    – user3051
    Jun 9, 2021 at 12:50

2 Answers 2

1
$\begingroup$

I've found this the way to work with autoplot (ggfortify's prcomp autoplot) - Use a df that has the field you need for color as the first column (to make things easy), then do:

autoplot(prcomp(df[,-1], data=df, color='column1_name')

If the color-by field is not the first column, do:

autoplot(prcomp(df[,!(colnames(df) == "color_by_column_name")], data = df, color = "color_by_column_name")

EDIT:

That being said, you're going to need to process your data so patients are along rows instead of column. Transpose your dataset so rows are columns and columns become rows; change patient ID to rownames and Gene name to colnames, then merge the patient category dataset with the expression dataset. Once you merge, your dataset should look like:

         Category    A      B      Z
Patient1 Healthy     324    431    232
Patient2 Disease     433    342    124
Patient3 Healthy     343    124    267

And then,

autoplot(prcomp(df[,!(colnames(df) == "Category")]), data = df, color = "Category")
$\endgroup$
3
  • $\begingroup$ Hi Ram,Thanks a lot for taking the time to answer me. I did exactly what you told me above and there is an error Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric Any ideas why is that ? I am not sure at all $\endgroup$
    – Athon
    Jun 14, 2021 at 14:08
  • $\begingroup$ Ensure that the columns you're passing to prcomp are all numeric. That's probably what's going wrong here. $\endgroup$
    – Ram RS
    Jun 14, 2021 at 23:08
  • $\begingroup$ sorted, thanks a lot :) $\endgroup$
    – Athon
    Jun 14, 2021 at 23:11
2
$\begingroup$

Can't say I've ever used autoplot, but this is fairly easy to achieve using "base" ggplot.

One thing autoplot does is put the axis on the same scale reflecting the percentage variation captured by each PC - I think this is a good thing, but most PCA plots you'll see don't do this.

As you haven't given us any info on your metadata, I'll just make a few assumptions - lets make it a data.frame called "colData" with a "Status" column.

prcomp stores the PC scores in "pca_res$x"

# generate some metadata
colData <- data.frame(row.names=rownames(pca_res$x),status=as.factor(sample(c("diseased","healthy"),nrow(pca_res$x),replace = T)))

# get % variation for each PC 
pca_res$percentVar <- pca_res$sdev^2/sum(pca_res$sdev^2)

# change scale to reflect % variation
d <- t(data.frame(t(pca_res$x)*pca_res$percentVar))

# merge metadata and modified PC scores by row.names (by=0)
d <- merge(d,colData,by=0) 

# make the base plot (PC1 vs PC2)
g <- ggplot(d,aes(x=PC1,y=PC2,colour=status))
# add the points
g <- g + geom_point(size=2)
# add a "wow" factor to the points
geom_point(size=4,alpha=0.5)
# display the plot
g 

There may be a way to do this with autoplot - e.g.

autoplot(pca_res,colour=as.numeric(colData$status))

gets you part of the way there, but I think using the standard ggplot syntax is easy enough anyway.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.