アメリエフのブログ

Amelieff Staff Blog

遺伝子シンボルから遺伝子領域ファイル(BED)生成

こんにちはー 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

を用意するのもあり。
ではあるが、人が見るものである以上、リストとして用意するのは見やすいように前者の「一行一遺伝子のファイル」であることが多い気がする。