哈希错误跟踪计数，运行时间长

Question

我正在使用 Perl 编写一个程序，但我的输出是错误的，并且需要很长时间才能处理。该代码旨在接收一个大型 DNA 序列文件，以 15 个字母增量 (kmers) 读取它，一次向前走 1 个位置。我应该将 kmer 序列输入到哈希中，它们的值是该 kmer 的出现次数——这意味着每个键应该是唯一的，当找到重复项时，它应该增加该特定 kmer 的计数。我从我的教授预期的输出文件中知道，我有太多的行，所以它允许重复并且不正确计数。它也运行了 5 分钟以上，所以我必须 Ctrl+C 才能逃脱。当我查看 kmers.txt 时，该文件至少已正确写入和格式化。

#!/usr/bin/perl

use strict;
use warnings;
use diagnostics;

# countKmers.pl
# Open file /scratch/Drosophila/dmel-2L-chromosome-r5.54.fasta
# Identify all k-mers of length 15, load them into a hash
# and count the number of occurences of each k-mer. Each
# unique k-mer and its' count will be written to file
# kmers.txt

#Create an empty hash
my %kMersHash = ();

#Open a filehandle for the output file kmers.txt
unless ( open ( KMERS, ">", "kmers.txt" ) ) {
    die $!;
}

#Call subroutine to load Fly Chromosome 2L
my $sequenceRef = loadSequence("/scratch/Drosophila/dmel-2L-chromosome-r5.54.fasta");

my $kMer      = 15;    #Set the size of the sliding window

my $stepSize  = 1;     #Set the step size

for (

    #The sliding window's start position is 0
    my $windowStart = 0;

    #Prevent going past end of the file
    $windowStart <= ( length($$sequenceRef) - $kMer );

    #Advance the window by the step size
    $windowStart += $stepSize

    )

{

    #Get the substring from $windowStart for length $kMer
    my $kMerSeq = substr( $$sequenceRef, $windowStart, $kMer );

#Call the subroutine to iterate through the kMers
    processKMers($kMerSeq);

}

sub processKMers {

    my ($kMerSeq) = @_;

    #Initialize $kCount with at least 1 occurrence  
    my $kCount = 1;

    #If the key already exists, the count is  
    #increased and changed in the hash
    if ( not exists $kMersHash{$kMerSeq} ) {

            #The hash key=>value is loaded: kMer=>count
            $kMersHash{$kMerSeq} = $kCount;
    }

    else {

            #Increment the count 
            $kCount ++;

            #The hash is updated 
            $kMersHash{$kMerSeq} = $kCount;
    }

    #Print out the hash to filehandle KMERS
    for (keys %kMersHash) {
            print KMERS $_, "\t", $kMersHash{$_}, "\n";
    }
}

sub loadSequence {

    #Get my sequence file name from the parameter array
    my ($sequenceFile) = @_;

    #Initialize my sequence to the empty string
    my $sequence = "";

    #Open the sequence file
    unless ( open( FASTA, "<", $sequenceFile ) ) {
            die $!;
    }

    #Loop through the file line-by-line
    while (<FASTA>) {

            #Assign the line, which is in the default 
            #variable to a named variable for readability.
            my $line = $_;

            #Chomp to get rid of end-of-line characters
            chomp($line);

            #Check to see if this is a FASTA header line
            if ( $line !~ /^>/ ) {

                    #If it's not a header line append it 
                    #to my sequence
                    $sequence .= $line;
            }

    }

    #Return a reference to the sequence
    return \$sequence;
}

Answer 1

以下是我将如何编写您的应用程序。 processKMers 子例程归结为只是增加一个哈希元素，所以我已经删除了它。我还更改了标识符以匹配在 Perl 代码中更常见的 snake_case，并且我没有看到 load_sequence 返回对序列的引用的任何点，因此我将其更改为返回字符串本身

use strict;
use warnings 'all';

use constant FASTA_FILE => '/scratch/Drosophila/dmel-2L-chromosome-r5.54.fasta';
use constant KMER_SIZE  => 15;
use constant STEP_SIZE  => 1;

my $sequence = load_sequence( FASTA_FILE );

my %kmers;

for (my $offset = 0;
        $offset + KMER_SIZE <= length $sequence;
        $offset += STEP_SIZE ) {

    my $kmer_seq = substr $sequence, $start, KMER_SIZE;

    ++$kmers{$kmer_seq};
}

open my $out_fh, '>', 'kmers.txt' or die $!;

for ( keys %kmers ) {
    printf $out_fh "%s\t%d\n", $_, $kmers{$_};
}

sub load_sequence {

    my ( $sequence_file ) = @_;

    my $sequence = "";

    open my $fh, '<', $sequence_file or die $!;

    while ( <$fh> ) {
        next if /^>/;
        chomp;
        $sequence .= $_;
    }

    return $sequence;
}

这是增加哈希元素的一种更简洁的方法，无需直接在哈希上使用++

my $n;

if ( exists $kMersHash{$kMerSeq} ) {
    $n = $kMersHash{$kMerSeq};
}
else {
    $n = 0;
}

++$n;
$kMersHash{$kMerSeq} = $n;

Answer 2

除了processKMers 之外，您的代码中的一切看起来都很好。主要问题：

• $kCount 在对 processKMers 的调用之间不持久，因此在您的 else 语句中，$kCount 将始终为 2

• 每次调用 processKMers 时都会打印，这会拖慢您的速度。经常打印会显着减慢您的进程，您应该等到程序结束再打印一次。

保持您的代码基本相同：

sub processKMers {

    my ($kMerSeq) = @_;

    if ( not exists $kMersHash{$kMerSeq} ) {
            $kMersHash{$kMerSeq} = 1;
    }
    else {
            $kMersHash{$kMerSeq}++;
    }
}

然后你想将你的打印逻辑移动到你的 for 循环之后。