β.Log [ Computational Biology News ]

Shotgun SILAC-MS data Analysis with just MaxQuant and Excel

2021-03-30T14:59:00.004+05:30

Novel insert in Spike glycoprotein of Severe acute respiratory syndrome coronavirus 2

2020-03-07T22:20:00.000+05:30

BLAST result alignment view for Spike glycoprotein of Coronavirus (sp|P0DTC2)

Fig. 1. BLAST result alignment view for Spike glycoprotein of Coronavirus (sp|P0DTC2)

One of the sequences in the pre-release proteome of 2019-nCoV Wuhan Coronavirus [1] is Spike glycoprotein (sp|P0DTC2). This sequence, when aligned with BLAST against nucleotide database NCBI [2], shows a peptide insert PRRA (Figure 1) with respect to Bat coronavirus RaTG13 [2]. Further, the last couple of residues in this insert appear as IL and VL in Bat SARS-like coronavirus [3] and Recombinant coronavirus clone [4] respectively. However, the SARS coronavirus https://www.ncbi.nlm.nih.gov/nucleotide/AY395003.1 has either LL ending or completely missing this peptide as in https://www.ncbi.nlm.nih.gov/nucleotide/DQ640652.1 GDH-BJH01. The furin-cleavage site of this protein seems to be crucial for Angiotensin-converting enzyme 2 receptor interaction and further entry into the host [5]. Further the PR as dipeptide is shown to be crucial for structure and function [6]. So this site might be something to look into in order to detect and possibly control infection of coronavirus. The current structure for the protein http://www.rcsb.org/structure/6VW1 doesn’t cover this location. However the proteome [7] released ftp.pride.ebi.ac.uk/pride/data/archive/2020/03/PXD017710/ shows the increasing amounts of this protein P0DTC2 at various timepoints (http://corona.papers.biochem2.com/) but there seems to be a batch effect (Figure 3) possibly due to the multiplexing (https://biochem2.com/images/IBCII-Images/Research_Groups_Images/Projects/PQC/TMT_scheme.xlsx) . This seems to be mostly captured in first two principal components (Figure 4).

Fig 2. Rising levels of P0DTC2 with timepoints 2,6,10 and,24 hours http://corona.papers.biochem2.com/

Fig 3. SARS-CoV-2 infected host cell proteomics reveal potential batch effect, hierarchical clustering based on the correlation-coefficient calculated over quantification results [7]

Fig 4. PCA over protein groups quantified across the samples in Figure 2

References

[1] 2019-nCoV Wuhan Coronavirus protein sequences, (n.d.). ftp://ftp.uniprot.org/pub/databases/uniprot/pre_release (accessed March 7, 2020).

[2] Bat coronavirus RaTG13, complete genome - Nucleotide - NCBI, (n.d.). https://www.ncbi.nlm.nih.gov/nucleotide/MN996532.1?report=genbank&log$=nuclalign&blast_rank=53&RID=675M66YC016 (accessed March 7, 2020).

[3] D. Hu, C. Zhu, L. Ai, T. He, Y. Wang, F. Ye, L. Yang, C. Ding, X. Zhu, R. Lv, J. Zhu, B. Hassan, Y. Feng, W. Tan, C. Wang, Genomic characterization and infectivity of a novel SARS-like coronavirus in Chinese bats, Emerg. Microbes Infect. 7 (2018) 154. https://doi.org/10.1038/s41426-018-0155-5.

[4] M.M. Becker, R.L. Graham, E.F. Donaldson, B. Rockx, A.C. Sims, T. Sheahan, R.J. Pickles, D. Corti, R.E. Johnston, R.S. Baric, M.R. Denison, Synthetic recombinant bat SARS-like coronavirus is infectious in cultured cells and in mice, Proc. Natl. Acad. Sci. U. S. A. 105 (2008) 19944–19949. https://doi.org/10.1073/pnas.0808116105.

[5] A.C. Walls, Y.-J. Park, M.A. Tortorici, A. Wall, A.T. McGuire, D. Veesler, Structure, function and antigenicity of the SARS-CoV-2 spike glycoprotein, BioRxiv. (2020) 2020.02.19.956581. https://doi.org/10.1101/2020.02.19.956581.

[6] L.R. Hayes, L. Duan, K. Bowen, P. Kalab, J.D. Rothstein, C9orf72 arginine-rich dipeptide repeat proteins disrupt karyopherin-mediated nuclear import, ELife. 9 (2020) e51685. https://doi.org/10.7554/eLife.51685.

[7] SARS-CoV-2 infected host cell proteomics reveal potential therapy targets | Research Square, (n.d.). https://www.researchsquare.com/article/rs-17218/v1 (accessed March 14, 2020).

PRRA peptide search https://www.uniprot.org/peptidesearch/ against Uniprot

Viruses (82 results)

Alphaherpesvirinae (18 results)

Human herpesvirus 1 (HHV-1) (Human herpes simplex virus 1) (3 results)

