Uniprot → SVG

My colleague @Solena recently asked me to write a tool that would help her to prepare a large number of figures for an article. The tool I wrote fetches an entry for a given Uniprot accession and creates a SVG diagram , editable in Inkscape. It is available at

http://lindenb.github.io/pages/uniprot/paintsvg.html

That's it,

Pierre

Translating a DNA sequence in Google Spreadsheet using #GoogleAppScript

Google has released Google App Script.
"Google Apps Script is a JavaScript cloud scripting language that lets you extend Google Apps and build web applications. Scripts are developed in Google Apps Script’s browser-based script editor, and they are stored in and run from Google's servers.
Google Apps Script is very versatile. Here are some examples of things you can do with Google Apps Script":

• Build custom functions in a Google Spreadsheet
• Extend certain Google Apps products by creating custom menus linked to scripts
• Create and publish web applications, which can run on their own or embedded within a Google Site
• Schedule tasks like report creation and distribution and run them on a custom schedule
• Automate workflows such as document or expense approvals, order fulfillment, time-tracking, and more

In the current post, I'll show how to create a custom javascript function that will

Translate a
DNA
to a
Protein
into a
Google Spreadsheet

Create a new Google Spreadhseet. Open the menu "Tools" > "Script Manager...". Click on New...

Click on "Create Blank Project".

