使用第二个文件评估 awk 中 if 语句的条件

Question

我正在尝试编写一个操作两个文件的 bash 脚本：

文件1

Region  Coords. RsId    Position    Alleles Disease PValue  OddsRatio   RegionID
1p13.2  1:113839149..114551845  rs2476601   114377568   G>A Alopecia Areata 8.90E-08    1.34    869
2q13    2:111444884..111809030  rs3789129   111698040   A>C Alopecia Areata 1.50E-08    0.76    871
2q33.2  2:204611195..204817281  rs3096851   204763882   A>C Alopecia Areata 3.58E-08    1.32    802
2q33.2  2:204611195..204817281  rs1024161   204721752   G>A Alopecia Areata 3.55E-13    1.44    802
2q33.2  2:204611195..204817281  rs231775    204732714   A>G Alopecia Areata 2.20E-20    1.39    802
4q27    4:122982314..123605528  rs7682241   123523875   C>A Alopecia Areata 4.27E-08    1.34    803
4q27    4:122982314..123605528  rs7682481   123524026   G>C Alopecia Areata 4.80E-09    1.23    803
5q31.1  5:131783213..132135372  rs848   131996500   C>A Alopecia Areata 4.80E-09    1.27    872

文件2

CHR_A         BP_A       SNP_A  CHR_B         BP_B       SNP_B           R2 
 2    204721752   rs1024161      2    204732714    rs231775     0.849535 
 2    204721752   rs1024161      2    204763882   rs3096851      0.68029 
 2    204732714    rs231775      2    204763882   rs3096851     0.739633 
 4    123523875   rs7682241      4    123524026   rs7682481            1

我想读取 file1，如果 file2 的第 3 列 (SNP_A) 或第 4 列 (SNP_B) 中不存在第 3 列 (RsId) 值，则将该行写入输出。我尝试了以下方法：

#this executable file is called filter_file.sh

#!/bin/bash
file1=$1
file2=$2
outfile=$3

while read CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2; do 
cat $file2 | awk "(\$3!~/$SNP_A|$SNP_B/) {print}"
done < $file1 > $outfile

./filter_file.sh file1 file2 out

当我单独测试 awk 语句时它有效，但是当我将它添加到 bash while 循环时，它打印了所有 file1，包括标题，四次。 这一步的代码有什么问题？

一旦这工作正常，如果 file1 的第 3 列 (RsId) 值存在于 file2 的第 3 列 (SNP_A) 或第 4 列 (SNP_B) 中，我想将具有最低值的行写入 file1 第 7 列的输出（P值）。

我不确定如何开始任务的第二部分。通过阅读其他 awk 问题，我想我可以尝试使用 if 语句设置如下内容：

#!/bin/bash
file=$1
file2=$2
outfile=$3

while read CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2; do 
cat $file2 | awk "{
if ((\$3!~/$SNP_A|$SNP_B/))
    print $0;
else
    #Statement that prints only the row with the lowest value for column 7 
}"
done < $file1 > $outfile

我可以使用哪些方法来执行此步骤？

如果人们可以指出一些可能对解决这类问题有帮助的教程，我们将不胜感激。

所需的输出文件将如下所示（顺序无关紧要）：

Region  Coords. RsId    Position    Alleles Disease PValue  OddsRatio   RegionID
1p13.2  1:113839149..114551845  rs2476601   114377568   G>A Alopecia Areata 8.90E-08    1.34    869
2q13    2:111444884..111809030  rs3789129   111698040   A>C Alopecia Areata 1.50E-08    0.76    871
2q33.2  2:204611195..204817281  rs231775    204732714   A>G Alopecia Areata 2.20E-20    1.39    802
4q27    4:122982314..123605528  rs7682481   123524026   G>C Alopecia Areata 4.80E-09    1.23    803
5q31.1  5:131783213..132135372  rs848   131996500   C>A Alopecia Areata 4.80E-09    1.27    872

Answer 1

无需使用 bash 循环文件，您可以在 awk 中完成所有操作：

$ awk 'NR == FNR && NR > 1 {snp_a[$3]; snp_b[$6]; next}
       FNR > 1 && !($3 in snp_a || $3 in snp_b)' file2 file1