Human herpesvirus 2 (strain HG52) (HHV-2) (Human herpes simplex virus 2) (4 results)

Varicellovirus (8 results)

Bovine alphaherpesvirus 1 (4 results)

Equid alphaherpesvirus 1 (Equine herpesvirus 1) (2 results)

Suid alphaherpesvirus 1 (2 results)

Betaherpesvirinae (7 results)

Human cytomegalovirus (HHV-5) (Human herpesvirus 5) (4 results)

Roseolovirus (3 results)

Epstein-Barr virus (strain GD1) (HHV-4) (Human herpesvirus 4) (3 results)

Alphavirus (arboviruses group A) (3 results)

Sleeping disease virus (SDV) (1 results)

Avian infectious bursal disease virus (IBDV) (Gumboro disease virus) (10 results)

Drosophila x virus (isolate Chung/1996) (DXV) (1 results)

Hepatitis C virus genotype 6 (2 results)

Langat virus (2 results)

Luteovirus (5 results)

Barley yellow dwarf virus (isolate MAV) (BYDV) (2 results)

Barley yellow dwarf virus (isolate MAV-PS1) (BYDV) (1 results)

Barley yellow dwarf virus (isolate PAV) (BYDV) (3 results)

Barley yellow dwarf virus (isolate P-PAV) (BYDV) (1 results)

Nodaviridae (4 results)

Betanodavirus (3 results)

Boolarra virus (BoV) (1 results)

Sapporo virus (2 results)

Tymovirales (2 results)

cellular organisms (2261 results)
Keywords: #Spike #glycoprotein #SARS-CoV-2 #COVID-19 #SARS-CoV-2 #HCoV-19 #SARS2, #2019-nCoV, #Wuhan #Coronavirus #SeafoodMarket #Human #pneumonia #virus

Creating real-time audio-visual effects using brain waves

2014-02-03T05:05:00.001+05:30

Finally got some quality time to spend with a birthday gift from Katharina. It is one of those dry EEG sensor containing devices which can capture the EEG signal just by placing the sensor on the forehead, no gel, no hell ;) To top it all, NeuroSky mobile, can connect to any Bluetooth device and through a good collection of APIs, EEG signals and be collected and processed. This felt like the easiest way to enter the amazing world of brain-computer interface :) and so I did...

Quick Google search showed that Windows+Python would be the fastest way to hack together a piece of code which can unleash some of the power of this device i have in my hands, or rather on my head! So off i went and downloaded python, pyserial (for bluetooth stream), numpy (nice array handling and some Fourier magic), matplotlib (plotting, though some dependencies to make this work are dateutil, pyparsing and six) along with NeuroPy, Sahil's python library to programmatically make NeuroSky function calls.

Once the setup was nice and working, i fired up the IDLE, python IDE, which comes with default install and stated coding. After rummaging through documentation of the above modules and extremely helpful suggestions available at stack overflow, i am finally able to capture and (dis)play frequencies in my EEG signal. The code is not so readable at the moment, apologies for that, but it is out there in my repo, if you feel like taking a look:

wget http://misccb.googlecode.com/git-history/519b92947fb874d3d3842e08518fbb42ed93adfa/neuroplay.py

Easiest way to run this code is via IDLE, otherwise just fire up the Windows command prompt <%windir%\system32\cmd.exe>:
cd "path to python install dir"
python "path to code download dir"\neuroplay.py

I will attempt to replace windsound with Olivier's pyo module and run this over Linux. So much looking forward to the next block of quality time to spend on random things :)

Batch script for running MaxQuant command line tool

2013-01-16T22:00:00.001+05:30

I had to process multiple Orbitrap raw files using the same parameter with Max Quant, was not able to find a simple tool for this, thus I wrote a batch script which can be downloaded from https://docs.google.com/file/d/0BxbjZeVL8S4EQW1zbVd4TzRYSDg/edit .

It needs a preconfigured parameter file to begin with. For that, i generally open the Max Quant for a dummy file called TestFile.raw and set the desired parameters and fasta file for search through its GUI. Then save this as testpar.xml in the data directory using the File->Save Parameters option of MaxQuant GUI.

The script needs to be modified for following variables (can open the mqrun.bat with notepad and modify accordingly):
set MAXQUANTDIR= full path of MaxQuant installation directory
set THREADS= number of threads to run"
set DATADIR= full path of the data directory
set PREFIXRAW= leave it blank if all the raw files need to be processes or use prefix of desired raw files
set PARAMFILE= name of the generated parameter file, in this case it is testpar.xml
set SEARCHTEXT= name of the raw file used for generating the parameter file, in this case it is TestFile

Save the script and run it via double clicking on it through windows explorer, or open a cmd terminal and cd to the directory containing it and type mqrun (or whatever name you used to save the file, but make sure it had .bat extension). It creates a text file in the data directory called tempfile.txt containing the name of raw files. The files will be sequentially processed using MaxQuantCmd.exe and the results, namely combined\txt\* + combined\*.apl files, will be stored in RES folder.

