-6

I want a simple perl script that can help me impute missing nucleotides in aligned sequences: As an example, my old_file contains the following aligned sequences:

seq1
ATGTC
seq2
ATGTC
seq3
ATNNC
seq4
NNGTN
seq5
CTCTN

So I now want to infer all Ns in the file and get a new file with all the Ns inferred based on the majority nucleotide at a particular position. My new_file should look like this:

seq1
ATGTC
seq2
ATGTC
seq3
ATGTC
seq4
ATGTC
seq5
CTCTC

A script with usage: "impute_missing_data.pl old_file new_file" or any other approach will be helpful to me. Thank you.

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 に答える