Bash/Awk 比较两个文件，在坐标之间打印值，否则打印 0答案

【问题标题】：Bash/Awk Compare two files, print value when it is between coordinates else print 0Bash/Awk 比较两个文件，在坐标之间打印值，否则打印 0
【发布时间】：2021-08-09 22:27:55
【问题描述】：

我有两个文件。如果“染色体”列在两个文件之间匹配，并且 File1 的位置在 File2 的 Start_position 和 End_position 之间，我想关联两个 cell_frac 值。如果 File2 中不存在 Gene（染色体 + 位置），我希望两个 cell_frac 值都等于 0。

文件 1：

Hugo_Symbol Chromosome  Position
Gene1       1           111111
Gene2       1           222222
Gene3       2           333333
Gene4       2           333337

文件2：

Chromosome  Start_Position  End_Position    cell_frac_A1    cell_frac_A2
1           222220          222230          0.12            0.01
2           333330          333340          0.03            0.25
3           444440          444450          0.01            0.01

期望的输出：

Hugo_Symbol Chromosome  Position    cell_frac_A1    cell_frac_A2
Gene1       1           111111      0               0
Gene2       1           222222      0.12            0.01
Gene3       2           333333      0.03            0.25
Gene4       2           333337      0.03            0.25

编辑：这是我现在使用的代码的开头（输出不正确）：

awk '
NR==FNR{ range[$1,$2,$3,$4,$5]; next }
FNR==1
{
  for(x in range) {
    split(x, check, SUBSEP);
    if($2==check[1] && $3>=check[2] && $3<=check[3]) { print $1"\t"$2"\t"$3"\t"check[4]"\t"check[5]}
  }
}
' File2 File1

但是，当基因不存在时，我没有设法将 0（与“else”）关联起来。我得到错误的行数。你能给我更多的提示吗？非常感谢。

【问题讨论】：

输入文件有多大（MBytes，行数）并且都已按 Chromosome 排序？
您应该以代码的形式展示您所拥有的东西——如果您没有随身携带，请前往您拥有它的地方。
文件有多大？正如@markp-fuso 暗示的那样，为每个文件 1 位置搜索所有文件 2 范围的蛮力方法有点微不足道。如果文件很大，那可能需要很长时间。更复杂的方法会更快，但是您还没有说足够的内容来知道这是否有必要。文件2范围是否重叠的另一件事。您在示例中暗示了这一点，但它总是正确的吗？
有效方法的问题是没有来自File1 的唯一字段组合可以作为File2 中任何字段的键，染色体值不是唯一的，因此不能用作索引（单独或与SUBSEP 和任何其他字段一起使用）。因此，除非您有更多未显示的列，否则您将陷入效率极低的蛮力方法（即使您将File1 读入任何数组并且split() 将每次检查的值与File2 相比较。
文件 1 大约 5MB（~900 行），文件 2 大约 KB（~50 行）。文件 1 大约有 500 列，文件 2 大约有 20 列。文件 1 不完全按染色体排序：我有两个按突变类型的连续系列（1->22，然后又是 1->22）。我没有可用的代码，我会尽快用它编辑我的帖子。文件 2 中的范围不重叠，我的示例暗示多个基因可以具有在同一范围内的位置。确实，由于栏目太多，我只展示了感兴趣的栏目。

标签： bash awk

【解决方案1】：

一个awk-唯一的想法...

注意：请参阅我的其他答案以了解假设/理解以及我的 file1 版本

awk '     # process file2
FNR==NR   { c=$1                            # save chromosome value
            $1=""                           # clear field #1
            file2[c]=$0                     # use chromosome as array index; save line in array
            next
          }

          # process file1
          { start=end=-9999999999           # default values in case
            a1=a2=""                        # no match in file2

            if ($2 in file2) {              # if chromosome also in file2
                split(file2[$2],arr)        # split file2 data into array arr[]

                startpos =arr[1]
                endpos   =arr[2]
                a1       =arr[3]
                a2       =arr[4]
            }

            # if not the header row and file1/position outside of file2/range then set a1=a2=0
            if (FNR>1 && ($3 < startpos || $3 > endpos)) a1=a2=0

            print $0,a1,a2
          }