I know it is kind of quick and dirty. Feel free to use it and let me know if you find some bugs :) If you modify the script, kindly update its git version at http://code.google.com/p/misccb/source/browse/mqrun.bat?spec=svnd6632cc2bdc3a102fec0883bdbf0a4691bfe7c8c&r=d6632cc2bdc3a102fec0883bdbf0a4691bfe7c8c (will need to login with Google account).

Why mistakes are often repeated?

2012-05-08T01:00:00.002+05:30

“Insanity is repeating the same mistakes and expecting different results.”

- Narcotics Anonymous

I have been wondering about this for a long long time... I feel our daily lives are often a mixture of random events and results of decisions we make as we go on living. Brain as an advanced pattern recognition device is constantly trying to make sense of this, attributing events to nature and/or our decision(s). Such phenomena can lead to perplexing behaviour, for example pigeons developing a causal relation between their actions and random awards. Slowly such perceived causal connections becomes well ingrained that despite the occurrences of counter events, brain fails to update these rules.

Probably at certain level of abstractness, these rules are presented as heuristics which enables the brain to make decisions in a fast and frugal way. Since diverse probability distributions converge to normal distribution, these heuristics can work like shortcuts, exploiting our intuitive sense of normal distribution and being more right than wrong on average. However as we move away from scale of our sensory perception, these intuitive sense starts failing us. For example we often hear that people find it hard to understand physics, especially mechanics of things at atomic and cosmological scales despite statistical techniques used being similar across domains.

When we reflect on our lives, we often see that on several occasions we failed to learn and keep repeating similar mistakes. Why it is so hard to learn from mistakes? Probably it is due to earlier mentioned heuristic-based decision making ingrained in our brain which may be very hard to replace since they generally work well enough. These are simple shortcuts exploiting our intuitive sense of normal distribution. But when things are distributed differently, say power-law, our estimations fails us, for example things like in top 1%, max and min are separated by a much larger factor than one 'normally' imagines.

In a Bayesian sense, it seems that priors don't really update. For example often tests result for relatively rare conditions are misinterpreted by patients as well as others. Conscious awareness of such mistake is often not enough and training in statistical inference seems essential in minimizing "risk illiterate". This might reflect how ingrained these decision making processes really are.

At the basic computational level, for example considering brain as a fully connected neural network (ANN), we can contemplate how this can happen. ANN like systems often show the non-linear relation between training-time and errors. Finding optimal solutions even in such a simplistic view of brain is hard as the global minimum is obfuscated by multiple local minima in the solution space and back-propagation has a tendency to converge in a local minimum. Since these local minima typically have different characters, probably the system favors the one which leads to less lethal decision selected by millions of years of evolutionary process. Retraining or (un)learning in this context makes sense only if the repeated mistake is lethal enough and cost of the retraining is bearable in the long run. Further the nature of such event should be of the updating existing connections rather than disruptive, thus accommodating already learned decisions and not deleting them.

Exploring the above ideas concretely to observe and explain the phenomenon of repeated-mistakes and our subtle tendencies for knee-jerk reactions rather than critically thought out decisions is what I am going to explore in coming posts. Since this is a computational blog, I plan to use ACT-R, SPM, pyNEST and Weka for this journey :)

Calculating expected distance between consecutive occurrence of a motif

2012-04-02T18:31:00.002+05:30

Recently Mohan asked me an interesting question about positions where "GATC" motif occurs in an E.coli genome(s). While figuring this out, I observed that mean of the distance between consecutive GATC is around ~250 for these genomes (242.6 for the linked one). This made us wonder if this is special or just the expected value? With a lot of help from Pragun and loads of generalization, we found it to be around ~255. Like all(!) probability problems, we might have found the right answer with the wrong approach, thus sharing with all to spot and correct the blunder(s :)

To calculate the theoretical expected distance (d) between consecutive occurrence of a motif consisting of individual nucleotide bases (A,T,G,C), lets assume uniform distribution of individual occurrence of these bases. With this assumption one can say that probability of occurrence of individual bases, P(A)=P(G)=P(T)=P(C)=1/n, where n is number of bases here being 4. So probability of a particular motif in general is 1/4 multiplied by the length of the motif (M). In this particular case it is GATC which is of length 4 making P(M)=1/4^4=1/256=1/N. Further lets assume that the genome is infinitely extended (on both sides to make calculations easier ;)

Abstracting the genome while focusing on a region with consecutive occurrence will look like xxxx(M)...d...(M)xxx.... and expectation of this probability (E) of this (?) is d*P(not finding M before)*P(M)*P(not finding M within distance d)*P(M) where:
P(not finding M before)=Summation_x=0..Inf[(1-1/N)^x]=N
(summing geometric sequence)
d*P(not finding M within distance d)=Summation_d=0..Inf[d*(1-1/N)^d]=(N-1)/N^2
(summing arithmetico geometric sequence)
making E=N-1=255 for GATC :)

To check out the results, just wrote a script to simulate genomes:

