Использование Perl хэши для обработки разделителями табуляции файлов

Question

У меня есть два файла:

• file_1 имеет три колонки (Marker (SNP), Хромосома, и позиция)
• file_2 имеет три колонки (хромосома, peak_start, и peak_end).

Все столбцы являются числовыми, за исключением столбца SNP.

Файлы упорядочены, как показано на скриншотах. file_1 имеет несколько сотен SNP в виде строк, а файл_2 имеет 61 пик. Каждый пик отмечен значком peak_start и peak_end. В любом файле может быть любая из 23 хромосом, а файл_2 имеет несколько пиков на хромосому.

Я хочу найти, если позиция SNP в файле_1 попадает в пик_старт и пик_ен в файле_2 для каждой соответствующей хромосомы. Если это так, я хочу показать, какой SNP падает в том пике (желательно записать вывод в файл с разделителями табуляции).

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

Вот пример моего кода. Он предназначен только для иллюстрации моего вопроса и до сих пор ничего не делает, поэтому подумайте об этом как о «псевдокоде».

#!usr/bin/perl 

use strict; 
use warnings; 

my (%peaks, %X81_05); 
my @array; 

# Open file or die 

unless (open (FIRST_SAMPLE, "X81_05.txt")) { 
    die "Could not open X81_05.txt"; 
} 

# Split the tab-delimited file into respective fields 

while (<FIRST_SAMPLE>) { 

    chomp $_; 
    next if (m/Chromosome/); # Skip the header 

    @array = split("\t", $_); 
    ($chr1, $pos, $sample) = @array; 

    $X81_05{'$array[0]'} = (
     'position' =>'$array[1]' 
    ) 
} 

close (FIRST_SAMPLE); 

# Open file using file handle 
unless (open (PEAKS, "peaks.txt")) { 
    die "could not open peaks.txt"; 
} 

my ($chr, $peak_start, $peak_end); 

while (<PEAKS>) { 
    chomp $_; 

    next if (m/Chromosome/); # Skip header 
    ($chr, $peak_start, $peak_end) = split(/\t/); 
    $peaks{$chr}{'peak_start'} = $peak_start; 
    $peaks{$chr}{'peak_end'} = $peak_end; 
} 

close (PEAKS); 

for my $chr1 (keys %X81_05) { 
    my $val = $X81_05{$chr1}{'position'}; 

    for my $chr (keys %peaks) { 
     my $min = $peaks{$chr}{'peak_start'}; 

     my $max = $peaks{$chr}{'peak_end'}; 

     if (($val > $min) and ($val < $max)) { 
      #print $val, " ", "lies between"," ", $min, " ", "and", " ", $max, "\n"; 
     } 
     else { 
       #print $val, " ", "does not lie between"," ", $min, " ", "and", " ", $max, "\n"; 
     } 
    } 
} 

Более удивительный код:

1. http://i.stack.imgur.com/fzwRQ.png
2. http://i.stack.imgur.com/2ryyI.png

Answer 1

Очки, поднятые @David, хороши; попробуйте включить их в свои программы. (Я взял большую часть кода из сообщения Дэвида.)

Одна вещь, которую я не понимал, - это то, зачем загружать как пиковые значения, так и позицию в хеш, так как загрузка будет достаточной. Поскольку каждая хромосома имеет более одной записи, используйте HoA. Мое решение основано на этом. Вам может потребоваться изменить cols и их позиции.

use strict; 
use warnings; 

our $Sep = "\t"; 
open (my $peak_fh, "<", "data/file2"); 
my %chromosome_hash; 

while (my $line = <$peak_fh>) { 
    chomp $line; 
    next if $line =~ /Chromosome/; #Skip Header 
    my ($chromosome) = (split($Sep, $line))[0]; 
    push @{$chromosome_hash{$chromosome}}, $line; # Store the line(s) indexed by chromo 
} 
close $peak_fh; 

open (my $position_fh, "<", "data/file1"); 

while (my $line = <$position_fh>) { 
    chomp $line; 
    my ($chromosome, $snp, $position) = split ($Sep, $line); 
    next unless exists $chromosome_hash{$chromosome}; 

    foreach my $peak_line (@{$chromosome_hash{$chromosome}}) { 
     my ($start,$end) = (split($Sep, $line))[1,2]; 

     if ($position >= $start and $position <= $end) { 
      print "MATCH REQUIRED-DETAILS...$line-$peak_line\n"; 
     } 
     else { 
      print "NO MATCH REQUIRED-DETAILS...$line-$peak_line\n"; 
     } 
    } 
} 
close $position_fh; 

