Блог

Я использую Perl для генерации списка уникальных экзонов (которые являются единицами генов).

Я сгенерировал файл в этом формате (с сотнями тысяч строк):

chr1 1000 2000 gene1

chr1 3000 4000 gene2

chr1 5000 6000 gene3

chr1 1000 2000 gene4

Позиция 1 — это хромосома, позиция 2 — начальная координата экзона, позиция 3 — конечная координата экзона, а позиция 4 — это название гена.

Поскольку гены часто состоят из разного расположения экзонов, у вас есть один и тот же экзон в нескольких генах (см. Первый и четвертый наборы). Я хочу удалить эти «дубликаты», то есть удалить gene1 или gene4 (не важно, какой из них будет удален).

Я часами бился головой о стену, пытаясь выполнить то, что (я думаю) является простой задачей. Кто-нибудь может указать мне правильное направление? Я знаю, что люди часто используют хэши для удаления повторяющихся элементов, но это не совсем дубликаты (поскольку названия генов разные). Важно, чтобы я также не потерял название гена. В противном случае это было бы проще.

Вот совершенно нефункциональный цикл, который я пробовал. В массиве «exons» каждая строка хранится как скаляр, отсюда и подпрограмма. Не смейтесь. Я знаю, что это не работает, но, по крайней мере, вы можете видеть (я надеюсь), что я пытаюсь сделать:

 for (my $i = 0; $i < scalar @exons; $i  ) {
my @temp_line = line_splitter($exons[$i]);                      # runs subroutine turning scalar into array
for (my $j = 0; $j < scalar @exons_dup; $j  ) {
    my @inner_temp_line = line_splitter($exons_dup[$j]);        # runs subroutine turning scalar into array
    unless (($temp_line[1] == $inner_temp_line[1]) amp;amp;           # this loop ensures that the the loop
            ($temp_line[3] eq $inner_temp_line[3])) {           # below skips the identical lines
                if (($temp_line[1] == $inner_temp_line[1]) amp;amp;   # if the coordinates are the same
                    ($temp_line[2] == $inner_temp_line[2])) {   # between the comparisons
                        splice(@exons, $i, 1);                  # delete the first one
                    }
            }
}
  

}

 my @exons = (
    'chr1 1000 2000 gene1',
    'chr1 3000 4000 gene2',
    'chr1 5000 6000 gene3',
    'chr1 1000 2000 gene4'
);

my %unique_exons = map { 
    my ($chro, $scoor, $ecoor, $gene) = (split(/s /, $_));
    "$chro $scoor $ecoor" => $gene
} @exons;

print "$_ $unique_exons{$_} n" for keys %unique_exons;
  

Это придаст вам уникальность, и будет включено последнее название гена. Это приводит к:

 chr1 1000 2000 gene4 
chr1 5000 6000 gene3 
chr1 3000 4000 gene2
  

Вы можете использовать хэш для быстрой дедупликации, но вам нужен способ объединить части, которые вы хотите использовать для обнаружения дубликатов, в одну строку.

 sub extract_dup_check_string {
    my $exon = shift;
    my @parts = line_splitter($exon);
    # modify to suit:
    my $dup_check_string = join( ';', @parts[0..2] );
    return $dup_check_string;
}

my %seen;
@deduped_exons = grep !$seen{ extract_dup_check_string($_) }  , @exons;
  

Вы можете использовать хэш для отслеживания дубликатов, которые вы уже видели, а затем пропустить их. В этом примере предполагается, что поля во входном файле разделены пробелами:

 #!/usr/bin/env perl                                                                                                                                                                                                                   

use strict;
use warnings;

my %seen;
while (my $line = <>) {

  my($chromosome, $exon_start, $exon_end, $gene) = split /s /, $line;
  my $key = join ':', $chromosome, $exon_start, $exon_end;

  if ($seen{$key}) {
    next;
  }
  else {
    $seen{$key}  ;
    print $line;
  }

}
  

Настолько просто, насколько это возможно. Я старался использовать как можно меньше магии.

 my %exoms = ();
my $input;

open( $input, '<', "lines.in" ) or die $!;

while( <$input> )
{
    if( $_ =~ /^(w s ){3}(w )$/ ) #ignore lines that are not in expected format
    {
        my @splits = split( /s /, $_ ); #split line in $_ on multiple spaces

        my $key = $splits[1] . '_' . $splits[2];

        if( !exists( $exoms{$key} ) )
        {
            #could output or write to a new file here, probably output to a file
            #for large sets.
            $exoms{$key} = @splits; 
        }
    }

}


#demo to show what was parsed from demo input
while( my ($key, $value) = each(%exoms) )
{
    my @splits = @{$value};
    foreach my $position (@splits)
    {
        print( "$position " );
    }
    print( "n" );

}