こんにちはー hr-kです。
今回はみなさんご存知のsnpEffの機能について書いていきますー
みなさんよく使うアノテーション機能ですが、よく知っているvcfへのアノテーション
$ java -jar snpEff.jar hg19 test.vcf > test_ann.vcf
やbedファイルへのアノテーション
$ java -jar snpEff.jar -i bed -o bed hg19 test.bed > test_ann.bed
以外にも、ちょっと面白い使い方ができます。
遺伝子名からのbedファイル生成
ゲノム解析をしているみなさんは「あれ？この遺伝子どこにあったっけ？」とか「解析したい遺伝子は決まってるけど領域に直すのか。。。辛いな。。。」と思うことが度々起きると思います（僕は半年に一回くらい定期的にありますね）。
そんなときに便利なのがsnpEffのgenes2bedを用いた「遺伝子シンボル ⇨ BED」の変換です！
例えばhg19のゲノムでNOTCH2とNOTCH3の位置を知りたい時は次のコマンドを実行します。
$ java -jar snpEff.jar genes2bed hg19 NOTCH2 NOTCH3 > NOTCH2-3.bed
するとーNOTCH2-3.bedに遺伝子に対するbedが出来上がるのです！
#chr start end geneName;geneId;strand 1 120454175 120612317 NOTCH2;NOTCH2;- 19 15270443 15311792 NOTCH3;NOTCH3;-
これで特定のゲノム配列について遺伝子がある領域を調べることができますね。
文法は
$ java -jar snpEff.jar genes2bed [genome] [領域を見たい遺伝子シンボルを半角スペース区切りで] > [出力するBEDファイル]
となっています。お試しあれー
ちょっと応用
たくさんの遺伝子についてのBEDを作りたい時は、見たい遺伝子シンボルに対してこんなファイル（genes.txt）を用意してー
NOTCH1 NOTCH2 NOTCH3 BRCA1 BRCA2 AKT
こんな感じで（※）渡しちゃいます！
$ java -jar snpEff.jar genes2bed hg19 $(cat genes.txt) > genes.bed
できたgenes.bedをのぞいてみましょー。
#chr start end geneName;geneId;strand 1 120454175 120612317 NOTCH2;NOTCH2;- 9 139388884 139440238 NOTCH1;NOTCH1;- 13 32889616 32973809 BRCA2;BRCA2;+ 17 41196311 41277500 BRCA1;BRCA1;- 19 15270443 15311792 NOTCH3;NOTCH3;-
はい、素晴らしいですねー。
何が素晴らしいって、入れた遺伝子シンボルの順じゃなくてソートして染色体順=>ポジション順に自動で並べ替えられて出てくるのです！
これで遺伝子シンボルから簡単にBEDファイルを作れそうですね〜（ただしsnpEffでdatabaseがあるゲノム、もしくは作れるゲノムに限る）
実は「遺伝子シンボル => bed」の変換は、RでもbiomaRt使ってできるんですが、snpEffの方が大分簡単にできます。
興味ある人はRでもやってみてね！（面倒だけどね！）
※： cat genes.txt の結果を変数として使用するが、変数にする際に改行がなくなる。こんな感じ。
$ echo $(cat genes.txt) NOTCH1 NOTCH2 NOTCH3 BRCA1 BRCA2 AKT
もちろん、初めから
NOTCH1 NOTCH2 NOTCH3 BRCA1 BRCA2 AKT
を用意するのもあり。
ではあるが、人が見るものである以上、リストとして用意するのは見やすいように前者の「一行一遺伝子のファイル」であることが多い気がする。
