2010-04-06 5 views
5

Sto cercando di eseguire un filtro basato su composizione su una vasta raccolta di stringhe (sequenze di proteine).
Ho scritto un gruppo di tre subroutine per occuparmene, ma mi trovo nei guai in due modi: uno minore, uno maggiore. Il guaio minore è che quando uso List::MoreUtils 'pairwise' ricevo avvertenze sull'utilizzo di $a e $b una sola volta e che non sono inizializzate. Ma credo che sto chiamando questo metodo correttamente (basato sulla voce di CPAN per esso e alcuni esempi dal web).
Il problema principale è un errore "Can't use string ("17/32") as HASH ref while "strict refs" in use..."Perché questa istruzione viene considerata come una stringa invece del suo risultato?

Sembra che questo può avvenire solo se il loop foreach in &comp dà i valori hash come una stringa invece di valutare l'operazione di divisione. Sono sicuro di aver commesso un errore da principiante, ma non riesco a trovare la risposta sul web. La prima volta che ho visto il codice perl era lo scorso mercoledì ...

use List::Util; 
use List::MoreUtils; 
my @alphabet = (
'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 
'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V' 
); 
my $gapchr = '-'; 
# Takes a sequence and returns letter => occurrence count pairs as hash. 
sub getcounts { 
my %counts =(); 
foreach my $chr (@alphabet) { 
    $counts{$chr} = ($_[0] =~ tr/$chr/$chr/); 
} 
$counts{'gap'} = ($_[0] =~ tr/$gapchr/$gapchr/); 
return %counts; 
} 

# Takes a sequence and returns letter => fractional composition pairs as a hash. 
sub comp { 
my %comp = getcounts($_[0]); 
foreach my $chr (@alphabet) { 
    $comp{$chr} = $comp{$chr}/(length($_[0]) - $comp{'gap'}); 
} 
return %comp; 
} 

# Takes two sequences and returns a measure of the composition difference between them, as a scalar. 
# Originally all on one line but it was unreadable. 

sub dcomp { 
my @dcomp = pairwise { $a - $b } @{ values(%{ comp($_[0]) }) }, @{ values(%{ comp($_[1]) }) }; 
@dcomp = apply { $_ ** 2 } @dcomp; 
my $dcomp = sqrt(sum(0, @dcomp))/20; 
return $dcomp; 
} 

Molto apprezzamento per qualsiasi risposta o consiglio!

+0

Si prega di consultare anche http://stackoverflow.com/questions/2261871/whats-the-point-of-use-vars-in-this-perl-subroutine/2261957#2261957 per ulteriori spiegazioni sul perché/come questa roba accade. +1 – DVK

+0

Su quale riga si trova l'errore di hash hash? – DVK

+0

Per il problema minore, vedere http://stackoverflow.com/questions/1490505/how-do-i-prevent-listmoreutils-from-warning-about-using-a-and-b-only-once. – FMc

risposta

2

%{ $foo } considera $foo come riferimento hash e dereferenziamento; allo stesso modo, @{} assegnerà riferimenti di matrice. Poiché comp restituisce un hash come elenco (gli hash diventano elenchi quando si passa ae da funzioni) e non un riferimento hash, lo %{} è errato. Potresti potenzialmente interrompere lo %{}, ma lo values è un modulo speciale e richiede un hash, non un hash passato come elenco. Per passare il risultato di comp a values, comp è necessario restituire un riferimento hash che viene quindi rimosso.

C'è un altro problema con il vostro dcomp, vale a dire che l'ordine di values (come il documentation dice) "vengono restituiti in un ordine apparentemente casuale", in modo che i valori passati al blocco pairwise non sono necessariamente per lo stesso carattere. Invece di values, puoi usare le fette di hash. Ora torniamo a comp restituendo un hash (come elenco).

sub dcomp { 
    my %ahisto = comp($_[0]); 
    my %bhisto = comp($_[1]); 
    my @keys = uniq keys %ahisto, keys %bhisto; 
    my @dcomp = pairwise { $a - $b } , @ahisto{@keys}, @bhisto{@keys}; 
    @dcomp = apply { $_ ** 2 } @dcomp; 
    my $dcomp = sqrt(sum(0, @dcomp))/20; 
    return $dcomp; 
} 

Questo non affronta cosa succede se un personaggio appare in uno solo dei $_[0] e $_[1].

uniq lasciato come esercizio per il lettore.

+0

Penso che le slice hash debbano essere racchiuse da @ {[]}? Ora conosco il pericolo dei valori(), grazie. –

1

re: problema minore

Questo è bene e un problema comune con (alcuni) dei moduli List::Util e List::MoreUtils.

Un modo per rimuovere gli avvisi è semplicemente dichiarare quelli special variables in anticipo, in questo modo:

our ($a, $b); 

Un'altra potrebbe essere quella precedere il pairwise con:

no warnings 'once'; 

Vedere perlvar per ulteriori informazioni circa $ a e $ b

/I3az/

+0

Vuoi fare quella dichiarazione nel sotto-blocco o nel, uh, spazio dei nomi globale (non sono così bravo nei termini corretti)? Mi sembra il primo se uso il metodo una volta e quest'ultimo se lo uso in molti posti? –

+0

