Help for PCA analysis
1
0
Entering edit mode
2.7 years ago
axe880 • 0

Hi everyone,

I am looking to generate some PCA plots for my current project. I am aiming to generate four different plots. This represents the two plant varieties I am using (drought tolerant and drought susceptible) as well as the two tissue types (leaf and root).

I am working with an example script that was provided to me. My issue is that I am unable to successfully edit this script. I am looking to plot all four replicate types (described by the "Time" category for each of the four listed categories above.

Any help would be appreciated.

Here is the script:

# remove genes with low expression
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# transform data before saving output for use later
vsd <- vst(dds, blind = TRUE)

# save
write.csv(as.data.frame(assay(vsd)) %>% rownames_to_column(var = "gene"),
          here("VSDtransformed"),
          row.names = FALSE)

# PCA
pca_data <- plotPCA(vsd, intgroup = c("Tissue", "Time", "Variety"), ntop = 1000, returnData = TRUE)

# Variance Explained
var_capture <- round(100 * attr(pca_data, "percentVar"), 1)

# Plot
p1 <- plot_ly(pca_data, x = ~PC1, y = ~PC2, text = ~paste("Sample:", name, "<br>Variety:", Variety),
              mode = "markers", color = ~Time, symbol = ~Tissue, width = 900)
p1 %>% layout(xaxis =list(title = paste0("PC1", var_capture[1], "% variance")),
              yaxis = list(title = paste0("PC2", var_capture[2], "% variance")))
PCA R • 889 views
ADD COMMENT
3
Entering edit mode

In order to make your post more clean, I suggest you using the format "Code Sample" to format your code lines in a better way. Also, for having better answers, it is usually a good idea to provide some example of the input data, and also, the specific error you are getting when running the code :).

ADD REPLY
1
Entering edit mode
2.7 years ago

I'm not sure why you don't want to put all your samples into one plot, but if you really don't want to do that, here's how you can make a smaller dds object at is a subset of yours. Something like:

dds_drop <- dds[ ,dds$Day == "Day0"]

Then make a vsd for that object, and recalculate PCAData for it.

ADD COMMENT
0
Entering edit mode

Thank you. I essentially have 48 samples:

12 for Tolerant (Root) 12 for Tolerant (Leaf) 12 for Susceptible (Root) 12 for Susceptible (Leaf)

Each 12 represents four time points with three replicates each. I am aiming to separate the dataset into these four groups so I can check the replicates are closely grouped for each time point.

Thanks again

ADD REPLY

Login before adding your answer.

Traffic: 2096 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6