使用Perl散列处理制表符分隔文件

Question

我有两个文件：

• file_1有三列（标记（SNP），染色体和位置）
• file_2有三列（染色体，peak_start，和peak_end）。

除SNP列外，所有列都是数字。

文件排列如屏幕截图所示。 file_1有几百个SNP作为行，而file_2有61个峰。每个峰都由一个peak_start和peak_end标记。可以有任何一个文件中的23条染色体，file_2每个染色体有几个峰。

我想查找file_1中SNP的位置是否落入file_2中的peak_start和peak_end，以确定每个匹配的染色体。如果是这样，我想显示哪个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引发的问题很好，尝试将这些纳入您的程序中。 （我从@David的帖子中借用了大部分代码。）

我不明白的一件事是，为什么加载散列值的峰值和位置，因为加载一个就足够了。由于每个染色体有多个记录，使用HoA。我的解决方案基于此。您可能需要更改列和他们的位置。

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循环，因为你期望找到一些SNP的第二不少。因此，通过你的%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，你应该使用3参数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：

我做了最后的循环，打印所有的东西，你指定（标签定界），但你没有指定的格式。如果你必须改变它。

#! /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来隐藏哈希散列的结构。这可以防止错误，例如认为开始窥视存储在START_PEEK而不是PEAK_START中。 Perl不会检测到这类错误的错误。因此，我更喜欢每当使用数组的数组或散列哈希时使用对象。

Answer 4

我用@tuxuday和@ David的代码来解决这个问题。这是最终的代码，做了我想要的。我不仅学到了很多，但我已经能够成功解决我的问题！荣誉家伙！

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;