UNIX for Bioinformatics

5. awk

awk は、CSV や TSV ファイルなどのようなデータを処理するときに使うスクリプト言語の一つである。awk はコマンドや bash と同時に使うことができ、非常に手軽である。データの解析に先立って、平均や分散を見積もったりするときに利用すると便利である。ただし、本格的な解析の場合は、awk を使わずに、Python や R などの強力なプログラミング言語を利用するとよい。このページでは awk の基本的な使い方を紹介する。

awk 基本文型

awk の基本文型はパターンとアクションを 1 組として並べた文型をとる。例えば、あるファイルの各行の内容に対して、その行が指定パターン（pattern）とマッチするならば、処理（action）を行う、といった処理は次のように書く。

pattern { action }

パターンと処理が複数あるとき、この基本文型を続けて書けばよい。例えば、 行の内容が pattern1 にマッチするならば処理 1（action1）を行い、pattern2 にマッチするならば処理 2（action2）を行う、といった処理を行う場合は、次のように書く。

xxxxxxxxxx
pattern1 { action1 }
pattern2 { action2 }

複数のパターンとアクションを使う場面として、例えば、"apple" を含むならば処理 1、"orange" を含むならば処理 2 のように、パターン分けしたいときに使う。

パターン

パターンは正規表現で指定する。例えば、ある行の内容全体が "apple" を含むならば処理 1、"orange" を含むならば処理 2 を行うためには、パータンを次のように指定する。パターンマッチの対象を / で囲む。

xxxxxxxxxx
$0 ~ /apple/ { action1 }
$0 ~ /orange/ { action2 }

なお、$0 は組み込み変数とよばれる特殊な変数である。$0 に 1 行の内容が保存されている。$0 ~ /apple/ は、その行のどこかに "apple" という単語が含むかどうかを確認している。そして、その行の中でどこかに "apple" という単語が見つかれば、action1 が実行される。なお、行全体（$0）に対するパターンマッチの場合、$0 ~ を省略することができる。例えば、以下のように書くこともできる。

xxxxxxxxxx
/apple/ { action1 }
/orange/ { action2 }

CSV あるいは TSV データには列（フィールド）の概念がある。このとき、例えば 2 列目のデータが "apple" を含むならば処理 1、2 列目のデータが "orange" を含むならば処理 2 を行うには、次のように 2 列目のみに対してパターンマッチを行う。$2 は組み込み変数と呼ばれる特殊な変数であり、2 列目を意味する。同様に、$3、$4 などはそれぞれ 3 列目や 4 列目を意味する。

xxxxxxxxxx
$2 ~ /apple/ { action1 }
$2 ~ /orange/ { action2 }

正規表現を使ったパターンマッチ以外に、パターンに行数指定ができる。例えば、11 行目以降のデータに対して処理を行う場合は、次のように書く。なお、NR は組み込み変数の一つで、現在ファイルの何行目を読んでいるのかの情報が保存されている。

xxxxxxxxxx
NR > 10 { action }

条件が複数ある場合、AND 演算子 && あるいは OR 演算子 || で条件をつなげることができる。例えば、3 列目に "apple" を含む行でかつ 10 行目以降の行に対して処理を行う場合、次のように書く。

xxxxxxxxxx
$3 ~ /apple/ && NR > 10 { action }

アクション

アクションには、加減乗除の計算処理や出力処理を行うことができる。また、必要に応じて、if、while、for などのフロー制御もできる。ただ、フロー制御まで使わなければならない処理については、awk に固執するよりも、Perl、Python や R で書いた方が断然効率がいい。そのため、ここでは、awk のフロー制御について詳細に紹介しない。

アクションの書き方を単独に紹介するよりも、実例を使った方がわかりやすい。以下で、実例を使いながら、アクションの使い方を紹介する。

使用例

列の抽出

awk を使用して、ある条件を満たした列のデータだけを抽出する例を説明する。この操作は、grep および cut コマンドでも実現できる。ここで、具体的に、data/field_data ディレクトリの中にある iris.txt ファイルに対して、"setosa" を含む行のうち、3 列目のデータを出力してみる。

​x
cd
cd unix4bi/data/field_data
grep setosa iris.txt | cut -f 3
## 3.5
## 3
## 3.2
## ...
## ...
​
awk '/setosa/ { print $3 } ' iris.txt
## 3.5
## 3
## 3.2
## ...
## ...

awk では、1 列目、2 列目、3 列目、・・・のデータはそれぞれ組み込み変数の $1、$2、$3、・・・の中に保存される。3 列目のデータを出力させたいので、アクションを print $3 のように書く。

次に、同じく setosa を含む列のうち、2 列目から 5 列目までのデータを出力する例を示す。この際、列と列の間をカンマで区切る必要がある。カンマで区切ることで、実際に出力される文字列は空白 1 つ分で区切られる。