Answer 2

Вам нужно только один for петлю, потому что вы ждете, чтобы найти некоторые из ОНП во второй партии. Следовательно, проведите через свой хэш-код %X81_05 и проверьте, соответствует ли он одному из %peak. Что-то вроде:

for my $chr1 (keys %X81_05) 
{ 
    if (defined $peaks{$chr1}) 
    { 
     if ( $X81_05{$chr1}{'position'} > $peaks{$chr1}{'peak_start'} 
      && $X81_05{$chr1}{'position'} < $peaks{$chr1}{'peak_end'}) 
     { 
      print YOUROUTPUTFILEHANDLE $chr1 . "\t" 
       . $peaks{$chr1}{'peak_start'} . "\t" 
       . $peaks{$chr1}{'peak_end'}; 
     } 
     else 
     { 
      print YOUROUTPUTFILEHANDLE $chr1 
       . "\tDoes not fall between " 
       . $peaks{$chr1}{'peak_start'} . " and " 
       . $peaks{$chr1}{'peak_end'}; 
     } 
    } 
} 

Примечание: Я не проверял код.

Глядя на скриншоты, которые вы добавили, это не сработает.

Answer 3

Несколько программных намеков на Perl:

Вы можете сделать это:

open (PEAKS, "peaks.txt") 
    or die "Couldn't open peaks.txt"; 

Вместо этого:

unless (open (PEAKS, "peaks.txt")) { 
    die "could not open peaks.txt"; 
} 

Это больше стандартного Perl, и это немного легче читать.

Говоря о стандартной Perl, вы должны использовать три аргумента open форму, и использовать скаляры для дескрипторов файлов:

open (my $peaks_fh, "<", "peaks.txt") 
    or die "Couldn't open peaks.txt"; 

Таким образом, если имя вашего файла просто происходит, чтобы начать с | или >, его будет по-прежнему работать.Использование переменных скаляров (переменные, начинающиеся с $) упрощают передачу файлов между функциями.

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

Теперь у меня есть 23 пары хромосом , но каждая из этих хромосом может иметь на нем тысячи SNP. Если вы используете хромосому таким образом, вы можете хранить только один SNP на хромосому. Это то, что вы хотите? Я замечаю, что ваши данные показывают все те же хромосомы. Это означает, что вы не можете управлять хромосомой. Я игнорирую это сейчас и использую свои собственные данные.

Я также заметил разницу в том, что вы сказали, что файлы содержали, и как ваша программа использует их:

Вы сказали: «Файл 1 имеет 3 колонки (SNP, Хромосома, и положение)» , но ваш код:

($chr1, $pos, $sample) = @array; 

Который я предполагаю, это хромосома, положение и SNP. Каким образом находится файл?

Вы должны прояснить, что именно вы просите.

В любом случае, это проверенная версия, которая выводится в формате с разделителями табуляции. Это немного более современный формат Perl. Обратите внимание, что у меня есть только один хеш по хромосоме (как вы указали). Сначала я прочитал peaks.txt. Если я найду в своем файле позиции хромосому, которая не существует в моем файле peaks.txt, я просто игнорирую ее. В противном случае, я добавлю в дополнительных хэшей для ПОЗИЦИИ и SNP:

Я делаю окончательный цикл, который выводит все из (вкладка delimitated), как было указано, но вы не указали формат. Измените его, если вам нужно.

#! /usr/bin/env perl 

use strict; 
use warnings; 
use feature qw(say); 
use autodie;  #No need to check for file open failure 
use constant { 
    PEAKS_FILE  => "peak.txt", 
    POSITION_FILE => "X81_05.txt", 
}; 

open (my $peak_fh, "<", PEAKS_FILE); 
my %chromosome_hash; 

while (my $line = <$peak_fh>) { 
    chomp $line; 
    next if $line =~ /Chromosome/; #Skip Header 
    my ($chromosome, $peak_start, $peak_end) = split ("\t", $line); 
    $chromosome_hash{$chromosome}->{PEAK_START} = $peak_start; 
    $chromosome_hash{$chromosome}->{PEAK_END} = $peak_end; 
} 
close $peak_fh; 

