Added code to reproduce the figure in the paper

Author: Digp

FigurePaper/example1/dependencies.R

@@ -0,0 +1,59 @@
+## "cache" saves trinucleotide backbones in RAM for latter RMSD calculations
+cache <-
+function(ntID, name, ntinfo) {
+    ## Trim trinucleotide backbone
+    pdb <- trimByID(ntID=ntID, ntinfo=ntinfo, prev=1, post=1,
+                    alt=c("uniq"), cutoff=0, justbackbone=TRUE,
+                    verbose=FALSE)
+    ## Save as an independent pdb file in RAM
+    assign(name, pdb, envir=.GlobalEnv)
+}
+
+## "pwrmsd" calculates pair wise RMSD between nucleotides presaved in RAM
+pwrmsd <-
+function(ntIDs) {
+    pdb1 <- get(paste("nt", ntIDs[1], sep=""))
+    pdb2 <- get(paste("nt", ntIDs[2], sep=""))
+    sel1 <- atom.select(pdb1)
+    sel2 <- atom.select(pdb2)
+    fit=fit.xyz(pdb1$xyz, pdb2$xyz,
+                fixed.inds=sel1$xyz, mobile.inds=sel2$xyz)
+    return(rmsd(pdb1$xyz[,sel1$xyz],fit[,sel2$xyz]))
+}
+
+## "filter_pyle" reproduces the iterative process of RMSD comparison to chose
+## an helical nucleotide of reference
+filter_pyle <- 
+function(ntinfo, helicalntID, RMSD, bandwidths=c(40,40), cutoff=0.85) {
+
+    twideMatrix <- t(acast(RMSD, nt1~nt2, value.var="rmsd"))
+    twideMatrix<-rbind(rep(NA,length(helicalntID)-1),twideMatrix)
+    twideMatrix<-cbind(twideMatrix,rep(NA,length(helicalntID)))
+    rownames(twideMatrix)[1]<-helicalntID[1]
+    colnames(twideMatrix)[length(helicalntID)]<-helicalntID[length(helicalntID)]
+
+    Matrix <- twideMatrix
+    while (length(helicalntID) > 1) {
+        Matrix <- Matrix[which(rownames(Matrix) %in% helicalntID),
+                            which(colnames(Matrix) %in% helicalntID)]
+
+        column_sums <- colSums(Matrix, na.rm=TRUE)
+        row_sums <- rowSums(Matrix, na.rm=TRUE)
+        RMSDaverages <- (column_sums + row_sums)/(length(row_sums) - 1)
+        index <- which.min(RMSDaverages)
+        lowestaverage <- RMSDaverages[index]
+        reference_nt <- helicalntID[index]
+        A <- Matrix[which(rownames(Matrix) == reference_nt), ]
+        names(A) <- colnames(Matrix)
+        A <- A[complete.cases(A)]
+        B <- Matrix[, which(colnames(Matrix) == reference_nt)]
+        names(B) <- colnames(Matrix)
+        B <- B[complete.cases(B)]
+        reference_nt_RMSD <- append(A,B); rm(list=c("A","B"))
+        helicalntID <- sort(append(reference_nt, as.numeric(
+                            names(reference_nt_RMSD)[which(
+                            reference_nt_RMSD < lowestaverage)])))
+    }
+
+    return(helicalntID)
+}

FigurePaper/example1/etatheta.R

@@ -0,0 +1,72 @@
+#!/usr/bin/Rscript
+## EXAMPLE 1
+## Description: From Leontis non-redundant RNA list, get eta-theta plot
+##              following Pyle's method.
+## Author: Diego Gallego
+
+## ----------------------------------------------------------------------------
+## Preprocess
+## ----------------------------------------------------------------------------
+## Load necessary packages for the example
+library(veriNA3d)
+library(bio3d)
+library(reshape2) ## Specific dependency for this example!
+
+## Load specific functions for the example
+source("dependencies.R")
+
+## ----------------------------------------------------------------------------
+## Get Raw data
+cat("Getting data\n")
+## ----------------------------------------------------------------------------
+
+## Get latest Leontist non-redundant list of RNA structures
+infochains <- getLeontisList(threshold="2.5A", as.df=TRUE)
+
+## Get structural data
+ntinfo <- pipeNucData(pdbID=infochains$pdb, model=infochains$model,
+                        chain=infochains$chain, cores=1)
+
+## ----------------------------------------------------------------------------
+## Clean data. Raw -> Tidy
+cat("Cleaning data\n")
+## ----------------------------------------------------------------------------
+
+## Subset north nucleotides
+north <- ntinfo[cleanByPucker(ntinfo, surenorth=TRUE), ]
+
+## Find non-broken canonical nucleotides
+id <- which(north$puc_valid==TRUE & north$Break==FALSE & north$big_b == FALSE & 
+            north$eta_valid == TRUE & north$theta_valid == TRUE &
+            grepl("^[GAUC]-[GAUC]-[GAUC]$", north$localenv, perl=TRUE))
+usefulnt <- north$ntID[id]
+
+## Filter helical nucleotides before ploting =================================
+
+## Create all the necessary trinucleotides for pair-wise RMSD calculations
+invisible(mapply(FUN=cache, ntID=usefulnt, name=paste("nt", usefulnt, sep=""),
+                    MoreArgs=list(ntinfo=ntinfo)))
+
+## Find helical nucleotides
+HDR <- findHDR(ntID=usefulnt, ntinfo=ntinfo, x="eta", y="theta", SD_DENS=15)
+
+## Compute pair-wise RMSD between helical nucleotides
+pairwise <- t(combn(HDR[[1]], 2))
+RMSD_calc <- apply(FUN=pwrmsd, MARGIN=1, X=pairwise)
+RMSD <- data.frame(nt1=pairwise[, 1], nt2=pairwise[, 2], rmsd=RMSD_calc)
+
+## Apply Pyle method
+ref <- filter_pyle(ntinfo=ntinfo, helicalntID=HDR[[1]], RMSD=RMSD)
+
+pairs <- data.frame(usefulnt=usefulnt, reference=rep(ref, length(usefulnt)))
+RMSD2 <- apply(FUN=pwrmsd, MARGIN=1, X=pairs)
+
+## ===========================================================================
+## ----------------------------------------------------------------------------
+## Exploratory analysis. Plot
+cat("Plotting data\n")
+## ----------------------------------------------------------------------------
+tiff("example1.tiff", height=10.5, width=16.5, units="cm", res=300)
+plot2D(ntinfo=ntinfo, x="eta", y="theta", etatheta=TRUE,
+        ntID=pairs[which(RMSD2 > 0.85), 1])
+dev.off()