参照ゲノムからコード領域の座標を抽出する
- specificな目的ですが、非モデル生物のエクソーム解析をスピーディに進めるうえで便利なので、Tipsとして紹介します。
- データベースが完備されているゲノムから、タンパクコード領域（CDS）の領域の座標を一括してリスト化する方法です。
SnpEff
- SnpEffに登録されている参照ゲノムを使います。2021年12月時点で、v5.0eが利用可能です。
- 例としてヒトに最も近縁な生物であるチンパンジー（学名Pan troglodytes）のCDSをリスト化してみます。
- まず、チンパンジーのデータベース名を確認します。
java -jar snpEff.jar databases > snpEff_database.list # 全データベースをリスト化
grep 'Pan_troglodytes' snpEff_database.list # チンパンジーを含む行を検索
$ CHIMP2.1.4.75 Pan_troglodytes
https://snpeff.blob.core.windows.net/databases/v5_0/snpEff_v5_0_CHIMP2.1.4.75.zip
$ Pan_tro_3.0.99 Pan_troglodytes
https://snpeff.blob.core.windows.net/databases/v5_0/snpEff_v5_0_Pan_tro_3.0.99.zip
- チンパンジーゲノムのデータベースは
CHIMP2.1.4.75とPan_tro_3.0.99の2種類があって、後者の方が最新なので、こちらをダウンロードします。
java -jar snpEff.jar download Pan_tro_3.0.99
-
genes2bed -cdsで、CDSを取り出してbed形式で保存できます。
java -Xmx8g -jar snpEff.jar genes2bed -cds Pan_tro_3.0.99 > Pan_tro_3.0.99_cds.bed
- 中身を見てみると、並びが必ずしもソートされておらず、領域の重複もあります。
head -6 Pan_tro_3.0.99_cds.bed
$ #chr start end geneName;geneId;transcriptId;exonRank;strand
$ 1 44314 44488 ENSPTRG00000047549;ENSPTRG00000047549;ENSPTRT00000090571;1;+
$ 1 57173 57284 ENSPTRG00000047549;ENSPTRG00000047549;ENSPTRT00000090571;9;+
$ 1 56011 56523 ENSPTRG00000047549;ENSPTRG00000047549;ENSPTRT00000090571;7;+
$ 1 55588 55704 ENSPTRG00000047549;ENSPTRG00000047549;ENSPTRT00000090571;5;+
$ 1 56732 56857 ENSPTRG00000047549;ENSPTRG00000047549;ENSPTRT00000090571;8;+
bedtools sort -i Pan_tro_3.0.99_cds.bed \ # ソートする
> Pan_tro_3.0.99_cds_sorted.bed
bedtools merge -i Pan_tro_3.0.99_cds_sorted.bed \ # 重複領域をマージする
> Pan_tro_3.0.99_cds_sorted_merged.bed
head -6 Pan_tro_3.0.99_cds_sorted_merged.bed
$ 1 44314 44515
$ 1 45175 45210
$ 1 49376 49567
$ 1 49653 49688
$ 1 50212 50314
- GATKなどの領域ファイルとして使いたい場合は事前にindex化しておきます。
gatk --java-options "-Xmx64G -Xms64G" IndexFeatureFile \
-I Pan_tro_3.0.99_cds_sorted_merged.bed
以上。（2022年1月4日、了。）