Perl Expectation Максимизация баллов шансов последовательностей ДНК - PullRequest
Новая
       18

Perl Expectation Максимизация баллов шансов последовательностей ДНК

0 голосов
/ 20 марта 2019

Моя общая цель этого кода - найти мотив (меньшую последовательность) в каждой последовательности ДНК, который будет сообщать максимальную оценку log-odds на основе матрицы оценки log-odds.Файл .txt, который я ищу, выглядит следующим образом: 

>Sequence 1
TCGACGATCAGACAG
>Sequence 2
TGGGACTTGCACG

.... и т. Д.

В данный момент я работаю над этапом максимизации своего кода, иЯ изо всех сил пытаюсь вычислить логарифмическую оценку мотива в моей последовательности ДНК.У меня есть код, который создает матрицу оценки лог-шансов - см. Следующее: 

#!/usr/bin/perl -w
#Above line used in Unix and Linux
#Grome Programming Assignment 2
#E-M
use strict;
use warnings;

# usage: script.pl {motif_width} {dnafile}
#USER SPECIFICATIONS
print "Please enter the filename of the fasta sequence data: ";
my $filename1 = <STDIN>;

#Remove newline from file
chomp $filename1;

#Open the file and store each dna seq in hash
my %id2seq = ();
my $id = '';
open (FILE, '<', $filename1) or die "Cannot open $filename1.",$!;
my $dna;
while (<FILE>)
{
    if($_ =~ /^>(.+)/)
    {
        $id = $1;
    }
    else
    {
        $id2seq{$id} .= $_;
    }
}
close FILE;
foreach $id (keys %id2seq)
{
    print "$id2seq{$id}\n\n";
}

#User specifies motif width
print "Please enter the motif width:\n";
my $width = <STDIN>;

#Remove newline from file
chomp $width;

#Default width is 3 (arbitrary number chosen)
if ($width eq '')
{
    $width = 3;
}
elsif ($width <=0)
{
    print "Please enter a number greater than zero:\n";
    $width = <STDIN>;
    chomp $width;
}

#User specifies number of initial random
#starting alignments
print "Please enter the number of initial random
starting alignments:\n";
my $start = <STDIN>;

#Remove newline from file
chomp $start;

#Default start is 50
if ($start eq '')
{
    $start = 50;
}
elsif ($start <=0)
{
    print "Please enter a number greater than zero:\n";
    $start = <STDIN>;
    chomp $start;
}

#User specifies number of iterations to
#perform expectation-maximization
print "Please enter the number of iterations for 
expectation-maximization:\n";
my $iteration = <STDIN>;

#Remove newline from file
chomp $iteration;

#Default iteration = 500
if($iteration eq '') 
{
    $iteration = 500;
}
elsif ($iteration <=0)
{
    print "Please enter a number greater than zero:\n";
    $iteration = <STDIN>;
    chomp $iteration;
}

#EXPECTATION
#Initialize counts for motif positions
#Incorporate pseudocounts initially
my %mot = map { $_ => [ (1) x $width ] } qw( A C G T );

# Initialize background counts
my %bg = map { $_ => 0 } qw( A C G T );

#Fill background and motif counts
foreach $id (keys %id2seq)
{
    #Generate random start site in the sequence
    #for motif to start from
    my $ms = int(rand(length($id2seq{$id})-$width));

    # Within a motif, count the bases at the positions
    for my $pos (0..length($id2seq{$id})-1)
    {
        my $base = substr($id2seq{$id}, $pos, 1);
        if ($pos >= $ms && $pos < $ms + $width) 
        {
            ++$mot{$base}[$pos-$ms]
                if exists($mot{$base});
        } 
        else 
        {
            ++$bg{$base}
            if exists($bg{$base});
        }
    }
}
#Print the background and motif counts
for my $base (qw( A C G T )) 
{
   print "$base @{$mot{$base}}\n";
}

print "\n";

for my $base (qw( A C G T )) 
{
   print "bg$base = $bg{$base}\n";
}

#Create frequency table of the motifs
#Get sum of the background
my $bgsum = 0;
for my $base (qw( A C G T))
{
    $bgsum = $bgsum + $bg{$base};
}
print "\n$bgsum\n\n";

#Create background frequency table
my %bgfreq = map { $_ => 0 } qw( A C G T );
for my $base (qw( A C G T))
{
    $bgfreq{$base} = $bg{$base} / $bgsum;
    print "bgfreq$base = $bgfreq{$base}\n";
}

#Get sum of each motif position
my @motsum = ( (0) x $width );
for my $base (qw( A C G T))
{
    for my $arrpos (0.. ($width-1))
    {
        $motsum[$arrpos] = $motsum[$arrpos] + @{$mot{$base}}[$arrpos];
    }
}

