☰

如何比较文件 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处的数组中。NR当读取第二个文件时，将数组中的值a和b给定索引FNR（文件的记录号）处的值与字段$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]));
}

这将两个文件名参数分别打开为文件句柄F1和。F2如果任一文件无法打开，它将退出并显示错误消息。

然后，虽然它们都没有到达 EOF（文件结尾），但它：

  • 一次从每个文件句柄读取一行到单独的数组中（@a对于 F1 和@bF2）。

    该函数从 array 最后一个元素中chomp()删除尾随换行符 ( ) ，以防止换行符出现在数组连接的任何输出行的中间。\n@a

  • 如果数组满足您的条件（$a[1]>=$b[3] 和 $a[2]<=$b[4]），它会将两个数组打印为一个输出行，并使用制表符作为字段分隔符连接。

将其另存为，例如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

相关内容