2

我有两个文件：

  • file_1 有三列（标记（SNP）、染色体和位置）
  • file_2 有三列（Chromosome、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
4

4 回答 4

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 (my $peaks_fh, "<", "peaks.txt") 
   or die "Couldn't open peaks.txt";

这样，如果您的文件名恰好以|or开头>，它仍然可以使用。使用标量变量（以 a 开头的变量$）可以更轻松地在函数之间传递文件句柄。

无论如何，只是为了确保我正确理解你：你说“我更喜欢......在染色体是关键的地方使用散列。 ”

现在，我有 23 对染色体，但每条染色体上可能都有数千个 SNP。如果以这种方式按染色体键控，则每个染色体只能存储一个 SNP。这是你想要的吗？我注意到您的数据显示所有相同的染色体。这意味着您不能按染色体键控。我暂时忽略了这一点，并使用我自己的数据。

我还注意到您所说的文件包含的内容以及您的程序如何使用它们有所不同：

您说：“文件 1 有 3 列（SNP、染色体和位置） ”，但您的代码是：

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

我假设是染色体、位置和 SNP。文件以哪种方式排列？

你必须明确你的要求。

无论如何，这是以制表符分隔格式打印的经过测试的版本。这是一种更现代的 Perl 格式。请注意，我只有一个按染色体的哈希（如您指定的那样）。我先读了peaks.txtin。如果我在我的位置文件中发现我的文件中不存在的染色体peaks.txt，我会忽略它。否则，我将为POSITION和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 存储在START_PEEK而不是PEAK_START. Perl 不会检测到这些类型的误键错误。因此，每当我处理数组数组或哈希值时，我更喜欢使用对象。
于 2012-05-14T03:56:39.943 回答
1

您只需要一个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'};
        }
    }
}

注意：我没有测试过代码。

查看您添加的屏幕截图，这是行不通的。

于 2012-05-13T23:27:54.140 回答
0

@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;
于 2012-05-14T07:33:02.233 回答
0

我使用@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;
于 2012-05-15T16:33:12.953 回答