### R code from vignette source 'Cheminfo.Rnw' ################################################### ### code chunk number 1: Cheminfo.Rnw:865-866 (eval = FALSE) ################################################### ## download.file(url="http://faculty.ucr.edu/~tgirke/HTML_Presentations/Manuals/Workshop_Dec_12_16_2013/Rcheminfo/Cheminfo.R", destfile="Cheminfo.R") ################################################### ### code chunk number 2: Cheminfo.Rnw:870-871 ################################################### library("ChemmineR") # Loads the package ################################################### ### code chunk number 3: Cheminfo.Rnw:874-877 (eval = FALSE) ################################################### ## library(help="ChemmineR") # Lists all functions and classes ## vignette("ChemmineR") # Opens PDF manual from R ## ?MW # Opens help for MW function ################################################### ### code chunk number 4: Cheminfo.Rnw:880-884 (eval = FALSE) ################################################### ## sdfset <- read.SDFset("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/Samples/sdfsample.sdf") ## sdfset # Returns summary of SDFset ## valid <- validSDF(sdfset) # Identifies invalid SDFs in SDFset objects ## sdfset <- sdfset[valid] # Removes invalid SDFs, if there are any ################################################### ### code chunk number 5: Cheminfo.Rnw:887-890 ################################################### data(sdfsample) sdfset <- sdfsample sdfset # Returns summary of SDFset ################################################### ### code chunk number 6: Cheminfo.Rnw:900-902 ################################################### write.SDF(sdfset[1:4], file="sub.sdf", sig=TRUE) list.files(pattern="sub.sdf") ################################################### ### code chunk number 7: Cheminfo.Rnw:906-907 ################################################### write.SDFsplit(x=sdfset, filetag="myfile", nmol=50) ################################################### ### code chunk number 8: Cheminfo.Rnw:911-913 ################################################### sdfsetsub <- read.SDFset("sub.sdf") sdfsetsub ################################################### ### code chunk number 9: Cheminfo.Rnw:923-927 ################################################### data(smisample); smiset <- smisample smiset smiset[[1]] view(smiset[1:2]) ################################################### ### code chunk number 10: Cheminfo.Rnw:931-933 ################################################### write.SMI(smiset[1:4], file="sub.smi", cid=TRUE) write.SMI(smiset[1:4], file="sub.smi", cid=FALSE) ################################################### ### code chunk number 11: Cheminfo.Rnw:937-940 (eval = FALSE) ################################################### ## library(ChemmineOB) # Requries OpenBabel ## mysdf <- smiles2sdf(smiset) ## mysmi <- sdf2smiles(mysdf) ################################################### ### code chunk number 12: Cheminfo.Rnw:957-961 (eval = FALSE) ################################################### ## p450 <- read.SDFset("p450.sdf") ## length(p450) ## write.SDF(p450[length(p450):1], file="inverse.sdf", sig=TRUE) ## write.SDFsplit(x=p450, filetag="p450", nmol=5) ################################################### ### code chunk number 13: Cheminfo.Rnw:1008-1014 (eval = FALSE) ################################################### ## view(sdfset[1:4]) # Summary view of several molecules ## length(sdfset) # Returns number of molecules ## sdfset[[1]] # Returns single molecule from SDFset as SDF object ## sdfset[[1]][[2]] # Returns atom block from first compound as matrix ## sdfset[[1]][[2]][1:4,] ## c(sdfset[1:4], sdfset[5:8]) # Concatenation of several SDFsets ################################################### ### code chunk number 14: Cheminfo.Rnw:1018-1020 (eval = FALSE) ################################################### ## grepSDFset("650001", sdfset, field="datablock", mode="subset") ## # To return index, set mode="index") ################################################### ### code chunk number 15: Cheminfo.Rnw:1030-1037 (eval = FALSE) ################################################### ## sdf <- sdfset[[1]] ## atomblock(sdf); sdf[[2]]; sdf[["atomblock"]] ## # All three methods return the same component ## header(sdfset[1:4]) ## atomblock(sdfset[1:4]) ## bondblock(sdfset[1:4]) ## datablock(sdfset[1:4]) ################################################### ### code chunk number 16: Cheminfo.Rnw:1041-1048 ################################################### sdfid(sdfset[1:4]) # Retrieves CMP IDs from Molecule Name field in header block. cid(sdfset[1:4]) # Retrieves CMP IDs from ID slot in SDFset. unique_ids <- makeUnique(sdfid(sdfset)) # Creates unique IDs by appending a counter to duplicates. cid(sdfset) <- unique_ids # Assigns uniquified IDs to ID slot ################################################### ### code chunk number 17: Cheminfo.Rnw:1063-1074 (eval = FALSE) ################################################### ## ## Task 1 ## cid(sdfset)[1:4] <- month.name[1:4] ## write.SDF(sdfset[1:4], file="sub.sdf", cid=TRUE) ## ## Task 2 ## bondList <- bondblock(sdfset) ## bondMA <- do.call("rbind", bondList) ## dim(bondMA) ## write.table(bondMA, file="bondMA.txt", quote=FALSE, sep="\t") ## ## Task 3 ## atomblock(sdfset)[1] <- atomblock(sdfset)[2] ## atomblock(sdfset[[1]]) == atomblock(sdfset[[1]]) ################################################### ### code chunk number 18: struct1 ################################################### data(sdfsample); sdfset <- sdfsample plot(sdfset[1:4], print=FALSE) # print=TRUE returns SDF summaries ################################################### ### code chunk number 19: struct2 ################################################### plot(sdfset["CMP1"], atomnum = TRUE, noHbonds=F , no_print_atoms = "", atomcex=0.8, sub=paste("MW:", MW(sdfsample["CMP1"])), print=FALSE) ################################################### ### code chunk number 20: struct3 ################################################### plot(sdfset[1], print=FALSE, colbonds=c(22,26,25,3,28,27,2,23,21,18,8,19,20,24)) ################################################### ### code chunk number 21: Cheminfo.Rnw:1130-1131 (eval = FALSE) ################################################### ## sdf.visualize(sdfset[1:4]) ################################################### ### code chunk number 22: Cheminfo.Rnw:1152-1157 (eval = FALSE) ################################################### ## ## Task 1 ## cid(p450) <- sdfid(p450) ## plot(p450[c("42631481", "42631375", "42631371", "42631260")]) ## ## Task 2 ## sdf.visualize(p450[c("42631481", "42631375", "42631371", "42631260")]) ################################################### ### code chunk number 23: Cheminfo.Rnw:1169-1171 ################################################### propma <- atomcountMA(sdfset, addH=FALSE) propma[1:4,] ################################################### ### code chunk number 24: Cheminfo.Rnw:1174-1176 ################################################### data(atomprop) atomprop[1:4,] ################################################### ### code chunk number 25: Cheminfo.Rnw:1186-1188 ################################################### MW(sdfset[1:4], addH=FALSE) MF(sdfset[1:4], addH=FALSE) ################################################### ### code chunk number 26: Cheminfo.Rnw:1192-1193 ################################################### groups(sdfset[1:4], groups="fctgroup", type="countMA") ################################################### ### code chunk number 27: Cheminfo.Rnw:1203-1209 ################################################### propma <- data.frame(MF=MF(sdfset, addH=FALSE), MW=MW(sdfset, addH=FALSE), Ncharges=sapply(bonds(sdfset, type="charge"), length), atomcountMA(sdfset, addH=FALSE), groups(sdfset, type="countMA"), rings(sdfset, upper=6, type="count", arom=TRUE)) propma[1:4,] ################################################### ### code chunk number 28: Cheminfo.Rnw:1220-1222 (eval = FALSE) ################################################### ## library(ChemmineOB) ## propOB(sdfset)[1:4,] ################################################### ### code chunk number 29: Cheminfo.Rnw:1232-1234 ################################################### datablock(sdfset) <- propma # Works with all SDF components datablock(sdfset)[1] ################################################### ### code chunk number 30: Cheminfo.Rnw:1237-1239 ################################################### blockmatrix <- datablock2ma(datablocklist=datablock(sdfset)) blockmatrix[1:2,1:12] ################################################### ### code chunk number 31: Cheminfo.Rnw:1249-1252 ################################################### bonds(sdfset[[1]], type="bonds")[1:4,] bonds(sdfset[1:2], type="charge") bonds(sdfset[1:2], type="addNH") ################################################### ### code chunk number 32: Cheminfo.Rnw:1262-1263 ################################################### (ringatoms <- rings(sdfset[1], upper=Inf, type="all", arom=TRUE, inner=FALSE)) ################################################### ### code chunk number 33: struct4 ################################################### atomindex <- as.numeric(gsub(".*_", "", unique(unlist(ringatoms$RINGS)))) plot(sdfset[1], print=FALSE, colbonds=atomindex) ################################################### ### code chunk number 34: Cheminfo.Rnw:1288-1298 ################################################### write.SDF(sdfset, "test.sdf") desc <- function(sdfset) { cbind(SDFID=sdfid(sdfset), MW=MW(sdfset), groups(sdfset), rings(sdfset, type="count", upper=6, arom=TRUE) ) } sdfStream(input="test.sdf", output="matrix.xls", fct=desc, Nlines=1000, silent=TRUE) read.delim("matrix.xls", row.names=1)[1:3,1:10] ################################################### ### code chunk number 35: Cheminfo.Rnw:1315-1325 (eval = FALSE) ################################################### ## ## Task 1 ## propma <- data.frame(MF=MF(p450, addH=FALSE), MW=MW(p450, addH=FALSE), ## Ncharges=sapply(bonds(p450, type="charge"), length), ## atomcountMA(p450, addH=FALSE), groups(p450, ## type="countMA"), rings(p450, upper=6, type="count", ## arom=TRUE)) ## ## Task 2 ## datablock(p450) <- propma ## ## Task 3 ## write.SDF(p450, file="p450mod.sdf", sig=TRUE) ################################################### ### code chunk number 36: Cheminfo.Rnw:1337-1339 ################################################### apset <- sdf2ap(sdfset) apset ################################################### ### code chunk number 37: Cheminfo.Rnw:1342-1346 (eval = FALSE) ################################################### ## showClass("APset") ## cid(apset) ## view(apset) ## as(apset, "list") ################################################### ### code chunk number 38: Cheminfo.Rnw:1349-1351 ################################################### (fpset <- desc2fp(apset)) view(fpset[1]) ################################################### ### code chunk number 39: Cheminfo.Rnw:1361-1362 ################################################### cmp.search(apset, apset["CMP1"], type=3, cutoff = 0.3, quiet=TRUE) ################################################### ### code chunk number 40: Cheminfo.Rnw:1365-1368 ################################################### fpset1024 <- names(rev(sort(table(unlist(as(apset, "list")))))[1:1024]) fpset <- desc2fp(apset, descnames=fpset1024, type="FPset") fpSim(fpset["CMP1"], fpset, method="Tanimoto", cutoff=0.2, top=6) ################################################### ### code chunk number 41: Cheminfo.Rnw:1378-1380 ################################################### fpset <- fp2bit(sdfsample, type=3) fpSim(fpset[1], fpset, method="Tanimoto", cutoff=0.2, top=6) ################################################### ### code chunk number 42: Cheminfo.Rnw:1383-1384 ################################################### fpSim(fpset[1], fpset[2]) ################################################### ### code chunk number 43: struct5 ################################################### library(fmcsR); data(fmcstest) # Loads library and test sdfset object test <- fmcs(fmcstest[1], fmcstest[2], au=2, bu=1) # Searches for MCS with mismatches plotMCS(test) # Plots both query compounds with MCS in color ################################################### ### code chunk number 44: Cheminfo.Rnw:1404-1405 ################################################### fmcsBatch(sdfset[[1]], sdfset)[1:2,] ################################################### ### code chunk number 45: Cheminfo.Rnw:1415-1417 (eval = FALSE) ################################################### ## compounds <- searchSim(sdfset[1]) ## compounds ################################################### ### code chunk number 46: Cheminfo.Rnw:1437-1447 (eval = FALSE) ################################################### ## ## Task 1 ## p450 <- read.SDFset("p450.sdf"); cid(p450) <- sdfid(p450) ## apset <- sdf2ap(p450) ## fpset1024 <- names(rev(sort(table(unlist(as(apset, "list")))))[1:1024]) ## fpset <- desc2fp(apset, descnames=fpset1024, type="FPset") ## pcset <- fp2bit(p450, type=3) ## ## Task 2 ## cmp.search(apset, apset[1], type=2, cutoff = 0.3, quiet=TRUE)[1:6] ## fpSim(fpset[1], fpset, method="Tanimoto", cutoff=0.2, top=6) ## fpSim(pcset[1], pcset, method="Tanimoto", cutoff=0.2, top=6) ################################################### ### code chunk number 47: Cheminfo.Rnw:1459-1463 ################################################### c1 <- cmp.cluster(fpset, cutoff=c(0.3, 0.6, 0.9), method="Tanimoto", quiet=TRUE) c1[1:8,] cluster.sizestat(c1, cluster.result=2) ################################################### ### code chunk number 48: Cheminfo.Rnw:1473-1474 ################################################### jarvisPatrick(nearestNeighbors(fpset, numNbrs=6), k=5, mode="a1a2b")[1:20] ################################################### ### code chunk number 49: Cheminfo.Rnw:1477-1479 ################################################### nnm <- nearestNeighbors(fpset,numNbrs=6) nnm$similarities[1:4,] ################################################### ### code chunk number 50: struct6 ################################################### simMA <- sapply(cid(fpset), function(x) fpSim(fpset[x], fpset, sorted=FALSE)) plot(cmdscale(as.dist(1-simMA))) ################################################### ### code chunk number 51: struct7 ################################################### library(gplots) hc <- hclust(as.dist(1-simMA), method="single") heatmap.2(1-simMA, Rowv=as.dendrogram(hc), Colv=as.dendrogram(hc), col=colorpanel(40, "darkblue", "yellow", "white"), density.info="none", trace="none") ################################################### ### code chunk number 52: Cheminfo.Rnw:1531-1545 (eval = FALSE) ################################################### ## ## Task 1 ## p450 <- read.SDFset("p450.sdf"); cid(p450) <- sdfid(p450) ## apset <- sdf2ap(p450) ## fpset1024 <- names(rev(sort(table(unlist(as(apset, "list")))))[1:1024]) ## fpset <- desc2fp(apset, descnames=fpset1024, type="FPset") ## cmp.cluster(fpset, cutoff=c(0.8, 0.85, 0.9), method="Tanimoto", quiet=TRUE) ## ## Task 2 ## jarvisPatrick(nearestNeighbors(fpset, numNbrs=6), k=5, mode="a1a2b") ## ## Task 3 ## simMA <- sapply(cid(fpset), function(x) fpSim(fpset[x], fpset, sorted=FALSE)) ## plot(cmdscale(as.dist(1-simMA))) ## ## Task 4 ## hc <- hclust(as.dist(1-simMA), method="single") ## heatmap.2(1-simMA, Rowv=as.dendrogram(hc), Colv=as.dendrogram(hc), col=colorpanel(40, "darkblue", "yellow", "white"), density.info="none", trace="none") ################################################### ### code chunk number 53: Cheminfo.Rnw:1553-1554 ################################################### sessionInfo()