Sul suo blog personale, Dalloliogm ha scritto un interessante articolo su come python può facilitare la vita dei bioinformatici . Il suo è un bel esempio di come, attraverso una serie di oggetti legati alla genetica, alle banche dati, ecc, è possibile creare una libreria sempre pronta per elaborazioni comode e veloci: un meta linguaggio specifico per le ricerche che il bioinformatico deve svolgere.
Diciamolo subito: amo python. Amo il fatto che non necessiti di parentesi e che lavori esclusivamente ad oggetti. Lo amo per la sua sinteticità! In poche righe è possibile fare grandi operazioni.

Vi propongo un breve esempio comparativo dello sviluppo di una medesima funzione prima in Perl e poi in Python. Tenete presente che io non sono un guru di nessuno dei due linguaggi, e sicuramente le routine che vi proporrò possono essere ottimizzate. Anzi: siete invitati ad ottimizzarle!
La questione è la seguente. La prima operazione che un bioinformatico interessato alla proteomica deve sviluppare è un piccolo algoritmo per svolgere una “digestione in silico delle proteine“. Ovvero: data la sequenza aminoacidica di una proteina, devo poterla spezzare in tanti peptidi attraverso un enzima. Per facilitare la cosa noi simuleremo l’uso della tripsina. Enzima che taglia (FRAMMENTA) la sequenza seguendo le seguenti regole:

1. Taglia la sequenza quando trova un aminoacido K o R.
2. Non lo taglia nel caso ci sia una P successiva, ovvero non nel caso KP e RP.
3. E’ possibile ci sia un “errore di taglio“, ovvero è possibile che l’enzima non funzioni nel 100% dei casi e quindi che lasci la sequenza parzialmente integra

Un esempio per capirci:
La sequenza della nostra proteina sarà: AAAAAKBBKPCCCCRDDDKEEEEPFFRKGGGG
la mia routine deve restituire la seguente lista:
AAAAAK
BBKPCCCR
DDDK
EEEEPFFR
K
GGGG
AAAAAKBBKPCCCR
BBKPCCCRDDDK
DDDKEEEEPFFR
EEEEPFFRK
KGGGG

Ed ecco come faccio la magia in Perl (in input la sequenza e il numero di possibili errori di tagli (miscleavages)). In sostanza la sub percorre la sequenza, taglia prima i singoli pezzi e poi dalla lista creata unisce 2 a 2 per avere le sequenze con miscleavage.
###########################################

sub tryp {
my $seq = shift @_;
my $miscleav = shift @_;
my %Hfrag =();
my @fragments = @fragReturn = ();
my @splittedSeq=split(//,$seq);
my @fragIndex = ();
my $start=0;
my $end=0;
for ($i=0; $i<$#splittedSeq; $i++){
if (($splittedSeq[$i] =~ /[KR]/) && ($splittedSeq[$i+1] =~ /[^P]/)) {
$end = $endf = $i;
@fragment = @splittedSeq[$start..$end];
push @fragments, join(“”,@fragment);
push @fragIndex, $start;
$start = $i+1;
}
}
if ($endf<$#splittedSeq) {
push @fragments, join (“”, @splittedSeq[($endf+1)..$#splittedSeq]);
push @fragIndex, $endf+1;
}
if ($miscleav != 0) {
for ($p=1; $p<=$miscleav; $p++){
for ($j=0; $j<=$#fragIndex-$p-1; $j++) {
@fragment = @splittedSeq[$fragIndex[$j]..$fragIndex[$j+$p+1]-1];
push @fragments , join(“”,@fragment);
}
@fragment = @splittedSeq[$fragIndex[$#fragIndex-$p]..$#splittedSeq];
push @fragments , join(“”,@fragment);
}
}
foreach my $frdef (@fragments) { $Hfrag{$frdef}=1;}
foreach my $frdef (keys %Hfrag) { push @fragReturn, $frdef; }
return @fragReturn;
}

######################################

la stessa cosa scritta in python:

class Protein(object):
def __init__(self,seq):
self.sequence=seq
def dig_tryp(self,misscleavages=0):
seqKP=re.sub(‘KP’,’1′,self.sequence)
seqRP=re.sub(‘RP’,’2′,seqKP)
seqModified=re.sub(‘K’,'K ‘,seqRP)
invseqKP=re.sub(’1′,’KP’,seqModified)
invseqRP=re.sub(’2′,’RP’,invseqKP)
listPeptides=invseqRP.split()
list=listPeptides
listTot=listPeptides
for num in range(1,misscleavages):
listTemp=[]
listTemp=[list[y]+listPeptides[y+num] for y,t in enumerate(list[:len(list)-1]) if listPeptides[y+num] ]
list=listTemp
listTot.extend(list)
return listTot

##################################

Ho semplicemente creato una classe Protein e un oggetto dig_tryp.
Come si può notare la sinteticità del secondo esempio facilita anche la lettura del codice. E anche quando lo si usa
basta a questo punto creare una istanza con qualcosa del tipo:
P1 = Protein(AAAAAKBBKPCCCCRDDDKEEEEPFFRKGGGG)
peptidi = P1.dig_tryp(1);
Naturalmente anche Perl può essere programmato ad oggetti. E’ possibile creare package pm che contengano le varie sub da richiamare in altri script, ma Perl non è nato come un linguaggio Object oriented e quindi ne eredita certe limitazioni.
Altro vantaggio non da poco del python: il 99.9% dei moduli che vi potranno mai servire sono integrati nel pacchetto di base. Non vi dovrete mai spezzare la testa per far funzionare il modulo DBI su un sistema SUN con Solaris 10 (per esempio) e connettervi al vostro bel DB mysql!
Una piccola nota positiva per i Perlisti convinti (giusto per accontentare anche loro. L’unico vero vantaggio non da poco che vedo in Perl per applicazioni bioinformatiche è l’uso comodo e veloce delle espressioni regolari, che in python non è così immediato.

Immaginate ora una libreria di librerie, scritte in python, adatte per creare classi e oggetti per ogni settore scientifico. Integrabili, portabili! Altro che Bioperl!

Related Posts

1. Protein ID mapping (come passar da una nomenclatura ad un’altra in due facili click)
2. Genome Commons (Farsi il genoma da soli)
3. I 7 peccati capitali commessi dai bioinformatici.