2

我对 Perl 和脚本非常陌生，但我需要一个代码来进行我的研究。我正在尝试计算存储在 multiFASTA 文件中的 DNA 序列中 11-mers 的频率。通过合并我找到的一些脚本，我写了这个：

#!/usr/bin/perl

$k = 11;  @bases = ('A','C','G','T');
@words = @bases; open FILE1, ">kmers.txt" or die $!;
for $i (1..$k-1)  {
   undef @newwords;
   foreach $w (@words)
   {
       foreach $b (@bases)
       {
          push (@newwords,$w.$b);
       }
   }
   undef @words;
   @words = @newwords;  
}
foreach $w (@words) {  
   print FILE1 "$w \n"; 
} 
close FILE1;   
my $input=$ARGV[0]; 
my $output=$ARGV[1];
open(IN,"<$input") || die ("Error opening $input $!"); 
open OUT, ">$output" or die $|; my $line = <IN>;  
print OUT $line; 
while ($line = <IN>) { 
   chomp $line; 
   if ($line=~m/^>/) { 
      print OUT  "\n",$line,"\n"; 
   } else { 
      print OUT $line; 
   } 
} 
print OUT "\n";

chomp $seq; chomp $k;
#obtain all distinct kmers open FILE2, ">out.txt" or die $!;

for $line (@lines) { 
   if ($line=~m/^>/) { next; } 
}
foreach($i=1; length($line) >= $k; $i++)    {   
   $line =~ m/(^.{$k})/;  
   $w{$1}{cnt}++;
   push @{$w{$1}{pos}}, $i;  
   $line= substr($seq, 1, length($line)-1);
   foreach $line (keys %kmers)    {
      print FILE2 "$kmers\n";
   }
   close FILE2; 
   close OUT;    
}

基本上，它读取文件，将所有序列行放在单独文件的一行中，记下所有 11mers 并创建一个“out.txt”文件，我希望他在其中存储具有 11-mer 频率的序列头。这是困难的部分（对我来说）：我如何告诉脚本编写序列标题以及每个序列的 11mer 频率？

4

3 回答 3

0

在对代码进行了一些修改之后，我产生了这个：

use strict;
use warnings;
my $in_file = $ARGV[0];
my $out_tvir = $ARGV[1];
my $k = $ARGV[2];
my %seq_hash; # key = seq_name, value = seq;
{
# redefine the record separator
local $/ = ">";
open IN, "<$in_file";
my $in_line = <IN>; # toss the first record
while ( $in_line = <IN> ) {
    chomp $in_line; # remove the ">" character in the end 
    my ( $seq_name, $seq ) = split( /\n/, $in_line, 2 );
    $seq =~ tr/ \t\n\r//d;    # Remove whitespace
    $seq_hash{$seq_name} = uc $seq;
}
close IN;
}

open OUT, ">$out_file";
open OUT2, ">$out_tvir";
foreach my $seq_name ( sort keys %seq_hash ) {
chomp $k;
%kmers = ();
while (length($seq_hash{$seq_name}) >= $k)
    {
    $seq_hash{$seq_name}=~ m/(^.{$k})/;
    $kmers{$1}++;
    $seq_hash{$seq_name}= substr($seq_hash{$seq_name}, 1,         length($seq_hash{$seq_name})-1);
    }
    $num_kmers = keys %kmers;
$px=();
$logpx=();
my $H=();
foreach $str (keys %kmers)
{
    my $px=$kmers{$str}/$num_kmers;
    $logpx=log($px);
    $H -= $px * log($px);
    if ($H <= 18) {print OUT2 ">$seq_name\t$H\n";}
}
}
close OUT;

...哪种工作，如果我省略最后一个“if（$ H ...”部分，只是让脚本通过列出与每个序列关联的所有 H 值来完成工作。我不知道为什么，尽管。

于 2013-02-01T16:47:19.287 回答
0

user2029917，您在未声明的变量方面遇到了一些问题，这会阻止脚本在use strict;打开的情况下运行；我做了一些修改并清理了一下。 

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

my $in_file = $ARGV[0];
my $out_tvir = $ARGV[1];
my $k = $ARGV[2];

my %seq_hash; # key = seq_name, value = seq;
{
    # redefine the record separator
    local $/ = ">";
    open IN, "<", $in_file or die "Can't open ${in_file}: $!";
    my $in_line = <IN>; # toss the first record
    while ( $in_line = <IN> ) {
        chomp $in_line; # remove the ">" character in the end 
        my ( $seq_name, $seq ) = split( /\n/, $in_line, 2 );
        $seq =~ tr/\t\n\r//d;    # Remove whitespace
        $seq_hash{$seq_name} = uc $seq;
    }
    close IN;
}