wget http://misccb.googlecode.com/git-history/c910718ad6d53d7491c82a00d0d94f75263ac0a3/gatcdistsim.pl

using the uniform distribution assumption.
The values are quite close to the calculated expectation :)

ash022@ash022-Latitude-E4200:~/misccb$ perl gatcdistsim.pl 10000000

100000 254.581632653061 167 40
200000 260.997382198953 175 7
300000 249.43 172 128
400000 249.615 179.5 29
500000 257.480680061824 177 74
600000 253.092405063291 179 32
^C

PS: for running the above code to get the statistics, one needs http://search.cpan.org/dist/Math-NumberCruncher/ which can be installed on Ubuntu using:
sudo apt-get install libmath-numbercruncher-perl

de-Bruijn assembler

2012-01-27T01:40:00.000+05:30

In the last post I talked about the overlap-layout-consensus (OLC) way of Genome assembly. The approach which is (really!) getting popular these days is the other one, de-Bruijn-graphs (DBG). It is based on the simple idea of converting the hard problem of finding Hamiltonian to relatively simpler Eulerian for assembling biological sequences. There is a beautiful tutorial like introduction co-authored by the 'father' of this idea, Prof. Pevzner.

To put it naively, in contrast with OLC approach where one represents the reads as vertex in graph which are connected if the end of the read representing the first vertex overlaps with the start of the second one, DBG approach has the overlapping parts of the reads as the vertices and edges representing the reads consisting of the those overlaps. Now all we have to do, well kind of, is to traverse this graph so that we visit edges once, which is what Eulerian path is, and it scales up much better than Hamiltonian making it a cool choice for the cqncers churning out massive amount of shorter reads .

I made a simple scripts to see this approach in action. It can be downloaded via:

wget http://misccb.googlecode.com/git-history/585ec7811a837a1f20ee37256abada7d6ee538a3/eulerovlpath.pl

To test is, one needs a simple genome (here 100 bases long), which can be done via PERL one-liner:

perl -e '@b=qw/A T G C/;while($l<100){print @b[int(rand(4))];$l++;}' > rangen.txt

generating 100 base pair uniformly random DNA string. The sequencing can be simulated via breaking up the above string into small (here 10 bases long) overlapping (here of 9 bases long) strings using PERL one-liner:

perl -ne 'for($c=0;$c < length;$c++){print substr($_,$c,10);print "\n"};' rangen.txt > rangen.txt.kmer

This reads file can be directly fed into the downloaded eulerovlpath.pl:

perl eulerovlpath.pl rangen.txt.kmer

generating the string:

GGGCGTGGGCATTCCAACGACTGAGTCTGCCCTAATTCGTTGAAGTACTATTCGCAACGACCGCCGTAGACAAATAAGAGTGTAGGGATTACCGGCATT

which is similar to what we started with :)

cat rangen.txt
AGGGCGTGGGCATTCCAACGACTGAGTCTGCCCTAATTCGTTGAAGTACTATTCGCAACGACCGCCGTAGACAAATAAGAGTGTAGGGATTACCGGCATT

Happy hacking !

Hamiltonian Assembler

2012-01-05T20:33:00.001+05:30

Just came across Andrew's Hamiltonian Cycle finder (HCF)

get:

 wget https://raw.github.com/ahh/ahh-toys/master/ham.sh

simple check:

 bash ham.sh 
 a b
 b c
 c a 
 [followed by Ctrl-D should produce] 

a b c

based purely on shell commands and thought of testing it out as a Genome assembler. The genome assembly problem is closely related to finding the shortest common superstring (S) of a given set of strings (s1, s2… sn). The superstring S corresponds to the genome and the set of strings being the short sequence reads produced by the sequencing machines.

The popular ways of approaching this problem are OLC and DBG (overlap-layout-consensus and de-bruijn-graph). OLC approach is similar to HCF or the popular textbook example, the traveling-salesman-problem :).

Typical bacterial genomes are strings of nucleotides (A,T,G,C) and circular. General approach to find it is to extract it from the organism, make multiple copies, shatter it into small manageable pieces, feed into a sequencer and generate list of overlapping reads (since we made multiple copies, it is highly likely that reads coming out will be overlapping and the sum-read-length/genome-size represents coverage). In OLC approach, one represents the reads as vertex in graph which are connected if the end of the read representing the first vertex overlaps with the start of the second one.

This graph is traversed to find a Hamiltonian Path and here is where one can use Adam's HCF :) Thus wrote a simple PERL script:

wget http://misccb.googlecode.com/git-history/9deb0b79350f92aa41735091dcc69688a9a45bf0/seqmerge.pl

and ran it:

perl seqmerge.pl 9 3 4 2

to create random small genome of a given length (9), break it into overlapping k-mers of given size (3) reflecting a given coverage (4), use the overlap (2) information to create list of connected k-mers. This is fed to the HCF and we compare the Simulated Genome to the HCF generated Genome to see if it works. For example:

