1

我有一个巨大的制表符分隔文件，如下所示：

contig04733 contig00012 77
contig00546 contig01344 12
contig08943 contig00001 14
contig00765 contig03125 88
等

我有一个单独的制表符分隔文件，其中只有这些重叠群对的一个子集，如下所示：

contig04733 contig00012
contig08943 contig00001
等

我想将第一个文件中与第二个文件中列出的行对应的行提取到一个新文件中。在这个特定的数据集中，我认为在两个文件中每对的哪一种方式应该是相同的。但也想知道是否说：

文件1 contig08943 contig00001 14

但在file2中

contig00001 contig08943

我仍然想要这种组合，是否也可以为此编写脚本？

我的代码如下。

use strict;
use warnings;

#open the contig pairs list
open (PAIRS, "$ARGV[0]") or die "Error opening the input file with contig pairs";

#hash to store contig IDs - I think?!
my $pairs;

#read through the pairs list and read into memory?
while(<PAIRS>){
    chomp $_; #get rid of ending whitepace
    $pairs->{$_} = 1;
}
close(PAIRS);

#open data file
open(DATA, "$ARGV[1]") or die "Error opening the sequence pairs file\n";
while(<DATA>){
    chomp $_;
    my ($contigs, $identity) = split("\t", $_);
    if (defined $pairs->{$contigs}) {
        print STDOUT "$_\n";
    }
}
close(DATA);
4

2 回答 2

0

将下面没有运行注释的代码拼凑在一起，以获得一个工作程序。我们从典型的正面内容开始，perl如果您犯了常见错误，这些内容会指示给您有用的警告。

#! /usr/bin/env perl

use strict;
use warnings;

在必要时向用户展示如何正确调用您的程序总是一个不错的选择。

die "Usage: $0 master subset\n" unless @ARGV == 2;

使用read_subset，程序读取命令行中命名的第二个文件。因为您的问题表明您不关心订单，例如，

contig00001 contig08943

相当于

contig08943 contig00001

代码同时增加$subset{$p1}{$p2}和$subset{$p2}{$p1}。

sub read_subset {
  my($path) = @_;

  my %subset;
  open my $fh, "<", $path or die "$0: open $path: $!";
  while (<$fh>) {
    chomp;
    my($p1,$p2) = split /\t/;
    ++$subset{$p1}{$p2};
    ++$subset{$p2}{$p1};
  }

  %subset;
}

在 Perl 程序中，使用散列来标记您的程序观察到的事件是非常常见的。事实上，Perl 常见问题解答中的许多示例都使用了名为 的哈希%seen，如“我见过这个”。</p>

通过用 删除第二个命令行参数pop，只留下主文件，让程序使用 . 轻松读取所有输入行while (<>) { ... }。填充后，代码将%subset每一行拆分为字段并跳过任何未标记为已看到的行。通过此过滤器的所有内容都打印在标准输出上。

my %subset = read_subset pop @ARGV;
while (<>) {
  my($f1,$f2) = split /\t/;
  next unless $subset{$f1}{$f2};
  print;
}

例如：

$猫文件1
contig04733 contig00012 77
contig00546 contig01344 12
contig08943 contig00001 14
contig00765 contig03125 88

$猫文件2
contig04733 contig00012
contig00001 contig08943

$ perl 提取子集文件 1 文件 2
contig04733 contig00012 77
contig08943 contig00001 14

要创建包含所选子集的新输出，请将标准输出重定向为

$ perl 提取子集文件 1 文件 2 >我的子集
于 2013-03-08T16:48:05.510 回答
0

尝试使用基于两个键的散列哈希（拆分后）

use strict;
use warnings;

#open the contig pairs list
open (PAIRS, "$ARGV[0]") or die "Error opening the input file with contig pairs";

#hash to store contig IDs - I think?!
#my $pairs;

#read through the pairs list and read into memory?
my %all_configs;
while(<PAIRS>){
    chomp $_; #get rid of ending whitepace
    my @parts = split("\t", $_); #split into ['contig04733', 'contig00012', 77]
    #store the entire row as a hash of hashes
    $all_configs{$parts[0]}{$parts[1]} = $_;
    #$pairs->{$_} = 1; 
}
close(PAIRS);

#open data file
open(DATA, "$ARGV[1]") or die "Error opening the sequence pairs file\n";
while(<DATA>){
    chomp $_;
    my ($contigs, $identity) = split("\t", $_);
    #see if we find a row, one way, or the other
    my $found_row = $all_configs{$contigs}{$identity} 
        || $all_configs{$identity}{$contigs};
    #if found, the split, and handle each part    
    if ($found_row) {
        my @parts = split("\t", $found_row);
        #same sequence as in first file
        my $modified_row  = $parts[0]."\t".$parts[1].(sqrt($parts[2]/100));
        #if you want to keep the same sequence as found in second file
        my $modified_row  = $contigs."\t".$identity.(sqrt($parts[2]/100));

        print STDOUT $found_row."\n"; #or
        print STDOUT $modified_row."\n";
    }
}
close(DATA);
于 2013-03-08T16:48:09.603 回答