查找并打印所有重叠的 k-mer

Question

我正在尝试编写一个 perl 程序，它读取一个 fasta 文件并打印出一个文本文件，其中包含来自序列 (fasta) 文件的所有可用（重叠）长度为 15 k-mers 的文本文件。当我搜索不重叠的 k-mer 时，该程序运行得非常好，但是当我对其进行编码以找到重叠的 k-mer 时，它需要永远执行，并且 Cygwin 最终在 12 小时后终止了程序。 （我将 match_count 留在那里计算总数，请随意忽略该行）

#!/usr/bin/perl
use strict;
use warnings;

my $k = 15;
my $input = 'fasta.fasta';
my $output = 'text.txt';
my $match_count = 0;

#Open File
unless (open(FASTA, "<", $input)){
    die "Unable to open fasta file", $!;
    }

    #Unwraps the FASTA format file
    $/=">";
    #Separate header and sequence
    #Remove spaces
unless (open(OUTPUT, ">", $output)){
die "Unable to open file", $!;
}

    while (my $line = <FASTA>){
            my($header, @seq) = split(/\n/, $line);
                    my $sequence = join '', @seq;

    while (length($sequence) >= $k){
        $sequence =~ m/(.{$k})/;
        print OUTPUT "$1\n";
        $sequence = substr($sequence, 1, length($sequence)-1);
    }
}

我要找的结果是：

A total of 20938309 k-mers printed in the text file when I use the wc -l command.

提前致谢！

Answer 1

不知道为什么你没有得到你想要的结果。

我想我会按照您的问题描述发布我使用过的 2 个程序。

第一个仅计算我用于测试的文件中的 kmers (fasta_dat.txt)。它不会将它们打印出来，而只是检查有多少 kmer。

#!/usr/bin/perl
use strict;
use warnings;
use Bio::SeqIO;

my $in  = Bio::SeqIO->new( -file   => "fasta_dat.txt" ,
                           -format => 'fasta');

my $count_kmers;
my $k = 15;
while ( my $seq = $in->next_seq) {
    $count_kmers += $seq->length - $k + 1;
}

print $count_kmers;

__END__
C:\Old_Data\perlp>perl t9.pl
18657

您可以看到计数（在__END__ 令牌之后），18657。当我使用您的代码打印出来时，这个计数与 kmers 的计数一致。

#!/usr/bin/perl
use strict;
use warnings;
use 5.014;
use Devel::Size 'total_size';

my $k = 15;
my $input = 'fasta_dat.txt';
my $output = 'kmers.txt';
my $match_count = 0;

#Open File
unless (open(FASTA, "<", $input)){
    die "Unable to open fasta file", $!;
    }

    #Unwraps the FASTA format file
    $/=">";
    #Separate header and sequence
    #Remove spaces
unless (open(OUTPUT, ">", $output)){
    die "Unable to open file", $!;
}

<FASTA>; # discard 'first' 'empty' record

my %seen;
while (my $line = <FASTA>){
    chomp $line;
    my($header, @seq) = split(/\n/, $line);
    my $sequence = join '', @seq;

    for my $i (0 .. length($sequence) - $k) {
        my $kmer = substr($sequence, $i, $k);
        print OUTPUT $kmer, "\n" unless $seen{$kmer}++;
    }
}
print total_size(\%seen);

更新 我运行的测试显示，哈希大小的内存增加了大约 100 倍。在我的测试中，kmers 的数量约为 18500。这导致哈希大小为 1.8MB。

对于您的数据，kmers 为 22M，将导致哈希大小约为 2.2GB。不知道这会不会超出你的记忆容量。