Commit a9bf129b authored by Digp's avatar Digp
Browse files

Deleted old folders

parent 979b28db
Pipeline #3091 skipped
## "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)
}
#!/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()
This source diff could not be displayed because it is too large. You can view the blob instead.
real 125m49.951s
user 113m31.407s
sys 0m21.852s
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment