1

私はこのような巨大なタブ区切りファイルを持っています：

contig04733 contig00012 77
contig00546 contig01344 12
contig08943 contig00001 14
contig00765contig0312588
など

そして、私は次のようなこれらのコンティグペアのサブセットのみを含む別のタブ区切りファイルを持っています：

contig04733
contig00012contig08943contig00001
など

2番目にリストされている行に対応する最初のファイルの行を新しいファイルに抽出したいと思います。この特定のデータセットでは、各ペアのどちらの方向が両方のファイルで同じである必要があると思います。しかし、次のように言うかどうかも知りたいです：

file1 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;

を使用すると、プログラムはコマンドラインで指定された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 FAQの例の多くは、%seen「これを見た」のように、という名前のハッシュを使用しています。</ p>

を使用して2番目のコマンドライン引数を削除するとpop、マスターファイルのみが残り、プログラムは。を使用してすべての入力行を簡単に読み取ることができますwhile (<>) { ... }。値が設定されると%subset、コードは各行をフィールドに分割し、表示済みとマークされていない行をスキップします。このフィルターを通過するものはすべて、標準出力に出力されます。

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

例えば：

$ cat file1
contig04733 contig00012 77
contig00546 contig01344 12
contig08943 contig00001 14
contig00765 contig03125 88

$ cat file2
contig04733 contig00012
contig00001 contig08943

$ perl extract-subset file1 file2
contig04733 contig00012 77
contig08943 contig00001 14

選択したサブセットを含む新しい出力を作成するには、次のように標準出力をリダイレクトします。

$ perl extract-subset file1 file2> my-subset
于 2013-03-08T16:48:05.510 に答える
0

2つのキーに基づくハッシュのハッシュを使用してこれを試してください（分割後）

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 に答える