☰

ファイル1の列2と3とファイル2の列4と5を比較する方法

tag-icon • tag-icon text-processing awk python bioinformatics
ファイル1の列2と3とファイル2の列4と5を比較する方法

タブ区切りのファイル1があります。

NC_025345       4569   4950   KX838946.2      
NC_025345       16546   17066   KJ641660.1      
NC_025345       11996   12085   KX932454.2

ファイル2：

NC_025345.1     RefSeq  gene    5690    7513    .       +       .       ID=gene-NZ82_gp4;Dbxref=GeneID:20964334;Name=NZ82_gp4;gbkey=Gene;gene_biotype=protein_coding;locus_tag=NZ82_gp4
NC_025345.1     RefSeq  gene    8016    10046   .       +       .       ID=gene-NZ82_gp5;Dbxref=GeneID:20964335;Name=NZ82_gp5;gbkey=Gene;gene_biotype=protein_coding;locus_tag=NZ82_gp5
NC_025345.1     RefSeq  gene    10337   16933   .       +       .       ID=gene-NZ82_gp6;Dbxref=GeneID:20964336;Name=NZ82_gp6;gbkey=Gene;gene_biotype=protein_coding;locus_tag=NZ82_gp6

ファイル1の列2と3をファイル2の列4と5と比較したいと思います。ファイル 1 の列 2 と列 3 がファイル 2 の列 4 と列 5 の間に重なっているか、該当する場合は、次の出力を使用して、ファイル 1 とファイル 2 の行全体を新しいファイルにマージしたいと思います。

NC_025345       11996   12085   KX932454.2     NC_025345.1     RefSeq  gene    10337   16933   .       +       .       ID=gene-NZ82_gp6;Dbxref=GeneID:20964336;Name=NZ82_gp6;gbkey=Gene;gene_biotype=protein_coding;locus_tag=NZ82_gp6

答え1

これは最も効率的なawkスクリプトではないかもしれません。

awk '{
  if (NR==FNR) {
    l[NR]=$0
    a[NR]=$2
    b[NR]=$3
  }
  else if (a[FNR]>=$4 && b[FNR]<=$5) {
    print l[FNR],$0
  }
}' file1 file2 > newfile

最初のファイル（）を読み取るときは、行とフィールド全体を配列NR==FNRの$0インデックス$2（レコード番号）に保存します$3。NR2番目のファイルを読み取ると、配列の値とa指定bされたインデックスFNR（ファイルのレコード番号）の値がフィールドsumと$4比較されます$5。

配列値が範囲内にある場合は、前の行と現在の行を印刷します。出力は新しいファイルに書き込まれますnewfile。

答え2

次のように実行できます。まず、ファイルを一緒に貼り付けてから、awkを実行して必要な行をフィルタリングします。

$ nf1=$(<file1 awk '{print NF;exit}') 

$ paste file1 file2 | awk -vnf1="$nf1" '$(4+nf1)<=$2 && $3<=$(5+nf1)'

答え3

#!/usr/bin/perl

use strict;

my $f1=shift;
open(F1,"<",$f1) || die "Couldn't open $f1: $!\n";

my $f2=shift;
open(F2,"<",$f2) || die "Couldn't open $f2: $!\n";

until (eof(F1) || eof(F2)) {
  my @a = split /\t/,<F1>;
  my @b = split /\t/,<F2>;
  chomp($a[@a-1]);

  # note: perl arrays start from 0, not 1
  print join("\t",@a, @b) if (($a[1] >= $b[3]) && ($a[2] <= $b[4]));
}

これにより、2つのファイル名引数がファイルハンドルF1とで開きます。F2どちらのファイルも開くことができない場合は、エラーメッセージで終了します。

その後、いずれもEOF（ファイルの終わり）に達しませんが、次のことを行います。

  • 各ファイルハンドルから一度に1行ずつ別々の配列（@aF1と@bF2の場合）に読み込みます。

    この関数は、配列の最後の要素から末尾の改行文字（）を削除しchomp()て、配列接続の出力行の途中に改行が表示されるのを防ぎます。\n@a

  • 配列が条件 ($a[1]>=$b[3] および $a[2]<=$b[4]) を満たす場合、タブ文字をフィールド区切り文字として使用して、2 つの配列を 1 つの出力行で印刷します。 。

たとえば、別の名前で保存してibk.pl実行可能にしたら、chmod +x ibk.pl次のように実行します。

$ ./ibk.pl file1 file2 
NC_025345       11996   12085   KX932454.2      NC_025345.1     RefSeq  gene    10337   16933   .       +       .       ID=gene-NZ82_gp6;Dbxref=GeneID:20964336;Name=NZ82_gp6;gbkey=Gene;gene_biotype=protein_coding;locus_tag=NZ82_gp6

関連情報