1p13.2  1:113839149..114551845  rs2476601   114377568   G>A Alopecia Areata 8.90E-08    1.34    869
2q13    2:111444884..111809030  rs3789129   111698040   A>C Alopecia Areata 1.50E-08    0.76    871
5q31.1  5:131783213..132135372  rs848   131996500   C>A Alopecia Areata 4.80E-09    1.27    872

第一个块应用于第一个文件 (file2) 并设置与感兴趣的两列对应的数组中的键。 next 跳过任何其他命令，因此脚本的其余部分仅适用于第二个文件 (file1)。当条件为真时打印行，即文件中的行号 (FNR) 大于 1，并且在两个数组中的任何一个中都找不到键。

---

对于问题的第二部分，事情变得有点复杂......希望 cmets 解释一下：

$ cat script.awk
# first file, save individual columns and pairs
NR == FNR && NR > 1 {snp_a[$3]; snp_b[$6]; pair[$3,$6]; next}

# second file, print first line
NR != FNR && FNR == 1

# second file, rest of lines
FNR > 1 {
    # print lines which aren't in either array
    if (!($3 in snp_a || $3 in snp_b)) {print}
    # save other lines and corresponding p-value
    else {s[$3] = $0; p[$3] = $8}
}
END {
    # loop through all lines
    for (i in s) {
        # empty the array f
        split("", f)
        # set initial line and min
        line = s[i]
        min = p[i]

        # locate associated lines
        if (i in snp_a) {
            for (j in snp_b) {
                # SUBSEP is a special variable used when combining keys
                # as in first block pair[$3,$6]
                if (i SUBSEP j in pair && j in s && p[j] < min) {
                    min = p[j]
                    line = s[j]
                }
            }
        }
        else if (i in snp_b) {
            for (j in snp_a) {
                if (j SUBSEP i in pair && j in s && p[j] < min) {
                    min = p[j]
                    line = s[j]
                }
            }
        }
        out[line]
    }
    for (i in out) print i
}
$ awk -f script.awk file2 file1
Region  Coords. RsId    Position    Alleles Disease PValue  OddsRatio   RegionID
1p13.2  1:113839149..114551845  rs2476601   114377568   G>A Alopecia Areata 8.90E-08    1.34    869
2q13    2:111444884..111809030  rs3789129   111698040   A>C Alopecia Areata 1.50E-08    0.76    871
5q31.1  5:131783213..132135372  rs848   131996500   C>A Alopecia Areata 4.80E-09    1.27    872
2q33.2  2:204611195..204817281  rs231775    204732714   A>G Alopecia Areata 2.20E-20    1.39    802
4q27    4:122982314..123605528  rs7682481   123524026   G>C Alopecia Areata 4.80E-09    1.23    803

我很确定END 块中的逻辑可以简化，但我想不出更好的方法来确定“对”。无论哪种方式，它都能达到您想要的输出。

Answer 2

您可以为此使用awk：

awk 'NR==FNR&&NR>1{r[NR]=$3$6;next}{p=1;for(i in r){if(r[i]~$3){p=0;break}}}p' \
  file2 file1

在包含 cmets 的多行版本中更容易理解：

# NR (row number) and FNR (current input file's row number) are equal
# only as long as we are processing the first input file.
# Means the following block only runs on file2 which is getting passed first.
# NR > 1 makes sure that we skip the header line.
NR==FNR && NR>1 {
    # Push $3 concatenated with $6 to an array call r (rows)
    r[NR]=$3$6
    # Stop processing this line and use the next. This makes sure
    # that the following block will only operate on file1
    next
}

# The following block runs on every line of file1
{
    # Init p flag
    p=1 
    # Iterate to the previously stored values
    for(i in r) {
        # If the entry of r is not matching field 3
        # Print the current line.
        if(r[i]~$3){
            # If a line matches we reset the p flag
            # and break the loop
            p=0
            break
        }
    }   
}

# Print if p flag isset (print is the default action in awk)
p

输出：

CHR_A         BP_A       SNP_A  CHR_B         BP_B       SNP_B           R2 
Region  Coords. RsId    Position    Alleles Disease PValue  OddsRatio   RegionID
1p13.2  1:113839149..114551845  rs2476601   114377568   G>A Alopecia Areata 8.90E-08    1.34    869
2q13    2:111444884..111809030  rs3789129   111698040   A>C Alopecia Areata 1.50E-08    0.76    871
5q31.1  5:131783213..132135372  rs848   131996500   C>A Alopecia Areata 4.80E-09    1.27    872