# JMDplots/grad_util.R
# utility functions (not exported) for gradox and gradH2O papers
# get metadata (location names and sequencing IDs) for a study
# extracted from mprep/mplot 20180312
get.mdata <- (studies, study, seqtype, remove.NA = ) {
samples <- studies[study]][1]]
((samples)) ("metadata for ", study, " study not available")
xlabels <- studies[study]]$xlabels
((xlabels)) xlabels <- samples
<- studies[study]]"group"]]
seqtype.for.ID <- seqtype
# change e.g. SRA_MGP to SRA_MG
seqtype.for.ID <- ("_MG.$", "_MG", seqtype.for.ID)
seqtype.for.ID <- ("_MT.$", "_MT", seqtype.for.ID)
IDs <- studies[study]][seqtype.for.ID]]
((IDs)) (seqtype.for.ID, " IDs not available for ", study, " study")
# remove NA IDs and corresponding samples, xlabels, and groups
(remove.NA) {
samples[is.na(IDs)] <-
samples <- (samples)
xlabels[is.na(IDs)] <-
xlabels <- (xlabels)
(() > 1) {
[is.na(IDs)] <-
<- ()
}
IDs <- (IDs)
}
abbrev <- studies[study]]"abbrev"]]
techtype <- studies[study]]"techtype"]]
dx <- studies[study]]"dx"]]
dy <- studies[study]]"dy"]]
((samples=samples, xlabels=xlabels, IDs=IDs, =, abbrev=abbrev, techtype = techtype, dx = dx, dy = dy))
}
# function to plot ZC of metagenomic DNA 20180215
# or ZC of metagenomic proteins 20180228
# NOTE: dataset should end in "_MG" or "_MT" (plot DNA and RNA compositions)
# or "_MGP" or "_MTP" (plot protein compositions)
## optional: vioplotx package is needed for split violin plots
plotMG <- (dataset="Guerrero_Negro_IMG_MG", plottype="bars",
samples=(10:1, =2, flag="0"), =(10:1, =2, flag="0"),
="mat", xlab="layer", =, abbrev=, dsDNA=, plot.RNA=,
taxid=, lwd=1, lty=2, lwd.bars=2, =, =, add.label=,
plot_real_x=, maxdepth=, H2O=, plot.it = , add.title = , yline = 2,
techtype = , dx = , dy = , datadir = ,
add = , all.labels = , = 19, = , srt = 45, ilabel = ) {
# samples: (used for suffixes on file names)
# labels: (used for labeling x-axis ticks)
# xlab: axis label: "layer", "depth", ...
isprotein <- ("_MGP$", dataset) | ("_MTP$", dataset)
# where to keep mean and high/lo (+/- SD) of ZC at each site
pI <- pI.SD <- GRAVY <- GRAVY.SD <- GC <- Xlo <- Xmean <- Xhi <- ()
# are we using user-supplied data?
user_data <-
((datadir)) {
# which paper is this dataset from?
((((), , dataset))) paper <- "gradox"
((((), , dataset))) paper <- "gradH2O"
# gradox or gradH2O data location in JMDplots package 20190928
datadir <- (("extdata/", paper), package = "JMDplots")
user_data <-
}
# set up for proteins or DNA
(isprotein) {
(user_data) filestart <- (datadir, "/MGP/", dataset)
{
rdsdata <- (paper, "_MGP")
filestart <- dataset
}
(H2O) Xfun <- nH2O Xfun <- Zc
(()) <- "darkgreen"
} {
# data directory for DNA
(user_data) filestart <- (datadir, "/MGD/", dataset, "D")
{
rdsdata <- (paper, "_MGD")
filestart <- (dataset, "D")
}
Xfun <- ZCnuc
(()) <- "red"
}
# keep sample values for violin plot 20180515
DNA <- (ZC=((samples)*100), =((samples)*100), stringsAsFactors=)
meancomp <-
# read the data from the RDS file 20191022
(!user_data) RDS <- (rdsdata, )
(i 1:(samples)) {
<- samples[i]
(!(taxid)) <- (filestart, "_", , "-", taxid, ".csv")
<- (filestart, "_", , ".csv")
# use NA if the file is missing 20180529
(user_data) fexists <- ()
fexists <- %in% (RDS)
(!fexists) {
Xmean <- (Xmean, )
Xlo <- (Xlo, )
Xhi <- (Xhi, )
(("missing file", ))
} {
(user_data) mycomp <- (, as.is=)
mycomp <- RDS[file]]
(isprotein) {
# calculate GRAVY 20191024
myGRAVY <- GRAVY(mycomp)
GRAVY <- (GRAVY, (myGRAVY))
GRAVY.SD <- (GRAVY.SD, (myGRAVY))
# calculate isoelectric point 20191027
mypI <- pI(mycomp)
pI <- (pI, (mypI))
pI.SD <- (pI.SD, (mypI))
myX <- Xfun(mycomp)
} {
# use base-paired (double-stranded) DNA
(dsDNA) mycomp <- make_dsDNA(mycomp)
# calculate GC ratio 20180309
(!isprotein) GC <- (GC, GCnuc(mycomp))
myX <- Xfun(mycomp, "deoxyribose")
}
Xmean <- (Xmean, (myX))
Xlo <- (Xlo, (myX) - (myX))
Xhi <- (Xhi, (myX) + (myX))
# initialize data frame for mean compositions
((meancomp)) {
meancomp <- mycomp[1:(samples), ]
meancomp] <-
(meancomp) <- samples
}
# get mean compositions 20180505
meancomp[i, ] <- (mycomp)
# keep sample values for violin plot 20180515
istart <- (i-1) * 100 + 1
iend <- istart + 99
DNA$ZC[istart:iend] <- myX
DNA$[istart:iend] <- ("X", [i])
}
}
# make plot
(plot.it) {
# replace "X" (ZC or nH2O) with GRAVY or pI 20191028
plotXmean <- Xmean
plotXhi <- Xhi
plotXlo <- Xlo
((, "GRAVY")) {
plotXmean <- GRAVY
plotXhi <- GRAVY + GRAVY.SD
plotXlo <- GRAVY - GRAVY.SD
}
((, "pI")) {
plotXmean <- pI
plotXhi <- pI + pI.SD
plotXlo <- pI - pI.SD
}
isdeep <- (())
((all.labels)) all.labels <-
(plot_real_x) {
(!(maxdepth)) isdeep <- > maxdepth
<- (!isdeep])
() <- ()
(!("Bison_Pool", dataset) & !("Menez_Gwen", dataset)) <- ()
atx <-
} {
nsamp <- (all.labels)
<- (1, nsamp)
() <- ()
(plottype=="violin") [1] <- [1] - 0.5
(plottype=="violin") [length()] <- [length()] + 0.5
# for polygon plots: find the correct x-values (may not be 1:nsamp because of missing data) 20191005
atx <- (, all.labels)
}
(!add & ((taxid) | (taxid, 0))) {
(0, 0, =, =, xlab=xlab, ylab=, xaxt="n")
plot.labels <- all.labels
(!(ilabel)) plot.labels-ilabel] <- ""
(!(plot.labels)) {
# rotate labels 20190113
(!(srt)) {
# https://www.r-bloggers.com/rotated-axis-labels-in-r-plots/
# modified to use offset calculated with strheight 20191004
(x=1:nsamp, y=()$usr[3]-1.5*("A"), =plot.labels, srt=srt, adj=1, xpd=)
# add tick marks
(1, at=1:nsamp, =)
} {
# loop over labels so that R doesn't omit any (because they're crowded)
(i 1:nsamp) (1, at=i, =plot.labels[i])
}
} {
# for numeric labels, lower gap.axis to show more labels 20181215 (requires R 3.6.0)
(() >= "3.6.0") (1, at=atx, =plot.labels, gap.axis=0.02)
(1, at=atx, =plot.labels)
}
# add y-axis: ZC (or nH2O 20181231)
((, "GRAVY")) ("GRAVY", side = 2, = yline, las = 0)
((, "pI")) ("pI", side = 2, = yline, las = 0)
(H2O) (((n)[H[2]*O]), side=2, =yline, las=0)
(((Z)[C]), side=2, =yline, las=0)
}
(plottype=="lines") {
(atx!isdeep], plotXlo!isdeep], =, lty=3)
(atx!isdeep], plotXmean!isdeep], =)
(atx!isdeep], plotXhi!isdeep], =, lty=3)
}
# if(plottype=="violin") {
# vioplotx(ZC~sample, DNA, add=TRUE, col = "palevioletred", plotCentre = "line", side = "left", pchMed = 21, colMed = "palevioletred4", colMed2 = "palevioletred2")
# }
(plottype=="bars") {
# error bar plots 20180515
# apply small offset to x-position to separate DNA and RNA
(!(taxid)) xx <- 0
xx <- ((("usr")[1:2])) / 120
# bars (whiskers) at one SD from mean (arrows trick from https://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars)
(atx-xx, plotXlo, atx-xx, plotXhi, = 0.03, angle = 90, code = 3, =, lwd=lwd.bars)
# remove NAs so we can draw lines between all sites 20180529
iNA <- (plotXmean)
((atx)!iNA & !isdeep]-xx, plotXmean!iNA & !isdeep], =, lty=lty, lwd=lwd)
}
((plottype, 1, 1)=="#") {
# polygon plots 20191005
# assemble top and bottom coordinates of polygon
((atx, (atx)), (plotXhi, (plotXlo)), = plottype, border = )
# remove NAs so we can draw lines between all sites 20180529
iNA <- (plotXmean)
((atx)!iNA & !isdeep], plotXmean!iNA & !isdeep], =, lty=lty, lwd=lwd)
((atx)!iNA & !isdeep], plotXmean!iNA & !isdeep], = , =)
}
# add points to show > 1% species abundance 20181118
(!(taxid) & !(taxid, 0)) {
# gradox data location in JMDplots package 20190928
datadir <- ("extdata/gradox", package = "JMDplots")
<- (datadir, "/one_percent/", dataset, ".csv")
dat <- ()
# which samples have this taxid with at least 1% abundance?
dat <- dat[dat$taxid==taxid, ]
dat <- dat[dat$percentage >= 1, ]
ioneperc <- samples %in% dat$
(atx[ioneperc], plotXmean[ioneperc], =, cex=0.5, =)
}
}
# now do it for RNA
(!isprotein) {
RNA_Xlo <- RNA_Xmean <- RNA_Xhi <- ()
# keep sample values for violin plot 20180515
RNA <- (ZC=((samples)*100), =((samples)*100), stringsAsFactors=)
# gradox or gradH2O data location in JMDplots package 20190928
(!user_data) {
datadir <- (("extdata/", paper), package = "JMDplots")
rdsdata <- (paper, "_MGR")
RDS <- (rdsdata, )
filestart <- (dataset, "R")
}
(i 1:(samples)) {
<- samples[i]
(user_data) {
<- (datadir, "/MGR/", dataset, "R_", , ".csv")
fexists <- ()
} {
<- (filestart, "_", , ".csv")
fexists <- %in% (RDS)
}
(fexists) {
(user_data) myRNA <- (, as.is=)
myRNA <- RDS[file]]
myX <- Xfun(myRNA, "ribose")
RNA_Xmean <- (RNA_Xmean, (myX))
RNA_Xlo <- (RNA_Xlo, (myX) - (myX))
RNA_Xhi <- (RNA_Xhi, (myX) + (myX))
}
# keep sample values for violin plot 20180515
istart <- (i-1) * 100 + 1
iend <- istart + 99
RNA$ZC[istart:iend] <- myX
RNA$[istart:iend] <- ("X", [i])
}
(plot.RNA & (RNA_Xmean) > 0 & plot.it) {
# subtract offset for ZC of RNA
dZC <- -0.28
(plottype=="lines") {
(atx, RNA_Xlo + dZC, ="blue", lty=3)
(atx, RNA_Xmean + dZC, ="blue")
(atx, RNA_Xhi + dZC, ="blue", lty=3)
}
# if(plottype=="violin") {
# RNA$ZC <- RNA$ZC + dZC
# vioplotx(ZC~sample, RNA, add=TRUE, col = "lightblue", plotCentre = "line", side = "right", pchMed = 21, colMed = "lightblue4", colMed2 = "lightblue2")
# }
(plottype=="bars") {
# apply small offset to x-position to separate DNA and RNA
xx <- ((("usr")[1:2])) / 200
(atx+xx, RNA_Xlo + dZC, atx+xx, RNA_Xhi + dZC, = 0.03, angle = 90, code = 3, ="blue", lwd=lwd.bars)
(atx+xx, RNA_Xmean + dZC, ="blue", lty=2, lwd=lwd)
}
}
# return ZC values
outval <- (DNA=Xmean, RNA=RNA_Xmean, GC=GC, =, meancomp=meancomp, abbrev=abbrev, techtype = techtype, dx = dx, dy = dy, H2O = H2O)
} {
outval <- (AA=Xmean, GRAVY=GRAVY, pI=pI, =, meancomp=meancomp, abbrev=abbrev, techtype = techtype, dx = dx, dy = dy, H2O = H2O)
}
# add title 20180225
(((taxid) | (taxid, 0)) & plot.it & add.title & !add) {
main <- dataset2main(dataset, abbrev)
(main=main, font.main=1)
# add in-plot label: MG or MT 20180829
(add.label) {
label <- "MG"
(("_MT.*", dataset)) label <- "MT"
# change position for some datasets
xfrac <- 0.1
(("OMZ", dataset)) xfrac <- 0.9
(("Negro", dataset)) xfrac <- 0.9
(("ALOHA", dataset)) xfrac <- 0.9
CHNOSZ::label.plot(label, xfrac=xfrac, yfrac=0.9)
}
}
# done!
(outval)
}
# parse dataset name to plot title 20180505
dataset2main <- (dataset, abbrev=) {
main <- dataset
main <- ("_SRA", "", main)
main <- ("_IMG", "", main)
main <- ("_MGRAST", "", main)
main <- ("_GenBank", "", main)
main <- ("_iMicrobe", "", main)
# strip _MG or _MGP
main <- ("_MG.*", "", main)
main <- ("_MT.*", "", main)
# replace underscore with space, and put space in Baltic Sea
main <- ("_", " ", main)
main <- ("BalticSea", "Baltic Sea", main)
main <- ("Mud Volcano", "SYNH Mud Volcano", main)
# add dataset abbreviation at end 20180829
((abbrev, "GS")) {
# replace "GS" (Global Survey of Baltic Sea) by "MG" or "MT"
(("_MTP", dataset)) abbrev <- "MT"
(("_MGP", dataset)) abbrev <- "MG"
}
(!(abbrev)) main <- (main, " (", abbrev, ")")
main
}
### utility functions for ZC, GC and nH2O calculations ###
# calculate GC ratio for given nucleobase compositions 20180309
GCnuc <- (nuccomp) {
# find columns with names for the nucleobases
isbase <- (nuccomp) %in% ("A", "C", "G", "T", "U")
# keep only these columns
nuccomp <- nuccomp, isbase]
# find GC columns and not-GC columns
isGC <- (nuccomp) %in% ("G", "C")
isnotGC <- (nuccomp) %in% ("A", "T", "U")
# calculate GC ratio
(nuccomp, isGC]) / (nuccomp)
}
# calculate ZC for given nucleobase compositions 20171223
# changed to nucleosides 20180512
ZCnuc <- (nuccomp, sugar="deoxyribose") {
## ZC and nC of nucleobases
#ZC_nuc <- c(A=2, C=1.5, G=2.4, T=0.8, U=1.5)
#nC_nuc <- c(A=5, C=4, G=5, T=5, U=4)
# the number of carbons of the nucleosides (same in DNA and RNA)
nC_nuc <- (A=10, =9, G=10, T=10, U=9)
(sugar=="deoxyribose") {
# ZC of the nucleosides in DNA
# DNAsides <- c("deoxyadenosine", "deoxycytidine", "deoxyguanosine", "deoxythymidine", "deoxyuridine")
# ZC_nuc <- ZC(info(info(DNAsides))$formula)
ZC_nuc <- (A=0.8, =4/9, G=1, T=0.2, U=4/9)
}
(sugar=="ribose") {
# ZC of the nucleosides in RNA
# RNAsides <- c("adenosine", "cytidine", "guanosine", "thymidine", "uridine")
# ZC_nuc <- ZC(info(info(RNAsides))$formula)
ZC_nuc <- (A=1, =2/3, G=1.2, T=0.4, U=2/3)
}
# find columns with names for the nucleosides
isbase <- (nuccomp) %in% ("A", "C", "G", "T", "U")
ibase <- ((nuccomp)[isbase], (ZC_nuc))
# calculate the nC from the counts
multC <- ((nuccomp, isbase]) * nC_nuc[ibase])
# multiply nC by ZC to get total formal charge on carbon (Z)
multZ <- ((multC) * ZC_nuc[ibase])
# calculate the Z and nC in each composition (row)
Ztot <- (multZ)
nCtot <- (multC)
# the ZC in each composition (row)
Ztot / nCtot
}
# function to convert ssDNA base counts to dsDNA 20180318
make_dsDNA <- (nuccomp) {
nuccomp, ("A", "C", "G", "T")] <- nuccomp, ("A", "C", "G", "T")] + nuccomp, ("T", "G", "C", "A")]
nuccomp
}