xxxxxxxxxx
awk '/setosa/ { print $2, $3, $4, $5 } ' iris.txt
## 5.1 3.5 1.4 0.2
## 4.9 3 1.4 0.2
## 4.7 3.2 1.3 0.2
## ...
## ...

一方で、カンマで区切らない場合は、次のように、列と列の間は空白がなく、連結した状態で出力される。また、各列の情報の他に、付け足したい文字列があれば、">" のようにダブルクォートで囲んで付ける。

xxxxxxxxxx
awk '/setosa/ { print ">" $2 $3 $4 $5 } ' iris.txt
## >5.13.51.40.2
## >4.931.40.2
## >4.73.21.30.2
## ...
## ...

平均の計算

次に、awk を使って平均の計算に挑戦してみる。data/field_data ディレクトリの中にある iris.txt ファイルに対して、"setosa" を含む行のうち、3 列目のデータのデータに対して、平均を計算してみる。平均を計算するため、setosa の各行のデータの和と行数を保持しておく必要があるので、ここではそれぞれ変数 sum と n に保存する。

xxxxxxxxxx
awk '/setosa/ { sum = sum + $3; n = n + 1 } END { print sum/n }' iris.txt
## 3.428
​
awk '/setosa/ { sum = sum + $3; n = n + 1 } END { print "sum=" sum/n }' iris.txt
## sum=3.428

/setosa/ の付属処理では sum = sum + $3 および n = n + 1 という二つの処理を行っている。これらは、setosa を含む行のデータに対して、前者は 3 列目のデータの合計を計算し、後者は行数を計上している。

また、ここでは END というパターンが新たに使われている。この END は特別な意味を持っている。それは、ファイルに対する処理がすべて完了してから、END の付属アクションを実行せよ、という意味である。この場合、iris.txt ファイルに対する処理を終わった後に、これまでに得られた変数 sum および n を使って、最後に sum/n を計算して、出力してることになる。

同様に、4 行目、5 行目についても計算を行ってみる。

xxxxxxxxxx
awk '/setosa/ { sum = sum + $4; n = n + 1 } END { print sum/n }' iris.txt
## 1.462
​
awk '/setosa/ { sum = sum + $5; n = n + 1 } END { print sum/n }' iris.txt
## 0.246

上の例では、各列のデータをそれぞれ独立に平均を求めてきた。すべての列を同時に扱いたい場合は、sum2、sum3、sum4 などのようにそれぞれの列の合計を保存するための変数を用意する必要がある。以下では、2 列目と 3 列目の平均を同時に求める例を示す。

xxxxxxxxxx
awk '/setosa/ { sum2 = sum2 + $2; sum3 = sum3 + $3; n = n + 1 } END { print sum2/n, sum3/n }' iris.txt
## 5.006 3.428

列が多くなると、この書き方では大変になる。そこで、for 構文を使って、1 行を処理すると後に 2 列目から 5 列目のデータを繰り返し変数に足す処理を行う。

xxxxxxxxxx
awk '/setosa/ { for(i=2; i<6; i++){ sum[i] = sum[i] + $i }; n = n + 1 } END { for(i=2; i<6; i++){ print sum[i] / n } }' iris.txt 
## 5.006
## 3.428
## 1.462
## 0.246

各列に対して平均値を求めるという単純な計算にもかかわらず、awk のコードがこのように複雑になってきた。このように if や for などのフロー制御構文を使うような場合は、Perl、Python や R などを使うと簡単である。

行番号で行を抽出

次に行数を指定して、該当する行を抽出する例を示す。例えば、iris.txt ファイルの 20 行目から 30 行目までのデータを抽出したいときは、次のようにする。この処理は head および tail コマンドでも実現できる。

xxxxxxxxxx
awk '20 <= NR && NR <= 30 { print $0 }' iris.txt
## 20 5.1 3.8 1.5 0.3 setosa
## 21 5.4 3.4 1.7 0.2 setosa
## 22 5.1 3.7 1.5 0.4 setosa
## 23 4.6 3.6 1 0.2 setosa
## 24 5.1 3.3 1.7 0.5 setosa
## 25 4.8 3.4 1.9 0.2 setosa
## 26 5 3 1.6 0.2 setosa
## 27 5 3.4 1.6 0.4 setosa
## 28 5.2 3.5 1.5 0.2 setosa
## 29 5.2 3.4 1.4 0.2 setosa
## 30 4.7 3.2 1.6 0.2 setosa

ここで使っているパターンは 20 <= NR && NR <= 30 であるので、行数が 20 以上 30 以下ということになる。これらの行に対して、アクション print $0 を行っている。この $0 は、上でも述べたように、組み込み変数の一つで、行全体の内容が保存されている。ここまでに出てきた例を確認すると、awk では、行全体のデータを $0 に保存し、続けて 1 列目、2 列目、3 列目、・・・のデータをそれぞれ $1、 $2、 $3、 ・・・の中に保存していることになる。

