Bio3D in R

 

Contribute

Page history last edited by Barry Grant 1 yr ago

Contributing to Bio3d

We are always interested in adding additional functionality to the bio3d library. If you have ideas, sugestions or code that you would like to distribute as part of this package, please contact us.

  

Writing your own functions 

One of the main attractions of using the R environment is the ease with which users can write their own programs and functions. The R programming syntax is extremely easy to learn, even for users with no previous programming experience. Once the basic R programming control structures are understood, users can use the language as a powerful calculation machine to perform complex analyses of almost any type of quantitative data.

 

As an example, below we construct a simple function to extract the 1-letter aminoacid sequence from a specified PDB file. Note the addational input argument 'seqfile' which if set to TRUE will cause a FASTA format sequence file to be writen to the current directory.

 

pdb2fasta <- function(file, seqfile=FALSE) {

  pdb <- read.pdb(file)

  seq.pdb <- aa321(pdb$atom[pdb$calpha,"resid"])

 

  if(seqfile) {

    name <- basename(file)

    write.fasta(seqs=seq.pdb, ids=name,

                file=paste(name,".fasta",sep=""))

  }

  return(seq.pdb)

}

 

 

If you are using R in the standard interactive fashion, then pasting the above function into your R sesion and then issuing the comand should print the sequence and produce a file called '1w62.pdb.fasta'.

 

pdb2fasta("http://www.rcsb.org/pdb/files/1bg2.pdb", TRUE)

 

 

Using your code as a script

If you would rather use your code as an exicuitable script then one option is to use 'Rscript' which should reside in the same directory as your R exicuitable. Below is a simple example script.

 

#!/usr/bin/env Rscript

library(bio3d)

got <- getArgs( c("pdb","seqfile") )

 

##-- Check cmd options

if(!got["pdb"])

  stop("Command line argument 'pdb=?' is required")

 

if(!got["seqfile"]) {

  warning("Assuming 'seqfile=TRUE': producing FASTA file")

  seqfile=TRUE

}

 

pdb2fasta <- function(file, seqfile=FALSE) {

  pdb <- read.pdb(file)

  seq.pdb <- aa321(pdb$atom[pdb$calpha,"resid"])

 

  if(seqfile) {

    name <- basename(file)

    write.fasta(seqs=seq.pdb, ids=name,

                file=paste(name,".fasta",sep=""))

  }

  return(seq.pdb)

}

 

pdb2fasta(pdb, seqfile)

 

 

 

Copying the above code to a file named pdb2fasta.r and issuing the following UNIX comands.

chmod +x pdb2fasta.r

./pdb2fasta.r pdb=http://www.rcsb.org/pdb/files/1bg2.pdb

 

 

Comments (0)

You don't have permission to comment on this page.