-6

我想要一个简单的 perl 脚本，它可以帮助我在对齐序列中估算缺失的核苷酸：例如，我的 old_file 包含以下对齐序列：

seq1
ATGTC
seq2
ATGTC
seq3
ATNNC
seq4
NNGTN
seq5
CTCTN

所以我现在想推断文件中的所有 N，并获得一个新文件，其中所有 N 都是根据特定位置的多数核苷酸推断出来的。我的 new_file 应该是这样的：

seq1
ATGTC
seq2
ATGTC
seq3
ATGTC
seq4
ATGTC
seq5
CTCTC

具有用法的脚本：“impute_missing_data.pl old_file new_file”或任何其他方法将对我有所帮助。谢谢你。

4

3 回答 3

2

这似乎做了所需要的

use strict;
use warnings;

use Fcntl 'SEEK_SET';

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

my @counts;

while (<$fh>) {
  next if /[^ATGCN\s]/;
  my $i = 0;
  $counts[$i++]{$_}++ for /[ATGC]/g;
}

for my $maj (@counts) {
  ($maj) = sort { $maj->{$b} <=> $maj->{$a} } keys %$maj;
}

seek $fh, 0, SEEK_SET;

while (<$fh>) {
  s/N/$counts[pos]/eg unless /[^ATGCN\s]/;
  print;
}

输出

seq1
ATGTC
seq2
ATGTC
seq3
ATGTC
seq4
ATGTC
seq5
CTCTC
于 2012-09-06T16:44:05.393 回答
0 
use warnings;
use strict;
my (@data, $counts, @max);
#read in the file
while (<>) {
  chomp;
  next if /seq/;
  my @sings = split //; 
  for (my $i = 0; $i < @sings; $i++) {
    $counts->[$i]{$sings[$i]}++ if $sings[$i] ne 'N';
  }
  push (@data, \@sings);
}
# get most freq letters
foreach my $col (@$counts) {
  my ($max, $maxk) = (0, '');
  foreach my $cell (keys %$col) {
    if ($col->{$cell} > $max) {
      ($max, $maxk) = ($col->{$cell}, $cell);
    }   
  }
  push (@max, $maxk);
}
# substitute Ns with most freq letters
foreach (my $i = 0; $i < @data; $i++) {
  my $row = $data[$i];
  for (my $i = 0; $i < @$row; $i++) {
    if ($row->[$i] eq 'N') {
      $row->[$i] = $max[$i];
    }   
  }
  print "seq".($i+1)."\n".join("", @$row), "\n";
}
于 2012-09-06T16:03:31.083 回答
-1

这是我评论中的脚本以更易读的形式：

#!/usr/bin/perl
use strict;
my @stat;
while(<>) {
  print and next if /^seq/;
  chomp;
  my @seq = split //;
  for my $i (0..$#seq){
    my ($e, %s) = ($seq[$i], %{$stat[$i]}); # read-only aliases
    if($e=~/N/){
      my $substitution = [sort {$s{$a} <=> $s{$b}} keys %s]->[-1];
      $seq[$i] = $substitution;
      warn "substituted N with $substitution in col $i, count $s{$substitution}\n";
    } else {
      $stat[$i]->{$e}++;
    }
  }
  print @seq, "\n"';
}

要抑制undefined警告，请将它们静音（坏）或初始化统计信息：

for my $i (0..4) {
  for my $c (qw(A C G T)) {
    $stat[$i]->{$c} = 0;
}

或者

my @stat = map +{map {$_ => 0} qw(A C G T)} 0..4;
于 2012-09-06T16:22:36.740 回答