偶数行を抽出

awk を実行しているとき、処理中の行の番号が組み込み変数 NR に保存されている。この NR が奇数か、偶数かを判断することで、奇数行のみあるいは偶数行のみを抽出できるようになる。偶数行のみを抽出したければ、NR を 2 で割った余りがゼロかどうかを判定すればよい。

xxxxxxxxxx
awk 'NR % 2 == 0 { print $0 }' iris.txt
## 2  4.9 3 1.4 0.2 setosa
## 4  4.6 3.1 1.5 0.2 setosa
## 6  5.4 3.9 1.7 0.4 setosa
## ...
## ...

空行を削除

空行とは、行の先頭と末尾の間に文字がない行のことである。このような行を表すパターンとして、「行は先頭と末尾が存在するが、それ以外の文字がない」というパターンを作ればよい。正規表現では、行の先頭は ^ で表され、末尾が $ で表される。したがって、行の先頭から末尾に文字が存在しないという正規表現は ^$ で記述される。このパターンを使って空行を削除してみる。

行内容の全体があるパターンにマッチするときは、次のように $0 ~ /pattern/ あるいは $0 を省略して /pattern/ のように記述できる。

xxxxxxxxxx
$0 ~ /pattern/ {action }
/pattern/ { action }

これに対して、あるパターンにマッチしないときは ~ の代わりに、!~ を使用する。例えば、空行パターン ^$ にマッチしない行に対して、出力を行う場合は次のようにする。

xxxxxxxxxx
awk '$0 !~ /^$/ { print $0 }' iris.txt

練習問題 A1

問題 A1-1

unix4bi/data/pdb のディレクトリに PDB からダウンロードしたテキストファイル 1wkp.cif がある。このファイルにはタンパク質の立体構造の座標データが含まれている。立体構図の x, y, z 座標データは、cif ファイルの中で、先頭の 4 文字が ATOM の行に含まれている。次に示している情報が、その一部である。

xxxxxxxxxx
ATOM   1    N N   . ARG A 1 9   ? 54.205 11.388  56.124 1.00 27.06 ? 6    ARG A N   1
ATOM   2    C CA  . ARG A 1 9   ? 54.608 10.170  55.335 1.00 27.60 ? 6    ARG A CA  1
ATOM   3    C C   . ARG A 1 9   ? 55.808 10.418  54.415 1.00 27.35 ? 6    ARG A C   1
ATOM   4    O O   . ARG A 1 9   ? 55.683 10.337  53.204 1.00 27.82 ? 6    ARG A O   1

ATOM の行の情報のうち、4 列目が原子の種類を示している。C は炭素、N は窒素である。また、CA はアルファ炭素のことを表している。続けて、11 列目、12 列目、13 列目はタンパク質を構成する各原子の立体座標 x、y、z のデータとなっている。そして、19 列目はチェーン名（A、B、・・・）になっている。 

チェーン A のアルファ炭素の立体座標データは次のように取得できる。ここで、awk の部分を修正して、チェーン A のアルファ炭素の重心を求めよ。

xxxxxxxxxx
awk '$1 ~ /ATOM/ && $4 ~ /CA/ && $19 ~ /A/ { print $11,$12,$13 }' 1wkp.cif
## 54.608 10.170 55.335
## 58.196 10.963 54.262
## 61.511 8.962 54.256
## ...
## ...

問題 A1-2

data/field_data/yield_data ディレクトリに複数の収穫量データが含まれている。これらのデータ中で、akitakomachi の収量の平均値（全県前年度分のデータの和を計算し、それを件数で割った値）を求めよ。

なお、すべてのファイルに対して一括に awk で処理したい場合は、次のように cat コマンドですべてのファイルの内容を出力し、それらの出力内容をパイプで awk コマンドに渡すことで、処理できる。

xxxxxxxxxx
cat *.txt | awk '/akitakomachi/ { print $0 }'
## Akita  akitakomachi  599
## Iwate  akitakomachi  552
## Ibaraki  akitakomachi  546
## ...
## ...

問題 A1-3

data/fastq ディレクトリにある leaf_1.fastq （あるいは leaf_1.fq）は FASTQ フォーマットである。awk を使って、この FASTQ ファイルを FASTA ファイルに変換せよ。なお、FASTQ ファイルの配列 ID の先頭に付く "@" を ">" に変換しなくても良い。

xxxxxxxxxx
​
​
​
​

問題 A1-4

projects/mouse ディレクトリにある mortazavi_count_table.txt ファイルはタブ区切りのファイルとなっている。このファイルの 2 列目の数値の合計を計算せよ。ただし、このファイルの 1 行目はヘッダー行となっていることに注意せよ。

xxxxxxxxxx
​
​
​
​