0

我有一个包含 DNA 碱基序列表的原始文件，带有行和列标签，以及一个单独的“位置”文件，其中列出了列标签的子集。我需要处理原始文件，对位置文件标识的列中的值执行转换。

示例原始文件：

name pos1 pos2 pos3 pos4 pos5 pos6 pos7
name1 AT TA CT GT CC TC TT
name2 AA TA TT GT TC TC TT
name3 AT TT CG AT CT TC TT
name4 GT TA CT TT CC TC TT

位置文件示例：

pos1
pos3
pos6
pos7

在每个选定的字段上，我需要执行这些翻译：

A to T
C to G
G to C
T to A

因此，通过基于提供的位置文件处理示例原始文件获得的输出将是：

name pos1 pos2 pos3 pos4 pos5 pos6 pos7
name1 TA TA GA GT CC AG AA
name2 TT TA AA GT TC AG AA
name3 TA TT GC AT CT AG AA
name4 CA TA GA TT CC AG AA

因此第一行未修改，在随后的每一行中，对应于列标签pos1、pos3、pos6和pos7的字段都进行了转换，而其他字段保持不变。

我知道如何使用awkapplygsub()来修改整个输入行或专门修改第 n个字段，但我只需要修改位置文件中列出的那些字段，由数据文件第一行上的列标签标识。我怎样才能实现awk呢？

4

2 回答 2

1

假设不必保留确切的字段（列）分隔符——也就是说，您可以自由地将每个列分隔符更改为固定字符串，例如单个空格——您可以gsub()在每个字段的基础上使用然后重建记录。这解决了将更改限制为特定字段的问题。

另一个问题是根据位置文件中的数据和列标题确定要修改的字段。这是一种方法：

  • 使用BEGIN块，从位置文件中读取每一行并将其内容记录为数组索引。您可以将其视为在哈希表中记录每一行的内容。

  • 通过遍历字段并检查它们在标签数组中的存在，将预读取的列标签与从主输入的第一行读取的列标签匹配。对于那些存在的，将字段编号记录为第二个数组中的索引

  • 对于每个后续行，从其组成字段重建记录，根据字段编号是否记录为需要转换的字段编号，在原始值和修改后的值之间进行选择。

  • 特别注意awk可以通过存储在变量中的字段编号来引用字段。因此，myfield = 2; print $myfield产生与print $2

执行所有操作的awk程序可能如下所示：

#!/usr/bin/awk

function baseswap(seq) {
  gsub(/A/, "X", seq)
  gsub(/T/, "A", seq)
  gsub(/X/, "T", seq)
  gsub(/C/, "X", seq)
  gsub(/G/, "C", seq)
  gsub(/X/, "G", seq)
  return seq
}

BEGIN  {
         while ((getline < "positions") == 1) {
           labels[$1] = 1
         }
       }

FNR==1 {
          for (i = 2; i <= NF; i++) {
            if (labels[$i]) {
              fields[i] = 1
            }
          }
          print
          next
       }

       {
         record = $1
         for (i = 2; i <= NF; i++) {
           record = record " " (fields[i] == 1 ? baseswap($i) : $i)
         }
         print record
       }
于 2016-06-10T16:53:36.820 回答
1 
$ cat tst.awk
BEGIN {
    split("A T C G G C T A",t)
    for (i=1;i in t;i+=2) {
        map[t[i]] = t[i+1]
    }
}
NR==FNR {
    fldNames[$1]
    next
}
FNR==1 {
    for (i=1;i<=NF;i++) {
        if ($i in fldNames) {
            targets[i]
        }
    }
}
FNR>1 {
    $0 = tolower($0)
    for (fldNr in targets) {
        for (old in map) {
            gsub(tolower(old),map[old],$fldNr)
        }
    }
    $0 = toupper($0)
}
{ print }

$ awk -f tst.awk positions original
name pos1 pos2 pos3 pos4 pos5 pos6 pos7
NAME1 TA TA GA GT CC AG AA
NAME2 TT TA AA GT TC AG AA
NAME3 TA TT GC AT CT AG AA
NAME4 CA TA GA TT CC AG AA
于 2016-06-10T18:14:10.530 回答