perl seqmerge.pl 9 3 4 2
Simulated Genome CATTTTGAT
Generated k-mers CAT ATT TTT TTT TTG TGA GAT ATC TCA CAT ATT TTT
Genome file Genome.L9.K3.C4.O2.fa
Connectivity file Genome.L9.K3.C4.O2.cg
Running shell Hamiltonian path finder and writing output to Genome.L9.K3.C4.O2.cg.ham
Collapsed ham.sh output from Genome.L9.K3.C4.O2.cg.ham written to Genome.L9.K3.C4.O2.fa.ham
Assembled Genome TTTGATCAT

and as we can see the circular genome CATTTTGAT is same as TTTGATCAT, so it did work :)

Getting out gene list from Ensembl web query results with WWW::Mechanize

2011-07-20T04:51:00.000+05:30

Seems like the Ensembl API is not providing a way to get list of results (?) for web query, thus decided to write a simple perl script. So to get list of gene hits for Ensembl query 'histone h2a histone cluster 1' one can simply run:

perl getensidmech.pl 'histone h2a histone cluster 1'

and the quick and dirty code getensidmech.pl using WWW::Mechanize is there at http://research-code-base-animesh.googlecode.com/svn/trunk/code_lic/getensidmech.pl :)

Cancer classification using Google prediction API

2011-05-29T18:12:00.000+05:30

Google's prediction API is now available to everyone with nice tutorial on how to get going :) so I decided to check its performance using gene expression values for 5 selected genes using online-feature-selection from the classic AML/ALL dataset (training/test).
The dataset needed to be formatted (training/test) and copied to Google storage to get the default scripts working.
Coming to the results, well the training step:

./oauth-train.sh mlani/top5_tr.csv

gave a performance of about 3/4:
./oauth-check-training.sh mlani/top5_tr.csv
{
"kind": "prediction#training",
"id": "mlani/top5_tr.csv",
"selfLink": "https://www.googleapis.com/prediction/v1.2/training/mlani/top5_tr.csv",
"modelInfo": {
"modelType": "classification",
"classificationAccuracy": 0.75
},
"trainingStatus": "DONE"
}not too bad, considering they are using some sort of statistical learning and sample is not the recommended 10x of feature set...
Coming to testing part:
./oauth-predict.sh mlani/top5_tr.csv "1122,178,847,33,1018"
{
"kind": "prediction#output",
"id": "mlani/top5_tr.csv",
"selfLink": "https://www.googleapis.com/prediction/v1.2/training/mlani/top5_tr.csv/predict",
"outputLabel": "ALL",
"outputMulti": [
{
   "label": "ALL",
   "score": 0.899588
},
{
   "label": "AML",
   "score": 0.100412
}
]
}
so seems like the first example went well, putting ALL as ALL with ~90% score :)
To get it over all the test examples, needed to write a simple PERL script:
perl check-pred-googleapi.pl top5_te.csv
ALLALLALLALLALLALLALLALLALLALLALLALLALLALLAML1
ALLAML2
ALLALLALLAMLAMLAMLAMLAMLALL3
AMLALL4
AMLAMLALL5
AMLAMLAMLGot 29 out of 34(accuracy:0.852941176470588)
which I guess is pretty good considering black box like usage :)
Happy predicting !

PERL one liner to generate 1000 base genome

2011-05-22T00:47:00.001+05:30

perl -e '@b=qw/A T G C/;print ">Genome\n";while($l<1000){print @b[int(rand(4))];$l++;}'

Finding coding region with Fourier transform

2011-03-29T19:30:00.001+05:30

Been a long time since I posted anything here, well I hardly had anything worth posting ...
Recently a friend of mine asked for implementation of Genescan code presented in "Prediction of probable genes by Fourier analysis of genomic sequences" / "Ab initio gene identification: prokaryote genome annotation with GeneScan and GLIMMER" and thus I went back looking into my old posts here... and strangely there were none about it! So decided to post for future reference.
Running the code waexome.pl is straight forward, for example running it over ins.fas:

perl waexome.pl ins.fas 3 300

will split sequence(s) into window of 300 and look for periodicity at 3, generate one .fted.txt for each sequence in the ins.fas and also ins.fted.txt which will report regions > p2nr threshold (hardcoded to 4), naively these are the coding regions according to the Fourier transformation threshold as they show peak periodicity at 3 (size of codon). Thus the generated file ins.fas.fted.txt reports the coding region for this particular example :).

Harnessing the power of multicore

2008-07-29T00:35:00.003+05:30

