Skip to content

Instantly share code, notes, and snippets.

@sgibb
Created September 19, 2016 08:24
Show Gist options
  • Save sgibb/068ad4028b3e23e5ec2bb26a0d667c32 to your computer and use it in GitHub Desktop.
Save sgibb/068ad4028b3e23e5ec2bb26a0d667c32 to your computer and use it in GitHub Desktop.
plot proteomes in a similar manner as Kuharev et al 2015 10.1002/pmic.201400396
library("synapter")
library("MSnbase")
load_all()
combMSnSet <- readRDS("refCombMSNSet.RDS")
plotProteomes <- function(x,
xlim=log2(range(exprs(x), na.rm=TRUE)),
ylim=c(-4, 4),
xlab=expression(log[2](B)),
ylab=expression(log[2](A:B)),
mapping=factor(rep(1:2, each=5), labels=c("A", "B")),
ratios=c("Homo sapiens"=1,
"Saccharomyces cerevisiae"=2,
"Escherichia coli"=1/4,
"misc"=NA),
col=c("Homo sapiens"="#1b9e77",
"Saccharomyces cerevisiae"="#d95f02",
"Escherichia coli"="#7570b3",
"misc"="#e6ab02"),
pch=c("Homo sapiens"=20,
"Saccharomyces cerevisiae"=20,
"Escherichia coli"=20,
"misc"=18),
pcex=0.75,
alpha=255/2,
f=2/3, ## lowess span
st=Inf, ## saturation threshold
...){
## fetch protein description
pd <- as.matrix(fData(x)[, grep("protein\\.Description", fvarLabels(x))])
## fetch organism name
sp <- gsub(".*OS=([A-z]+ [A-z]+).*", "\\1", pd)
## a few proteins are not identified (and contain NA)
## or are identified as different species
## we choose the "mean" species
species <- sp[,1]
## misidentified?
mi <- rowSums(sp[,1] != sp)
mi <- is.na(mi) | mi > 0
species[mi] <- apply(sp[mi,, drop=FALSE], 1, function(xx)names(which.max(table(xx)))[1])
.rowMeansGroups <- function(x, groups) {
na <- !is.na(x)
mode(na) <- "numeric"
na <- t(rowsum(t(na), groups, na.rm=TRUE))
rs <- t(rowsum(t(x), groups, na.rm=TRUE))
rs/na
}
m <- .rowMeansGroups(exprs(x), mapping)
m[, 1] <- m[,1]/m[,2]
colnames(m)[1] <- "A/B"
m <- log2(m)
## misc colour for species different from Humans, Yeast or Ecoli
scol <- col[species]
scol[!species %in% names(col)] <- col["misc"]
scol <- paste0(scol, as.hexmode(round(alpha)))
## misc pch for species different from Humans, Yeast or Ecoli
spch <- pch[species]
spch[!species %in% names(pch)] <- pch["misc"]
plot(NA, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, ...)
abline(h=log2(ratios), col=col[names(ratios)], lty=2)
points(m[,2], m[,1], col=scol, pch=spch, cex=pcex)
legend("topright", legend=names(col), col=col, pch=pch, cex=pcex, bty="n")
abline(v=log2(st), lty=2, col="#808080")
text(x=log2(st), y=-4,
labels=paste0("saturation threshold = ", prettyNum(st)),
col="#808080", cex=pcex, pos=ifelse(mean(xlim) < log2(st), 2, 4))
if (!is.na(f) && !is.null(f)) {
for (nm in names(col)[names(col) != "misc"]) {
s <- m[species == nm, 2:1]
s <- s[rowSums(is.na(s) | is.infinite(s)) == 0L,]
lines(lowess(s, f=f), col=col[nm], lty=6, lwd=2)
}
}
}
msnSetSum <- requantify(combMSnSet, method="sum", saturationThreshold=Inf, onlyCommonIsotopes=FALSE)
msnSetRef <- requantify(combMSnSet, method="reference", saturationThreshold=Inf)
msnSetTh <- requantify(combMSnSet, method="th.mean", saturationThreshold=Inf, requantifyAll=FALSE)
xl <- c(10, 22)
pcex <- 0.5
pdf("compare.pdf", width=12, height=7, paper="a4r")
par(mfrow=c(2, 2))
plotProteomes(combMSnSet, main="none", xlim=xl, pcex=pcex)
plotProteomes(msnSetSum, main="sum", xlim=xl, pcex=pcex)
plotProteomes(msnSetRef, main="reference", xlim=xl, pcex=pcex)
plotProteomes(msnSetTh, main="th.mean", xlim=xl, pcex=pcex)
par(mfrow=c(1, 1))
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment