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が変換されますが、他のフィールドは変更されずに保持されます。

awkapply を使用してgsub()入力行全体を変更する方法、または n番目のフィールドを具体的に変更する方法は知っていますが、変更する必要があるのは、データ ファイルの最初の行の列ラベルで識別されるように、位置ファイルにリストされているフィールドのみです。どうすればそれを実装できawkますか?

4

2 に答える 2

1

正確なフィールド (列) 区切り記号を保持することが必須ではない場合、つまり、すべての列区切り記号を単一のスペースなどの固定文字列に自由に変更できる場合gsub()、フィールドごとに使用できます。その後、レコードを再構築します。これにより、特定のフィールドへの変更を制限するという問題が解決されます。

もう 1 つの問題は、位置ファイルのデータと列見出しに基づいて、変更するフィールドを特定することです。これを行う1つの方法は次のとおりです。

  • ブロックを使用してBEGIN、位置ファイルから各行を読み取り、その内容を配列インデックスとして記録します。これは、各行の内容をハッシュ テーブルに記録していると考えることができます。

  • フィールドをループし、ラベルの配列に存在するかどうかを確認することにより、事前に読み取られた列ラベルをメイン入力の最初の行から読み取られた列ラベルと一致させます。存在するものについては、フィールド番号をインデックスとして 2 番目の配列に記録します。

  • 後続の行ごとに、その構成フィールドからレコードを再構築し、元の値と、フィールド番号が変換が必要なフィールドの 1 つとして記録されているかどうかに基づいて変更された値のいずれかを選択します。

  • 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 に答える