NMF/SUM/RPA Pyroseq comparison script

pyrocomparison.R

library(microbiome)
hitchip <- list()
for (method in c("nmf", "sum", "rpa")) {

  tab <- read.profiling(level = "L1", method = method, data.dir = "./", log10 = TRUE)
  colnames(tab) <- toupper(gsub("\\.", "", colnames(tab)))
  hitchip[[method]] <- tab 

}

tab <- read.csv("L1-pyroseq.csv", sep = "\t", row.names = 1)
colnames(tab) <- gsub("\\.", "", toupper(gsub("_.fas", "", gsub("Sort_", "", colnames(tab)))))
hitchip[["pyro"]] <- log10(1 + tab)

# Shared targets
comsr <- intersect(rownames(hitchip[["sum"]]), rownames(hitchip[["pyro"]]))
comsc <- intersect(colnames(hitchip[["sum"]]), colnames(hitchip[["pyro"]]))

pdf("~/L1data-SUM-NMF-RPA-Pyroseq.pdf")
par(mfrow = c(2,2))
for (method in c("sum", "nmf", "rpa")) {
  v1 <- as.vector(unlist(hitchip[[method]][comsr,comsc]))
  v2 <- as.vector(unlist(hitchip[["pyro"]][comsr,comsc]))
  plot(v1, v2, main = paste(method, "vs pyro; corr:", round(cor.test(v1, v2, method = "spearman")$estimate, 3)), ylab = "pyroseq", xlab = method)
}
plot(1)
dev.off()

# -----------------------------------------------------------------------


pdf("~/L1correlations-NMF-RPA-Pyroseq.pdf")
plot(diag(cor(t(hitchip[["pyro"]][comsr,comsc]), t(hitchip[["nmf"]][comsr,comsc]))), diag(cor(t(hitchip[["pyro"]][comsr,comsc]), t(hitchip[["rpa"]][comsr,comsc]))), ylab = "RPA vs. Pyro; Pearson correlations", xlab = "NMF vs. Pyro; Pearson correlations", main = "L1-group correlations NMF/RPA vs. Pyroseq"); abline(0,1)
dev.off()