' file2 file1

这会生成：

Hugo_Symbol Chromosome  Position cell_frac_A1 cell_frac_A2
Gene2       1           222222 0.12 0.01
Gene3       2           333333 0.03 0.25
Gene1       1           111111 0 0
Gene4       2           333337 0.03 0.25
Gene5       4           444567 0 0

将最后一行更改为' file2 file1 | column -t 生成：

Hugo_Symbol  Chromosome  Position  cell_frac_A1  cell_frac_A2
Gene2        1           222222    0.12          0.01
Gene3        2           333333    0.03          0.25
Gene1        1           111111    0             0
Gene4        2           333337    0.03          0.25
Gene5        4           444567    0             0

通过将最后一行更改为' file2 <(head -1 file1; tail -n +2 file1 | sort -k2,3) | column -t 对file1 和Chromosome 和Position 进行预排序，生成：

Hugo_Symbol  Chromosome  Position  cell_frac_A1  cell_frac_A2
Gene1        1           111111    0             0
Gene2        1           222222    0.12          0.01
Gene3        2           333333    0.03          0.25
Gene4        2           333337    0.03          0.25
Gene5        4           444567    0             0

一个大问题（与我的其他答案相同）...在处理总共 519 列时，实际代码可能会变得笨拙，尤其是在需要穿插很多列的情况下；否则 OP 可能能够使用一些 for 循环来更轻松地打印列范围。

【讨论】：

【解决方案2】：

也许是 sql 而非 awk 的工作？

tr -s ' ' '|' <File1 >file1.csv
tr -s ' ' '|' <File2 >file2.csv
(
   echo 'Hugo_Symbol|Chromosome|Position|cell_frac_A1|cell_frac_A2'
   sqlite3 <<'EOD'
.import file1.csv t1
.import file2.csv t2
select distinct
   t1.hugo_symbol,
   t1.chromosome,
   t1.position,
   case
      when t1.position between t2.start_position and t2.end_position
      then t2.cell_frac_a1
      else 0
   end,
   case
      when t1.position between t2.start_position and t2.end_position
      then t2.cell_frac_a2
      else 0
   end
from t1 join t2 on t1.chromosome=t2.chromosome;
EOD
   rm file[12].csv
) | tr '|' '\t'

【讨论】：

只是好奇...对于大量数据...sqlite'smart'足以动态创建（临时）索引（或散列/连接键) 加快加入速度？如果没有，出于连接性能的原因，您可能需要考虑添加一个索引（在内表上）
@markp-fuso 是的，只是想添加一个索引，但解释似乎表明它已经这样做了：SCAN TABLE t2，SEARCH TABLE t1 USING AUTOMATIC COVERING INDEX (Chromosome=?)

【解决方案3】：

基于样本数据和 OP cmets 的假设和/或理解 ...

file1不是按染色体排序
file2是按染色体排序的
两个文件中的共同标头拼写相同（例如，file1:Chromosome 与 file2:Chromosom）
如果file1 中存在染色体，但file2 中不存在不，那么我们保留来自file1 的行，来自file2 的列留空
两个文件都比较小（file1：5MB，900 行；file2：几 KB，50 行）
注意：从繁琐的编码角度来看，列数（file1：500 列；file2：20 列）可能会有问题...稍后会详细介绍...

示例输入：

$ cat file1                           # scrambled chromsome order; added chromosome=4 line
Hugo_Symbol Chromosome  Position
Gene1       1           111111
Gene3       2           333333
Gene2       1           222222
Gene4       2           333337
Gene5       4           444567        # has no match in file2