open (my $position_fh, "<", POSITION_FILE); 

while (my $line = <$position_fh>) { 
    chomp $line; 
    my ($chromosome, $position, $snp) = split ("\t", $line); 
    next unless exists $chromosome_hash{$chromosome}; 

    if ($position >= $chromosome_hash{$chromosome}->{PEAK_START} 
      and $position <= $chromosome_hash{$chromosome}->{PEAK_END}) { 
     $chromosome_hash{$chromosome}->{SNP} = $snp; 
     $chromosome_hash{$chromosome}->{POSITION} = $position; 
    } 
} 
close $position_fh; 

# 
# Now Print 
# 

say join ("\t", qw(Chromosome, SNP, POSITION, PEAK-START, PEAK-END)); 
foreach my $chromosome (sort keys %chromosome_hash) { 
    next unless exists $chromosome_hash{$chromosome}->{SNP}; 
    say join ("\t", 
     $chromosome, 
     $chromosome_hash{$chromosome}->{SNP}, 
     $chromosome_hash{$chromosome}->{POSITION}, 
     $chromosome_hash{$chromosome}->{PEAK_START}, 
     $chromosome_hash{$chromosome}->{PEAK_END}, 
    ); 
} 

Несколько вещей:

• Оставлять пространства вокруг круглых скобок с обеих сторон. Это облегчает чтение.
• Я использую круглые скобки, если другие нет. Нынешний стиль не должен их использовать, если вам не нужно. Я использую их для всех функций, которые принимают более одного аргумента. Например, я мог бы сказать open my $peak_fh, "<", PEAKS_FILE;, но я думаю, что параметры начинают теряться, когда у вас есть три параметра для функции.
• Извещение Я использую use autodie;. Это приводит к завершению работы программы, если она не может открыть файл. Вот почему мне даже не нужно проверять, открыт ли файл.
• Я бы предпочел использовать объектно-ориентированный Perl, чтобы скрыть структуру хэша хэшей. Это предотвращает такие ошибки, как мысль о том, что начальный peek хранится в START_PEEK, а не PEAK_START. Perl не будет обнаруживать эти ошибки miskeyed. Поэтому я предпочитаю использовать объекты всякий раз, когда я делаю массивы массивов или хешей хешей.

Answer 4

Я использовал @tuxuday и код Дэвида для решения этой проблемы. Вот последний код, который сделал то, что я хотел. Я не только многому научился, но и смог успешно решить мою проблему! Кудос, ребята!

use strict; 
use warnings; 
use feature qw(say); 

# Read in peaks and sample files from command line 
my $usage = "Usage: $0 <peaks_file> <sample_file>"; 
my $peaks = shift @ARGV or die "$usage \n"; 
my $sample = shift @ARGV or die "$usage \n"; 

our $Sep = "\t"; 
open (my $peak_fh, "<", "$peaks"); 
my %chromosome_hash; 

while (my $line = <$peak_fh>) { 
    chomp $line; 
    next if $line =~ /Chromosome/; #Skip Header 
    my ($chromosome) = (split($Sep, $line))[0]; 

    push @{$chromosome_hash{$chromosome}}, $line; # Store the line(s) indexed by chromosome 
} 
close $peak_fh; 

open (my $position_fh, "<", "$sample"); 

while (my $line = <$position_fh>) { 
    chomp $line; 
    next if $line =~ /Marker/; #Skip Header 
    my ($snp, $chromosome, $position) = split ($Sep, $line); 

    # Check if chromosome in peaks_file matches chromosome in sample_file 
    next unless exists $chromosome_hash{$chromosome}; 

    foreach my $peak_line (@{$chromosome_hash{$chromosome}}) { 

     my ($start,$end,$peak_no) = (split($Sep, $peak_line))[1,2,3]; 

     if ($position >= $start and $position <= $end) { 

      # Print output 
      say join ("\t", 
       $snp, 
       $chromosome, 
       $position, 
       $start, 
       $end, 
       $peak_no, 
      ); 
     } 
     else { 
      next; # Go to next chromosome 
     } 
    } 
} 
close $position_fh;