Mettendo 'il nostro ($ a, $ b)' all'inizio del programma va bene in questo caso. $ a & $ b sono speciali e quindi * non dovrebbero * essere usati per altri tipi di barre e moduli come List :: Util. Il problema è che "warnings once" non tiene conto di variabili speciali $ a & $ b ad eccezione di 'sort' – draegtun

4

Ci sono alcuni errori nel codice. In primo luogo, si noti da perldoc perlop:

Poiché la tabella di traslitterazione è costruito a tempo di compilazione, né il SEARCHLIST né il REPLACEMENTLIST sono sottoposti ad una doppia citazione interpolazione.

Quindi il metodo di conteggio non è corretto. Credo anche che tu stia abusando di pairwise. È difficile valutare quale sia l'uso corretto perché non si forniscono esempi di quale output si dovrebbe ottenere con alcuni input semplici.

In ogni caso, vorrei riscrivere lo script come (ci sono alcune dichiarazioni di debug cosparsi in):

#!/usr/bin/perl 

use List::AllUtils qw(sum); 
use YAML; 

our ($a, $b); 
my @alphabet = ('A' .. 'Z'); 
my $gap = '-'; 

my $seq1 = 'ABCD-EFGH--MNOP'; 
my $seq2 = 'EFGH-ZZZH-KLMN'; 

print composition_difference($seq1, $seq2); 

sub getcounts { 
    my ($seq) = @_; 
    my %counts; 
    my $pattern = join '|', @alphabet, $gap; 
    $counts{$1} ++ while $seq =~ /($pattern)/g; 
    warn Dump \%counts; 
    return \%counts; 
} 

sub fractional_composition_pairs { 
    my ($seq) = @_; 
    my $comp = getcounts($seq); 
    my $denom = length $seq - $comp->{$gap}; 
    $comp->{$_} /= $denom for @alphabet; 
    warn Dump $comp; 
    return $comp; 
} 

sub composition_difference { 
    # I think your use of pairwise in the original script 
    # is very buggy unless every sequence always contains 
    # all the letters in the alphabet and the gap character. 
    # Is the gap character supposed to factor in the computations here? 

    my ($comp1, $comp2) = map { fractional_composition_pairs($_) } @_; 
    my %union; 
    ++ $union{$_} for (keys %$comp1, keys %$comp2); 

    my $dcomp; 
    { 
     no warnings 'uninitialized'; 
     $dcomp = sum map { 
      ($comp1->{$_} - $comp2->{$_}) ** 2 
     } keys %union; 
    } 

    return sqrt($dcomp)/20; # where did 20 come from? 
} 
+0

Senza entrare nella biologia (non sembra il posto!), Le lacune di sequenza non costituiscono differenze nella "composizione" delle sequenze, perché non sono una lettera dell'alfabeto (dei 20 aminoacidi). Il 20 deve (arbitrariamente) normalizzare $ dcomp da 0 a 1. Ottengo la mappa() adesso ... grazie! –

2

solo passando attraverso il codice che hai fornito, questo è come l'avrei scritto . Non so se funzionerà come volevi che funzionasse.

use strict; 
use warnings; 
our($a, $b); 

use List::Util; 
use List::MoreUtils; 

my @alphabet = split '', 'ARNDCQEGHILKMFPSTWYV'; 
my $gapchr = '-'; 
# Takes a sequence and returns letter => occurrence count pairs as hash. 
sub getcounts { 
    my($sequence) = @_; 
    my %counts; 
    for my $chr (@alphabet) { 
    $counts{$chr} =() = $sequence =~ /($chr)/g; 
    #() = forces list context 
    } 
    $counts{'gap'} =() = $sequence =~ /($gapchr)/g; 
    return %counts if wantarray; # list context 
    return \%counts; # scalar context 
    # which is what happens inside of %{ } 
} 

# Takes a sequence and returns letter => fractional composition pairs as a hash 
sub comp { 
    my($sequence) = @_; 
    my %counts = getcounts($sequence); 
    my %comp; 
    for my $chr (@alphabet) { 
    $comp{$chr} = $comp{$chr}/(length($sequence) - $counts{'gap'}); 
    } 
    return %comp if wantarray; # list context 
    return \%comp; # scalar context 
} 

# Takes two sequences and returns a measure of the composition difference 
# between them, as a scalar. 
sub dcomp { 
    my($seq1, $seq2) = @_; 
    my @dcomp = pairwise { $a - $b } 
    @{[ values(%{ comp($seq1) }) ]}, 
    @{[ values(%{ comp($seq2) }) ]}; 
    # @{[ ]} makes a list into an array reference, then dereferences it. 
    # values always returns a list 
    # a list, or array in scalar context, returns the number of elements 
    # ${ } @{ } and %{ } forces their contents into scalar context 

    @dcomp = apply { $_ ** 2 } @dcomp; 
    my $dcomp = sqrt(sum(0, @dcomp))/20; 
    return $dcomp; 
} 

Una delle cose più importanti da sapere sono le differenze tra contesti scalari, elenco e vuoto. Questo perché tutto ciò che è si comporta in modo diverso, nei diversi contesti.

+0

Ah, sì, vedo i contesti. 'se wantarray' è un trucco pulito. –

Problemi correlati