$ cat file2
Chromosome  Start_Position  End_Position    cell_frac_A1    cell_frac_A2
1           222220          222230          0.12            0.01
2           333330          333340          0.03            0.25
3           444440          444450          0.01            0.01

第一个问题是按Chromosome 和Position 对file1 进行排序，同时保留标题行：

$ (head -1 file1; tail -n +2 file1 | sort -k2,3)

Hugo_Symbol Chromosome  Position
Gene1       1           111111
Gene2       1           222222
Gene3       2           333333
Gene4       2           333337
Gene5       4           444567

我们现在可以join 基于Chromosome 列的2 个文件：

$ join -1 2 -2 1 -a 1 --nocheck-order <(head -1 file1; tail -n +2 file1 | sort -k2,3) file2

Chromosome Hugo_Symbol Position Start_Position End_Position cell_frac_A1 cell_frac_A2
1 Gene1 111111 222220 222230 0.12 0.01
1 Gene2 222222 222220 222230 0.12 0.01
2 Gene3 333333 333330 333340 0.03 0.25
2 Gene4 333337 333330 333340 0.03 0.25
4 Gene5 444567

地点：

-1 2 -2 1 - 加入Chromosome 列：-1 2 == 文件#1 列#2； -2 1 == 文件 #2 列 #1
-a 1 - 保留文件 #1 中的列（已排序 file1）
--nocheck-order - 禁用验证输入按连接列排序；可选的;如果语言环境认为 1 应该在 Chromosome 之前排序，则可能需要
注意： 对于示例输入/输出，我们不需要特殊的输出格式，因此我们可以跳过 -o 选项，但这可能需要根据 OP 对 519 列的输出要求（但它也可能变得笨拙）

从这里 OP 可以使用bash 或awk 进行比较（列#4/#5 之间的列#3）；一个awk想法：

$ join -1 2 -2 1 -a 1 --nocheck-order <(head -1 file1; tail -n +2 file1 | sort -k2,3) file2 | awk 'FNR>1{if ($3<$4 || $3>$5) $6=$7=0} {print $2,$1,$3,$6,$7}'

Hugo_Symbol Chromosome Position cell_frac_A1 cell_frac_A2
Gene1 1 111111 0 0                                         # Position outside of range
Gene2 1 222222 0.12 0.01
Gene3 2 333333 0.03 0.25
Gene4 2 333337 0.03 0.25
Gene5 4 444567 0 0                                         # no match in file2; if there were other columns from file2 they would be empty

为了匹配 OP 的示例输出（似乎是固定宽度要求），我们可以将其传递给 column：

$ join -1 2 -2 1 -a 1 --nocheck-order <(head -1 file1; tail -n +2 file1 | sort -k2,3) file2 | awk 'FNR>1{if ($3<$4 || $3>$5) $6=$7=0} {print $2,$1,$3,$6,$7}' | column -t

Hugo_Symbol  Chromosome  Position  cell_frac_A1  cell_frac_A2
Gene1        1           111111    0             0
Gene2        1           222222    0.12          0.01
Gene3        2           333333    0.03          0.25
Gene4        2           333337    0.03          0.25
Gene5        4           444567    0             0

注意：请记住，对于 OP 的 519 列，这可能是站不住脚的，特别是如果散布的列包含空格/空格（即，column -t 可能无法正确解析输入）

问题（除了不正确的假设和之前的注意事项）：

对于相对较小的文件，join | awk | column 的性能应该足够了
对于较大的文件，所有这些代码都可以合并到单个 awk 解决方案中，尽管内存使用可能是小型机器上的问题（例如，awk 的一个想法是通过数组将 file2 加载到内存中所以内存需要足够大以容纳所有 file2 ... 可能不是问题，除非 file2 的大小达到 100/1000 MB）
对于总共 519 列，awk/print 将变得笨拙，尤其是在需要移动/散布大量列时

【讨论】：