Most of us are working on multi - core processors , but the code we generally write do not utilize the extra cpu given to us. I do lot of "Align EST and genomic DNA sequences" using est2genome from EMBOSS. This is a task which is embarrassingly amenable to parallelization. So how do we go about doing this using PERL? Well following a posting at PerlMonks, I got hold of Parallel::ForkManager:
wget http://search.cpan.org/CPAN/authors/id/D/DL/DLUX/Parallel-ForkManager-0.7.5.tar.gz
and extracted to:
/scratch/misc/parallel/Parallel
[if you have root access to your machine, I would recommend to install this module instead of going for 'use lib'].
Once this is done, all you need is follow the examples and get going :)
For example, I have EST in fasta format files ending with .fas. The genome I am aligning to is NC_010336.fna. For this I can write the following code:


  use strict;
  use warnings;
  use lib '/scratch/misc/parallel/Parallel';
  use LWP::Simple;
  use Parallel::ForkManager;
  system("ls -1 *.fas > list.tmp");
  my $command="est2genome";
  my $genome="NC_010336.fna";
  my @tasks;
  open(F,"list.tmp");
  while(< F >){chomp;push(@tasks,$_);}
  close F;  
  my $tasksize= @tasks;
  print "There are #  $tasksize \n";
  my $pm = new Parallel::ForkManager($tasksize); 
  foreach my $task (@tasks) {
    $pm->start and next; # do the fork
    system("$command $task $genome $task.$genome.out");
    $pm->finish; 
  }
  $pm->wait_all_children;

and store is as est2genome_parallel.pl and run it as:
perl est2genome_parallel.pl
above code can be downloaded from http://sharma.animesh.googlepages.com/est2genome_parallel.pl, enjoy :)

Blast /usr/share/dict/linux.words

2008-07-28T23:47:00.006+05:30

After reading 'Waq’s Words and World', where the author tries to find anagrams by scanning through dictionary, I got interested in checking which is the biggest word in /usr/share/dict/linux.words which can give a significant hit with blasp against swissprot. So I vimmed dict.pl and wrote:


#!/usr/bin/perl
while(<>){chomp;split(/\s+/);$c++;$dictast{length(@_[0])}.=">s.$c\n@_[0]\n";}
foreach $w (sort {$b<=>$a} keys %dictast){
   open(FI,">temp.blast.in");
   print "Blasting $w length word(s) file\n";
   print FI"$dictast{$w}";
   close FI;
   system("blastcl3 -p blastp -d swissprot -i temp.blast.in -o temp.blast.out");
   open(FO,"temp.blast.out");
   while(<>){
           print $_;
           if($_=~/^Sequences producing significant/){
                   close FO;  
                   die;
           }
   }
}

The program died aligining SPECTROCOLORIMETRY with SPDCAERCGIMRLMDTRY from ">sp|Q0U3Y6|DDI1_PHANO DNA damage-inducible protein 1":


Query: 1          SPECTROCOLORIMETRY 18
                  SP+C   C + R+M+TRY
Sbjct: 229        SPDCAERCGIMRLMDTRY 246

To replicate or explore-further, the code can be downleaded from http://sharma.animesh.googlepages.com/dict.pl ad run as:
perl dict.pl /usr/share/dict/linux.words
[ /usr/share/dict/linux.words is the dictionary file coming with fedora 9, you may have it in different location, to find, do 'find / -name dict' ].

clipped from idazuwaika.wordpress.com

Anyway, to avoid getting rusty with programming (as I have with football), I wrote a Java program to scan through list of English words and find anagram pairs. Below are the interesting ones.

antagonist stagnation

ascertain sectarian

bacterial calibrate

coordinate decoration

courteous outsource

eroticism isometric

excitation intoxicate

prettiness persistent

satirical racialist

shattering straighten

supersonic percussion

Kyoto Prize goes to Karp

2008-07-27T01:40:00.001+05:30