#Create motif frequency table
my %motfreq = map { $_ => [ (0) x $width ]} qw( A C G T );
for my $base (qw( A C G T))
{
    for my $arrpos (0.. ($width-1))
    {
        $motfreq{$base}[$arrpos] = $mot{$base}[$arrpos] / $motsum[$arrpos];
    }
    print "motfreq$base @{$motfreq{$base}}\n";
}

#Create odds table of motifs
my %odds = map { $_ => [ (0) x ($width) ]} qw( A C G T );

for my $base (qw( A C G T))
{
    for my $arrpos (0.. ($width-1))
    {
        $odds{$base}[$arrpos] = $motfreq{$base}[$arrpos] / $bgfreq{$base};
    }
    print "odds$base @{$odds{$base}}\n";
}

#Create log-odds table of motifs
my %logodds = map { $_ => [ (0) x ($width) ]} qw( A C G T );
for my $base (qw( A C G T))
{
    for my $arrpos (0.. ($width-1))
    {
        $logodds{$base}[$arrpos] = log2($odds{$base}[$arrpos]);
    }
    print "logodds$base @{$logodds{$base}}\n";
}
#####################################################
sub log2
{
    my $n = shift;
    return log($n)/log(2);
}

Теперь мне нужно вычислить балл лог-шансов для мотива в каждой последовательности ДНК.Затем я переберу все возможные позиции в последовательности и найду максимальную оценку для каждой последовательности.Дальнейшая работа требует от меня вспомнить мотив с максимальным количеством баллов, но я еще не пытался это сделать (просто хотел выделить область для этих максимальных баллов).

Стратегия: я собираюсь создать хешиз баллов log-odds и максимальных баллов, чтобы удерживать максимальные баллы каждой последовательности как итеративные.Чтобы вычислить оценку log-odds, я посмотрю, какой элемент матрицы оценки log-odds соответствует элементам в мотиве.

#MAXIMIZATION
#Determine location for each sequence that maximally
#aligns to the motif pattern
#Calculate logodds for the motif

#Create hash of logodds scores and hash of maxscores
#so each id has a logodds score and max score
 my %loscore = map { $_ => [ (0) x (length($id2seq{$id})-$width) ]} qw( $id ); #Not sure if $id is correct, but I want a loscore for each $id
 my %maxscore = map { $_ => [ (0) x (length($id2seq{$id})-$width) ]} qw( $id ); #Not sure if $id is correct, but I want a maxscore for each $id
 foreach $id (keys %loscore, %maxscore, %id2seq)
 {
 my $len = length($id2seq{$id});
 for my $base (qw( A C G T ))
 {
     for my $pos (0..$len-1)
     {
         if ($id2seq{$id}[$pos] = $mot{$base})
         {
             for my $motpos (0..$width-1)
             {
                 $loscore{$id} = $loscore{$id} + $logodds{$base}[$motpos];
                 if ($loscore{$id} > $maxscore{$id})
                 {
                     $maxscore{$id} = $loscore{$id};
                 }
             }
         }
     }
 }
 print "$id2seq{$id} logodds score: $maxscore{$id}\n";
 }

#####################################################
sub log2
{
    my $n = shift;
    return log($n)/log(2);
}

Когда я раскомментирую шаг максимизации в коде и запускаю его, из секции максимизации ничего не печатается.Я не получаю никаких ошибок, но ничего нового печатает.Я понимаю, что мой шаг ожидания может быть упрощен (я все уберу после того, как он сработает), но сначала я сосредоточусь на этом шаге максимизации.Я знаю, что в коде максимизации есть много недостатков, особенно при попытке создать оценку лог-шансов и максимальную оценку хэшей.Любой ввод помогает!Заранее благодарю за терпение и советы.Дайте мне знать, если вам нужны какие-либо разъяснения.

1 Ответ

0 голосов
/ 20 марта 2019

foreach $id (keys %loscore, %maxscore, %id2seq) не перебирает все три хэша. Вы, вероятно, имели в виду foreach $id (keys %loscore, keys %maxscore, keys %id2seq)

if ($id2seq{$id}[$pos] = $mot{$base}) присваивает значение $id2seq{$id}[$pos]. Чтобы проверить на eqaulity вам нужно if ($id2seq{$id}[$pos] eq $mot{$base}). Но, вероятно, это тоже неправильно, поскольку $id2seq{$id} должна быть строкой, а $mot{$base} - массивом.

Не ясно, как должен работать ваш код, и трудно найти ошибки в таких длинных примерах кода.

...