open OUT, ">", $out_tvir or die "Can't open ${out_tvir}: $!";
foreach my $seq_name ( sort keys %seq_hash ) {
    chomp $k;
    my %kmers;
    while (length($seq_hash{$seq_name}) >= $k) {
        $seq_hash{$seq_name}=~ m/(^.{$k})/;
        $kmers{$1}++;
        $seq_hash{$seq_name}= substr($seq_hash{$seq_name}, 1, length($seq_hash{$seq_name})-1);
    }
    my $num_kmers = keys %kmers;
    my $px;
    my $logpx;
    my $H;
    foreach my $str (keys %kmers) {
        my $px=$kmers{$str}/$num_kmers;
        $logpx=log($px);
        $H -= $px * log($px);
        if ($H <= 18) {print OUT ">$seq_name\t$H\n";}
    }
}

close OUT;

exit;

它现在应该运行，但我不确定这个脚本是否会产生您想要的输出。例如，对于给定的 k-mer，它将打印它出现的每个 FASTA 条目的 H' 值（无论 FASTA 条目如何，它始终是相同的值，因为它是用总出现次数和总数计算的k-mers)。目前，它不打印所指的是哪个 k-mer。这可以通过将最后一位更改为 来解决print OUT ">$seq_name\t$str\t$H\n";，但我不确定这是否是您所追求的行为。如果您可以提供有关所需输出的更多详细信息，我们可能会提供更多帮助。

于 2013-02-18T15:32:29.767 回答
0

你不需要创建一个中间文件来强制每个序列是一行，你的循环计数 11mers 可以更简单，我不明白你为什么要创建带有所有可能的 11mers 的 kmers.txt，因为你没有'不要使用它。此外，还有许多语法错误和未使用的变量。而且您不需要 substr 的第三个参数。如果您将其关闭，则默认值将结束。

除了这些问题，请注意可能存在 4^11 个可能的 11mer，这几乎是 420 万个可能性。您的定义将是巨大的（取决于您正在分析的序列的长度）。我猜一个典型的基因会超过一千个 11mer，除非你正在分析重复序列。您可能会考虑在定义中仅包含任意数量的最丰富的 11mer（除非您计划以编程方式处理输出 - 但即便如此，有这么长的行可能会有问题）。

您提交的答案有一些额外的意图，这不在您的问题中，但暂时搁置一旁，这就是我将如何编写脚本以在您的序列中包含前 5 个 11mer 频率（任意选择第 5 位的任何关系）。我不会解决其他人建议您的正确编码实践问题 - 但您应该注意这些建议。

my $input=$ARGV[0]; 
my $output=$ARGV[1];

my $defline = '';
my $seq = '';
my $topkmers = '';

open(INPUT,$input);
open(OUTPUT,">$output");
select(OUTPUT);

while(<INPUT>)
  {
    chomp;
    if(/^>/)
      {
        if($seq ne '')
          {
            $topkmers = getTopKMers($seq,11,5);
            print("$defline $topkmers\n$seq\n");
          }
        $defline = $_;
        $seq = '';
      }
    else
      {$seq .= $_}
  }
#Take care of the last record
if($seq ne '')
  {
    $topkmers = getKMers($seq);
    print("$defline $topkmers\n$seq\n");
  }

close(INPUT);
close(OUTPUT);

sub getTopKMers
  {
    my $seq = uc($_[0]);
    my $size = $_[1];
    my $top = $_[2] - 1;   #Submit a 0 to get all kmers
    my $hash = {};

    #Create the abundance hash
    for(my $p = 0;$p < (length($seq) - $size);$p++)
      {push(@{$hash->{substr($seq,$p,$size)}},$p}

    #Sort by abundance
    my @sorted = sort {scalar(@{$hash->{$b}}) <=> scalar(@{$hash->{$a}})} keys(%$hash);

    #Get the top few most abundant kmers
    my @toplist = $top > -1 ? @sorted[0..$top] : @sorted;

    #Creates a string like "ATGCATGCCAA[20]=1,2,... CGTAGCTCTAG[18]=6,23,..."
    my $str = join(' ',
                   map {
                        "$_\[" .
                        scalar(@{$hash->{$_}}}) .
                        "]=" .
                        join(',',@{$hash->{$_}})
                       } @toplist);

    return($str);
  }

这可以消除排序并合并几个步骤以使其更高效，但是为了更易于阅读代码，还有一些话要说。

注意：我没有运行此代码，所以请原谅我忽略的任何错误。

于 2016-03-31T21:18:56.710 回答