Prof. Richard M. Karp [ http://www.cs.berkeley.edu/~karp/biography.html ] won the prestigious Kyoto Prize this year in Information Science category [ http://en.wikipedia.org/wiki/List_of_Kyoto_Prize_winners ].Good to see a Bioinformatics hero getting something like Nobel.

clipped from www.bizjournals.com

UC-Berkeley professor Richard Karp has won the 2008 Kyoto Prize, Japan's equivalent of the Nobel Prize, in recognition of his lifetime achievements in the field of computer theory.

Karp is one of three laureates named Friday, entitling him to a cash gift of about $460,000, a diploma and a gold medal recognizing his work in defining the field of theoretical computer science.

Switiching via common structural framework

2008-07-27T01:17:00.001+05:30

Josh explores [ http://www.thinkgene.com/researchers-reveal-insights-into-hidden-world-of-protein-folding/ ] an article from Booth et al. titled 'Mechanism of lid closure in the eukaryotic chaperonin TRiC/CCT' [ http://www.nature.com/nsmb/journal/v15/n7/abs/nsmb.1436.html or http://www.scribd.com/doc/3321492/Mechanism-of-lid-closure-in-the-eukaryotic-chaperonin-TRiCCCT ].

clipped from www.thinkgene.com

Virtually all proteins have to be folded-some in complex configurations-in order to function properly, and many are known to require a molecule called a chaperone to fold them. Frydman estimates that perhaps 10 percent of the proteins needing chaperones must have one that, like TRiC, is part of the subset called chaperonins. Other work done in Frydman’s lab has shown that proteins that have very complex folds seem to require chaperonins.

Biopython 1.47

2008-07-27T01:01:00.003+05:30

Hope you guys are convinced to explore Python for you job after reading http://www.bioinformaticszen.com/2008/02/bioinformatics-zen-faq/ and more recently Mark Bieda's resentment at http://markbieda.wordpress.com/2008/06/10/i-wish-i-had-started-with-python-earlier/ and some reasons and his ways at http://markbieda.wordpress.com/2008/06/18/python-for-perl-programmers-and-bioinformatics-people/ on how to go about it. I have one more reason for you to follow the wave, the biopython's release of version 1.4.7 [ http://biopython.org/wiki/Download ] ... I am also noticing a recent rise in interest for the language LittleB based on lisp [story via Eureka].

clipped from biopython.org

The latest release is Biopython 1.47, released on 5 July 2008. Get it from our Download page.

BZ: Bioinformatics Career Survey 2008

2008-07-27T00:48:00.001+05:30

Mike is running a survey about Bioinformatics Career [ http://www.bioinformaticszen.com/2008/07/creating-a-picture-of-different-careers-in-bioinformatics/ ]. Although it is running for quite some time, but I just came back to blog-world. I would appeal to all the people in the domain to fill this. The more the samples, the better the survey and more over this is time of data driven science [ http://www.wired.com/science/discoveries/magazine/16-07/pb_theory ] though Deepak has some issues with this development [ http://mndoci.com/blog/2008/06/25/chris-anderson-you-are-wrong/ ].

clipped from www.bioinformaticszen.com

When I wrote a post about my opinion of doing a career in bioinformatics I got the impression from the comments that this was something many people wanted more information about. Pedro and I had some discussion, and thought it might be interesting to create an online survey, to get current researchers’ opinions of working in the field of bioinformatics. So, the below survey begins today (July 1) and is filled out by as many people as possible over the course the next month. The data is then released into the public domain at the start of next month (August 1). Anyone who is then interested can contribute back analysis of this data, so that on September 1 hopefully we can compile together lots of interesting statistics and graphs into a handy document discussing the highs and lows of career in bioinformatics.

A Book Review from Adaptive Complexity: Microcosm E. coli and The New Science of Life

2008-07-27T00:14:00.001+05:30

After reading this review [ http://www.scientificblogging.com/adaptive_complexity/e_coli_as_biologys_decoder_ring? ] from Michael White [ http://www.scientificblogging.com/mwhite74/feed ] on a book from Carl Zimmer [ http://blogs.discovermagazine.com/loom/ ] titled 'Microcosm: E. coli and The New Science of Life' [ http://www.amazon.com/gp/search?ie=UTF8&keywords=Zimmer%20e.%20Coli&tag=funnierthanyo-20&index=books&linkCode=ur2&camp=1789&creative=9325 ], I am quite eager to read this book.
Presently I am reading a book recommended by Myers [ http://scienceblogs.com/pharyngula/ ] titled 'Developmental Plasticity and Evolution' [ http://www.amazon.com/gp/product/0195122356 ] by Mary Jane West-Eberhard [ http://www.stri.org/english/scientific_staff/staff_scientist/scientist.php?id=35 ]. This along with the article referred in the above review 'Evolutionary Dynamics of Prokaryotic Transcriptional Regulatory Networks' [ http://dx.doi.org/10.1016/j.jmb.2006.02.019 ] by Madan Babu and Eisen talking about one of the heroes of bioinformatics at http://phylogenomics.blogspot.com/2008/06/connection-between-video-games-and.html , I am pretty much convinced that Evolutionary biology has lot of opportunity for Bioinformatics guys to be busy for eternity, not just "500 years of exciting problems" [ http://tex.loria.fr/historique/interviews/knuth-clb1993.html ]!

clipped from www.scientificblogging.com

For such a slender book, Microcosm covers a wide-ranging selection of science. Zimmer begins by recapping the key events in the history of molecular biology, events in which E. coli was frequently a central player. Once scientists realized that E. coli had genes just like animals and plants, this gut bacterium gradually became one of the favorite model research systems in the then hot, new science of molecular biology. Experiments in E. coli revealed how genes are structured, how DNA is replicated, and how genes are controlled. Marshall Nirenberg and his colleagues cracked the essentially universal genetic code using cellular components from E. coli. Joshua Lederberg used bacterial 'sex' to bring the formidable tools of genetics to bacterial studies. Today, E. coli is the most well-mapped organism on the planet.

Last Lecture

2008-07-26T21:58:00.001+05:30

Randy Pausch [ http://en.wikipedia.org/wiki/Randy_Pausch ] passed away yesterday morning. I admire him a lot and whenever I feel down, I watch his Last Lecture [ http://www.youtube.com/watch?v=ji5_MqicxSo ]. May his soul rest in peace ...

clipped from www.thelastlecture.com

454

2008-06-20T14:56:00.004+05:30

I have been trying to understand the pyrosequencing technology, popularly known as 454. The very first question which came to my mind is WHY it is called 454? Is it because .454 kasool Handgun is sleeker then Shotgun ( Cf. Whole Genome Shotgun Sequence ) and slowly becoming as powerful ... or is it because of the fact that "In one instrument run sequence a minimum of 20 million base pairs in 4.5 hours" ... just a conjecture, if you guys have a proof, please share !
Anyways coming back to the technology, what I understand (below is a crude attempt to visualize it) is that there are 4 containers containing A,T,G and C which sequentially pass over the bowl containing a mixture of Polymerase, Sulfurylase, Luciferase, Apyrase and the DNA fragment to be sequenced. The addition of the nucleotide will light up ( thanks to lucifer-ase ) the CCD (a thing close to what we have in a digital camera and the costliest part of the instrument). Since we know the culprit nucleotide container, peak telling how many of that nucleotide is added, we can sequentially proceed and keep adding letters to the reads... this is happening in a massively parallel fashion and that brings out the power of 454...

Biocomputing

2008-06-14T18:14:00.001+05:30

Haynes et al. in the article "Engineering bacteria to solve the Burnt Pancake Problem" [ http://www.jbioleng.org/content/2/1/8#IDAJQH0E ] demonstrate how one can take advantage of direction sense embedded in DNA by coding 5'-3' as 0 and 3'-5' as 1 binary coding to convert cells into computers.

clipped from www.msnbc.msn.com

Ron Weiss, an assistant professor of electrical engineering and molecular biology at Princeton University, said the new study provides a “nice demonstration that there’s some capability to instruct cells to carry out computational tasks. And I think this certainly brings a new aspect to what’s been demonstrated before.” Weiss, who wasn’t involved with the research, said the effort is another sign that the field is progressing, despite its relative infancy. “We really are at the beginning, at the vacuum tube stage or something like that, if you use the comparison with computer electronics,” he said.

International Conference on Functional Programming 2008

2008-06-07T21:44:00.002+05:30

The 11th ICFP ( http://icfpcontest.org/ ) Programming Contest is from Friday, July 11, 2008 to Monday, July 14, 2008. More details on what the contest is like, check out http://cpoucet.wordpress.com/2008/06/11/icfp-contest-2008/ .
Mail queries to Tim Sheard at sheard@cs.pdx.edu.

clipped from lambda-the-ultimate.org

The ICFP Programming Contest is one of the most advanced and
prestigious programming contests, as well as being a chance to show
off your programming skills, your favorite languages and tools, and
your ability to work as a team. The contest is affiliated with the
International Conference on Functional Programming. Teams consisting
of one or more participants, from any part of the world, using any
programming language, may enter.

Unofficial Google Shell

2008-06-03T18:45:00.003+05:30

I have been waiting for google shell. I thought they will name it gshell, but it turns out:
goosh.org - the unofficial google shell ( http://goosh.org/ ).
It is pretty neat and useful in its neonatal state itself, like to read a feed, you can type, for eg:

r http://www.google.com/reader/public/atom/user/11136726768885096694/label/cbn-roll

and gets default 4 latest feeds, or you can check out a ncbi accession by typing, for eg:

NP_001108376

and get to the Fugu fish refseq link.

The most useful thing is the fact that it it text based and quite easy to parse visually as well as via scripts.
Another cool feature for me right now is the shell translation, like to get the english word for the norwegian word 'hvor', all I have to do is type:

t no en hvor

which correctly says:
translating "hvor" from "no" to "en":

"where"

More discussion going on at http://tech.slashdot.org/article.pl?sid=08/06/02/222234&from=rss .
Enjoy!

clipped from goosh.org

Goosh goosh.org 0.4.3-beta #1 Mon, 02 Jun 08 22:28:01 UTC Google/Ajax

help

command aliases parameters function
web (search,s,w) [keywords] google web search
lucky (l) [keywords] go directly to first result
images (image,i) [keywords] google image search
wiki (wikipedia) [keywords] wikipedia search
clear (c) 
clear the screen
help (man,h,?) [command] displays help text
news (n) [keywords] google news search
blogs (blog,b) [keywords] google blog search
feeds (feed,f) [keywords] google feed search
open (o) <url> open url in new window
go (g) <url> open url
more (m) 
get more results
in (site) <url> <keywords> search in a specific website
load 
<extension_url> load an extension
video (videos,v) [keywords] google video search
read (rss,r) <url> read feed of url
place (places,map,p) [address] google maps search
lang 
<language> change language
addengine 

add goosh to firefox search box
translate (trans,t) [lang1] [lang2] <words> google translation

Conserved Domain Database (CDD) has been updated

2008-06-01T16:59:00.001+05:30

The CDD [ http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml ] and the its search tool [ http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi ] has been recently updated. Also check out the latest approach for "automatically assigning subcellular locations" to protein from Newberg and Murphy [ http://pubs.acs.org/cgi-bin/abstract.cgi/jprobs/asap/abs/pr7007626.html ], more at http://www.sciencedaily.com/releases/2008/05/080512145232.htm .

clipped from beckerinfo.net