A new editor is opened. Copy the following javascript code (https://gist.github.com/3716137) into the editor. Save the javascript projet.

Close the script, go back to the spreasheet. You can now use your new function =translateDNA(dna):

That's it,

Pierre

Drawing a protein (Biostar #6172)

This post is my answer for this question on Biostar:Drawing a protein:
Dear all I often find protein's image like this (...) Do you know if there's a program to draw them (I mean circles with letters).

I wrote a Java-Swing application named WirePeptide displaying a draggable peptide. This application is available on github at https://github.com/lindenb/jsandbox/blob/master/src/sandbox/WirePeptide.java. The user can save the image as PNG, SVG and HTML+Canvas.

Compile & Run

cd jsandbox
ant wirepeptide
java -jar dist/wirepeptide.jar

Result (Canvas)

That's it,

Pierre

Translating a DNA to a Protein using server-side javascript and C: my notebook

In my previous post , I used Node.js to translate a DNA to a protein on the Server-side, using javascript. In the following post, I again will translate a DNAn but this time by calling a specialized C program on the server side.

Source code

The C program

The C program reads a DNA string from stdin a translate it using the standard genetic code:
Compilation:

gcc -o /my/bin/path/translate translate.c

The Node.js script

When the Node.js server receive a DNA parameter, it spawns a new process to the C program and we write the DNA to this process via 'stdin'.
Each time a new 'data' event (containing the protein) is received, it is printed to the http response. At the end of the process, we close the stream by calling 'end()'.

test

> node-v0.2.5/node translate.js
Server running at http://127.0.0.1:8080

> curl -s "http://localhost:8080/?dna=ATGATGATAGATAGATATAGTAGATATGATCGTCAGCCATACG"
MMIDRYSRYDRQPY

That's it,

Pierre

Server-side javascript: translating a DNA with Node.js

(wikipedia) Node.js is an evented I/O framework for the V8 JavaScript engine on Unix-like platforms. It is intended for writing scalable (javascript-based) network programs such as web servers.

In the following post I will create a javascript server translating a DNA to a protein.

Installing Node.js

I've downloaded the sources for Node.js from http://nodejs.org/#download. It compiled (configure+make) and ran without any problem.

The script

The following script contains a class handling a GeneticCode and the server TranslateDna translating the DNA to a protein, it handles both the POST and the GET http methods. It no parameter is found it displays a simple HTML form, else the form data are decoded and the DNA is translated. The protein is returned as a JSON structure.

Running the server

> node-v0.2.5/node translate.js
Server running at http://127.0.0.1:8080

Test

> curl "http://localhost:8080/"
<html><body><form action="/" method="GET"><h1>DNA</h1><textarea name="dna"></textarea><br/><input type="submit" value="Submit"></form></body></html>

> curl "http://localhost:8080/?dna=ATGAACTATCGATGCTACGACTGATCG"
{"protein":"MNYRCYD*S","query":"ATGAACTATCGATGCTACGACTGATCG"}

That's it,

Pierre

Where are the alternative reading frames in the Human Genome ?

The following post was inspired by a question asked recently on Biostar :"Do exons ever have different reading frames in spliced variants?".

To find those alternative reading frames I've used the table KnownGene available at UCSC from : http://hgdownload.cse.ucsc.edu/goldenPath/hg18/database/knownGene.txt.gz. This file contains the positions of the exons for each transcript in the human genome:

mysql -h genome-mysql.cse.ucsc.edu -A -u genome -D hg18 -e 'select * from knownGene limit 10\G'

(...)
*************************** 7. row ***************************
name: uc009vis.1
chrom: chr1
strand: -
txStart: 4268
txEnd: 6628
cdsStart: 4268
cdsEnd: 4268
exonCount: 4
exonStarts: 4268,4832,5658,6469,
exonEnds: 4692,4901,5805,6628,
proteinID:
alignID: uc009vis.1
*************************** 8. row ***************************
name: uc009vit.1
chrom: chr1
strand: -
txStart: 4268
txEnd: 9622
cdsStart: 4268
cdsEnd: 4268
exonCount: 9
exonStarts: 4268,4832,5658,6469,6720,7095,7777,8130,8775,
exonEnds: 4692,4901,5810,6628,6918,7605,7924,8229,9622,
proteinID:
alignID: uc009vit.1

The following java program creates an array of bytes having a length greater than the length of the human chromosome chr1. This array is initialized with the constant 'NIL'. Then for each chromosome and each transcript, we loop over each exon and we record what was the reading frame (0, 1 or 2) at a given position. If this position was already flagged with another frame, a warning is printed to stdout.

Compilation & Execution

javac BioStar3034.java
java BioStar3034

Result

(...)
chr1:53286155-53286156 (+)
chr1:53286156-53286157 (+)
chr1:53286157-53286158 (+)
chr1:53286158-53286159 (+)
chr1:53286159-53286160 (+)
chr1:53286160-53286161 (+)
chr1:53286161-53286162 (+)
chr1:53286162-53286163 (+)
(...)
java BioStar3034 | sort | uniq | wc -l
300696

(Image from UCSC/OpenWetWare)

That's it
Pierre

'R' = dna.translate("AGG") . A custom C function for R, My notebook.

In the following post, I will show how I've implemented a custom C function for R. This C function will translate a DNA to a protein. I'm very new to 'R' so feel free to make any comment about the code.

C code

The data in 'R' are stored in an opaque structure named 'SEXP'. A custom C function receives some SEXP arguments and returns another SEXP. So, the declaration of my function translate_dna is:

SEXP translate_dna(SEXP sequences)
{
(...)
}

We check if the argument 'sequences' has a type='character'.

if(!isString(sequences))
error("argument is not a string[]");

The input consists of one or more DNA, so our return value will be an SEXP array of type 'character' (STRSXP) with a size equals to the number of DNAs. This array is allocated:

SEXP array_of_peptides = allocVector(STRSXP,length(sequences));

We loop over each DNA sequence

for(in i=0;i<length(sequences);++i) ...

and we allocate some memory for the peptide:

const char *dna=CHAR(STRING_ELT(sequences, i));
int dna_length=strlen(dna);
char* peptide =malloc(dna_length/3+1);

The peptide is then filled with its amino acids:

for(j=0;j+2 < dna_length;j+=3)
{
peptide[n++]=_translate(dna[j],dna[j+1],dna[j+2]);
}
peptide[n]=0;

And we put this new string/peptide in the returned value

SET_STRING_ELT(array_of_peptides,i,Rf_mkChar(peptide));

Here, I am not sure how 'R' handles its memory. Should I care about the way R runs its garbage manager ? how should I use the macros PROTECT and UNPROTECT ?

Compilation

The code 'translate.c' is compiled as a dynamic library ' libtranslate.so':

gcc -fPIC -I -g -c -Wall -I ${R_HOME}/include translate.c
gcc -shared -Wl,-soname,libtranslate.so.1 -o libtranslate.so translate.o

R code

on the R side , the previous dynamic library is loaded:

dyn.load(paste("libtranslate", .Platform$dynlib.ext, sep=""))

A function dna.translate is declared: it forces the array to be an array of string and it invokes the C function 'translate_dna'

storage.mode(dna) <- "character"
.Call("translate_dna", dna)

We can now invoke the R function dna.translate:

peptides <- dna.translate(
c( "ATGGAGAGGCAGAAACGGAAGGCGGACATCGAGAAAGGGCTGCAGTTCATTCAGTCGACACTAC",
NULL,
"CCCAAAAGCAAGAAGAATATGAGGCCTTTCTGCT",
"CAAACTGGTGCAGAATCTGTTTGCTGAGGGCAATGA",
NULL,
"GGCAGATCAGGGAACTTCTAATGGATTGGGGTCC",
"GGATAACTGCACCTTCGCCTACCATCAGGAGGA",
"GGGTCCCAGGCAGCGCTGCCTGGGGGCTGGGGAG",
"TTGCGACAGGCTCCAGAAGGGCAAAGCCTGCCCAGAT",
"GCACCCCCCTCCTCTCCACCCTACCTTCCATCAACCA",
"AGG"
))

print(peptides)

Result:

[1] "MERQKRKADIEKGLQFIQSTL" "PKSKKNMRPFC" "QTGAESVC*GQ*"
[4] "GRSGNF*WIGV" "G*LHLRLPSGG" "GSQAALPGGWG"
[7] "LRQAPEGQSLPR" "APPSSPPYLPST" "R"

Full source code

:
C code:

#include <stdio.h>
#include <ctype.h>
#include <R.h>
#include <Rinternals.h>
/* the genetic code */
static const char* STANDARD_GENETIC_CODE="FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG";

static char _translate(char a,char b,char c);


/** translates DNA to protein
* @arg sequences one or more DNA sequence
* @return an array of peptides
*/
SEXP translate_dna(SEXP sequences)
{
int i;
SEXP array_of_peptides;
//check input has type= characters
if(!isString(sequences))
error("argument is not a string[]");
//prepare a new list for the translated sequence
array_of_peptides = allocVector(STRSXP,length(sequences));

//loop over the input
for(i=0;i< length(sequences);++i)
{
int j=0,n=0;
//transform the sequence into a ptr* of char
const char *dna=CHAR(STRING_ELT(sequences, i));
//get the length of this sequence
int dna_length=strlen(dna);
//alloc the protein sequence
char* peptide =malloc(dna_length/3+1);
if(peptide==NULL) error("out of memory");
//loop over the codons
for(j=0;j+2 < dna_length;j+=3)
{
//set the amino acid at 'n'
peptide[n++]=_translate(dna[j],dna[j+1],dna[j+2]);
}
//add EOS
peptide[n]=0;
//put a copy of this peptide in the vector
SET_STRING_ELT(array_of_peptides,i,Rf_mkChar(peptide));
//free our copy
free(peptide);
}
//return the array of petptide
return array_of_peptides;
}


static int base2index(char c)
{
switch(tolower(c))
{
case 't': return 0;
case 'c': return 1;
case 'a': return 2;
case 'g': return 3;
default: return -1;
}
}

static char _translate(char a,char b,char c)
{
int base1= base2index(a);
int base2= base2index(b);
int base3= base2index(c);
if(base1==-1 || base2==-1 || base3==-1)
{
return '?';
}
else
{
return STANDARD_GENETIC_CODE[base1*16+base2*4+base3];
}

}

Makefile:

run:translate.so
${R_HOME}/bin/R --no-save < translate.R

translate.so:translate.c
gcc -fPIC -I -g -c -Wall -I ${R_HOME}/include translate.c
gcc -shared -Wl,-soname,libtranslate.so.1 -o libtranslate.so translate.o

R code:

#load the dynamic library translate
dyn.load(paste("libtranslate", .Platform$dynlib.ext, sep=""))
#declare the function dna.translate
dna.translate <- function(dna)
{
#force dna to be type=character
storage.mode(dna) <- "character"
#call the C function
.Call("translate_dna", dna)
}

peptides <- dna.translate(
c( "ATGGAGAGGCAGAAACGGAAGGCGGACATCGAGAAAGGGCTGCAGTTCATTCAGTCGACACTAC",
NULL,
"CCCAAAAGCAAGAAGAATATGAGGCCTTTCTGCT",
"CAAACTGGTGCAGAATCTGTTTGCTGAGGGCAATGA",
NULL,
"GGCAGATCAGGGAACTTCTAATGGATTGGGGTCC",
"GGATAACTGCACCTTCGCCTACCATCAGGAGGA",
"GGGTCCCAGGCAGCGCTGCCTGGGGGCTGGGGAG",
"TTGCGACAGGCTCCAGAAGGGCAAAGCCTGCCCAGAT",
"GCACCCCCCTCCTCTCCACCCTACCTTCCATCAACCA",
"AGG"
))

print(peptides)

That's it

Pierre