ヒトの染色体ごとのファイルを処理する？
ヒトのゲノム情報等のファイルは一般的に非常に巨大になるので、染色体ごとに処理したくなるのが人情ですよね。
Linuxをほとんど触ったことが無かった時に、染色体ごとにファイルを処理することは様々なLinuxのブレース展開を覚える良い勉強になったので、まとめてみました。
以下は、bashが使用できる環境（Linux, MacまたはWindowsのWindows Subsystem for Linux 2）を想定しています。
Windows Subsystem for Linux 2の導入については、こちらを参照しました。感謝感謝。
ファイルの繰り返し処理
ヒトのゲノムに存在するバリアントは、全ゲノムシークエンスデータからだと１個人あたり数百万個も検知されます。
7,609人の全ゲノムシークエンスデータからだと、約8000万弱のバリアントが検知されています。ﾋﾄﾊﾐﾅﾁｶﾞｳ。
上述したバリアントデータには頻度情報も含まれますが、データの大きさは約1.8GBにも上ります。
ちなみに、データはこちらからダウンロード可能。
これらのように染色体ごとに分割されたファイルに対して同じような処理を適用することを考えます。
# ヒトの染色体の番号を表示する
ヒトのゲノムの染色体には、常染色体は1番から22番、性染色体はX、Y、ミトコンドリアが存在しています。
一般的には、染色体番号は1~22、X、Y、MTと表記されますが、SNPの解析で良く用いられるplinkでは X⇒23、Y⇒24、XのPseudo-autosomal region⇒25、MT⇒26 なんてこともありますし、chr1~chr22, chrX, chrY, chrMTと表示されることもあります。。
今、染色体の番号が1~26までのファイルがあり、これに対してブレース展開というものを使うと、とても便利です。
$ echo {1..26}
と打つと
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
となる。
便利！
これを使って、for文で以下のように繰り返し処理すれば、
for chr in {1..26}
do
cat なんちゃらファイル.chr${chr}.vcf
done
これで、なんちゃらファイル.chr1.vcf から なんちゃらファイル.chr26.vcf までのファイルを対象に、同じコマンド（ここでは、catコマンド）を実行できます。
# ヒトの染色体は 1~22、X、Y、MT じゃろがい
そうですよね。染色体の番号が1~26までなんて都合が良いことは無いですよね。
ということで、何とか、1~22、X、Y、MTを出力するように頑張りましょう。
1~22までは
$ echo {1..22}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
で良いとして、これの後に X、Y、MTを出力することを考えます。
ブレース展開で、{X,Y,MT} と書くと
$ echo {X,Y,MT}
X Y MT
となります。ここまでくれば出来たも同然。二つのブレース展開をつなぐと、
$ echo {1..22} {X,Y,MT}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X Y MT
パチパチパチ。
よって、
for chr in `echo {1..22} {X,Y,MT}`
do
cat なんちゃらファイル.chr${chr}.vcf
done
と、書いてあげると 1番染色体から22番染色体、X、Y、MTのファイルにコマンドを実行（ここでは cat）出来ることになります。
↓は単純な応用例で、上記の出力の先頭にchrが付くはずです。
for chr in `echo chr{1..22} chr{X,Y,MT}`
do
cat なんちゃらファイル.chr${chr}.vcf
done
このようにブレース展開を利用すると、解析のプログラムもすっきり書けるようになると思います。
それでは良いゲノム解析ライフを！