アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
フォーマットもいろいろ
バイオインフォマティクスでは多くのファイル形式が使われますが、GFFとGTFは名前も似ていてややこしいですね。ということで、今回はGFFとGTFの違いに触れながら、フォーマットの説明をしたいと思います。

GFF(General Feature Format)はゲノムの配列に関連した特徴を示した9列からなるタブ区切りのフォーマットです。
chr1    hg19_rmsk       transcript      10001   10468   .       +       .       ID=(CCCTAA)n;geneID=(CCCTAA)n
chr1 hg19_rmsk exon 10001 10468 1504.00 + . Parent=(CCCTAA)n
chr1 hg19_rmsk transcript 10469 11447 . - . ID=TAR1;geneID=TAR1

それぞれの列の説明です。
1. seqname:染色体またはスキャフォールドの名前
2. source:データベースやプロジェクト名
3. feature:配列のタイプ(例:repeat, exon, promoter, etc)
4. start:配列の開始位置
5. end:配列の終了位置
6. score:任意のスコアが入る。ない場合は '.'
7. strand:+(forward)、-(reverse)または '.'
8. frame:0から2の数値で、翻訳を開始する塩基の位置を示す(0なら配列の1番目がコドンの1番目の塩基になる)。コーディング領域でない場合は '.'
9. attribute:グループ

GTF(Gene Transfer Format)はGFFと8列目までは同じ項目ですが、GTFでは9列目にgene id とtranscript idを含み、セミコロンとスペース区切りで連続した複数の情報が付加されています。
chr1    hg19_rmsk       exon    16777161        16777470        2147.000000     +       .       gene_id "AluSp"; transcript_id "AluSp";
chr1 hg19_rmsk exon 25165801 25166089 2626.000000 - . gene_id "AluY"; transcript_id "AluY";
chr1 hg19_rmsk exon 33553607 33554646 626.000000 + . gene_id "L2b"; transcript_id "L2b";


ここまで簡単に2つのフォーマットの説明をしましたが、もっと詳しく知りたい方は、下記のリンクをご参照ください。
UCSC FAQ
Sanger Institute
The Sequence Ontology Project
Brent Lab Homepage
| onouek | バイオインフォマティクス | 14:33 | comments(0) | - |
Python デバッグ
Pythonに限った話ではありませんが、一生懸命プログラムを書いて、いざ実行すると、いくつものバグに遭遇します。エラーが出る、あるいは、エラーは出ないのに、プログラムが意図したとおりに動いてくれないということも。

こういう場合、きっとどこかで条件文を書き間違えたとか、整数同士で割り算をしたら結果の小数点が切り捨てられて整数になってしまったとか、そんな些細なミスを犯しているに違いありません。
しかし、どこで何がおかしくなっているのか、何行にもわたる長いプログラムを眺めていても、たいていの場合は解決しません。
そんなとき助けてくれるのが、デバッガです。

以前「Perl デバッグ」という記事にて、Perlのデバッガについてご紹介しています。
今回の記事では、同じように動くPythonの標準のデバッガをご紹介します。
大体同じような内容ですので、Perlの記事と併せて読んでいただけましたら幸いです。

Pythonを実行するときに -m pdb と付けることで利用できます。
(pdbって、python debugの略でしょうか?)

$ python -m pdb sample.py
> /home/work/directory/sample.py(9)()
-> """
(Pdb) コマンドを入力します

コマンドはPerlとほとんど同じです。

list(l):実行が止まっている行のソースが前後あわせて10行表示されます。

next(n):一行実行します。関数やモジュールは飛ばします。

step(s):nextと違い、関数やモジュールの中まで呼び出して一行ずつ実行します。

print(p):「print 変数」と実行することで、実行時点での変数の中身を表示できます。変数に意図した値が入っていない、といった問題を確認できます。

continue(c):ブレークポイントまでプログラム実行を進めます。ブレークポイントが設定されていないとき(あるいは、ブレークポイントが条件文の中に設定されていて、その条件が一切満たされないとき)は、プログラムの最後まで実行します。

break(b):ブレークポイント(continueコマンドの実行を止める場所。長いコードの、目的の場所の前や、ループの中などに設定すると便利です)を設定します。ブレークポイントは複数設定することができます。

clear(cl):設定したブレークポイントを削除します。

quit(q):途中で終了します。あるいは「「ctrl+D」でも終了できます。

他にもいろんなコマンドがあるので、気になった方は調べてみてください。
Pythonを書いてみたはいいけど、なぜかうまく動かない……とお悩みの方、これでバグをつぶしていきましょう!
| kubo | バイオインフォマティクス | 15:25 | comments(0) | - |
無害なメッセージと有害なメッセージ
hatです。私は昔から自分の体調を気にしすぎてしまうところがあります。

頭が痛いといっては脳腫瘍を疑い、お腹が痛いといっては盲腸を疑い、指のささくれが膿んだといっては指が壊死するのではないかと悶々とし、しかし結果的には何でもなくてここまで元気に生きてきました。

10年ほど前でしょうか、足の親指にある黒い塊を皮膚がんではと気に病み、死んだ後見られて困るものなどを身辺整理した後、勇気を出して皮膚科を診察した結果、「大きなホクロです」と診断されたことがあります。
せっかくなので「でも変な形をしているんです」などと少し食い下がってみましたが、「ただの変な形のホクロです」ときっぱりと言われました。

そう言われると確かにそれはホクロにしか見えず、安堵するやら恥ずかしいやらの複雑な心境で帰ってきました。そして、一目で良性か悪性か判断できるとは、お医者さんはさすがプロだなあと思いました。

ソフトウェアを実行していて警告やエラーが出た時に、気にするべきものか無視していいものか判断が迷うことがあります。

GATK(v1.6)を実行していて、次のような警告が出ることがあります。
WARN 09:58:00,492 RestStorageService - Error Response: PUT '/GATK_Run_Reports/xxx.report.xml.gz' -- ResponseCode: 403, ResponseStatus: Forbidden,

GATKは実行時にログをBroadInstituteまたはAmazonクラウドへ送るようにしているのですが、なんらかの理由でログ送信が失敗した場合にこの警告が出るようです。
他にエラーが無く、結果ファイルが出力されていれば、実行自体はうまくいっているので気にしなくて大丈夫です。

ちなみにログを送らないようにするには、https://www.broadinstitute.org/gatk/request-key からGATKのキーファイルを申請し、GATK実行コマンドに「-et NO_ET -k キーファイル」オプションを付けて実行する必要があるそうです。

参考
http://gatkforums.broadinstitute.org/discussion/1250/what-is-phone-home-and-how-does-it-affect-me
http://gatkforums.broadinstitute.org/discussion/1721/warnings-and-errors-during-score-calibration

| hat | バイオインフォマティクス | 16:12 | comments(0) | - |
chimerascanで融合遺伝子を検出する
融合遺伝子検出ソフトウェアは数多くありますが、ベストなソフトがないのが現状だと思います。
TopHat-FusionやdeFuse(deFuseの記事)などが有名ですが、今回はchimerascanというソフトの使い方を紹介したいと思います。

1. アノテーションデータの準備


・UCSCからヒトリファレンスゲノム(hg19)をダウンロードして、解凍します。染色体のFASTAファイルのみ(chr1.fa〜chrY.fa )を1ファイルにマージします。

・UCSCからヒトトランスクリプトームのデータをダウンロードして、解凍します。chimerascanのダウンロードリストからダウンロードすることもできます。


2. ソフトウェアのインストール


chimerascanからソースコードをダウンロードして解凍します。

・chimerasanのディレクトリにあるsetup.pyを使ってコンパイルします。

・chimerascan_index.pyでchimerascanのインデックスを作成します。
インデックスの作成にはbowtie-buildを使うのでパスを通しておく必要があります。

これで、ソフトウェアの準備ができました。

3. ソフトウェアの実行


chimerascanのディレクトリには、ヒト前立腺癌のRNA-Seqデータ(Fastqファイル)があり、融合遺伝子 TMPRSS2-ERG が含まれています。
これを用いた場合、以下のように実行できます。
$ python /usr/local/bin/chimerascan_run.py -v --quals solexa /home/genome/hg19/chimerascan/ /usr/local/src/chimerascan-0.4.5/tests/vcap_pe_53bp/TMPRSS2-ERG_1.fq /usr/local/src/chimerascan-0.4.5/tests/vcap_pe_53bp/TMPRSS2-ERG_2.fq output_dir


outputに生成されたchimerascan.bedpeからhtmlを作成することもできます。
$ python /usr/local/src/chimerascan-0.4.5/chimerascan/tools/chimerascan_html_table.py --read-throughs -o chimeras.html chimeras.bedpe


chimerascanは実行までの準備が比較的簡単だと思うので、ぜひ試してみてください。
| onouek | 次世代シーケンサー解析 | 16:15 | comments(2) | - |
BED、VCFをスッキリと染色体番号順にソートする方法
こんにちは、あんドーナツ好きで有名な久保(kubor)です。
ブログのネタを日頃から探しているのですが、なかには、過去に先輩が取り上げていそうで、意外にも紹介していない話題が結構あります。

例えば、sortコマンド。
このコマンドは、言わずと知れた基本的なUnixコマンドですが、
社内でもあまり知られていない、素敵なオプションがあります。

このオプションがあればBEDファイルやVCFファイルを綺麗にソートすることが出来ますよ。

では、実際にBEDファイルをソートしてみます。
sort -k 1,1 target.bed
chr1 172113674 172113784 MIR199A2
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr2 37869024 37899342 CDC42EP3
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 161581163 161583014 AGPAT4-IT1
chr6 39282473 39290330 KCNK16
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
-kオプションで一列目を指定していますが並び順が変ですね。
ポジションもソートされていませんね。
まずはポジションをソートしましょう。
sort -k 1,1 -k 2,2n target.bed
chr1 172113674 172113784 MIR199A2
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
-k 2,2n オプションで2列目を数値扱いでソートしました。
ポジションはバッチリですね。
あとは、染色体です。「-V」オプションを追加します。
sort -V -k 1,1 -k 2,2n target.bed
chr1 172113674 172113784 MIR199A2
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
綺麗にソートされましたね。
文句なしです。

本来「-V」オプションは、ソフトウェアのバージョン番号をソートするためのオプションですが、染色体番号にも使用可能です。

※ 残念ながらOS Xの場合は、付属のBSD版のsortコマンドに、「-V」オプションがありませんので、GNU版のsortが入っているCoreutils(http://www.gnu.org/software/coreutils/coreutils.html)をインストールする必要があります。

ちなみに、VCFファイルも1列目が染色体番号、2列目がポジションと、BEDファイルの例と同じですので、まったく同じオプションでソート可能です。
ですが、オプションを毎回書いてると面倒なので、エイリアスしておくといいですよ。

sortBedはbedtoolsですでにあるコマンド(http://bedtools.readthedocs.org/en/latest/content/tools/sort.html)ですし、sortFineで、登録しちゃいましょう。

bashをお使いの方は、~/.bashrc
zshをお使いの方は、~/.zshrc
などへ以下を記述します。
alias sortFine="sort -V -k 1,1 -k 2,2n"
必要であれば設定ファイルを再読み込みします。
source ~/.bashrc
or
source ~/.zshrc

それでは、使ってみましょう。
sortFine target.bed
chr1 172113674 172113784 MIR199A2
chr2 37869024 37899342 CDC42EP3
chr2 149804343 149883273 KIF5C
chr2 203141153 203141241 SNORD70
chr3 14693252 14714166 CCDC174
chr3 49058050 49058142 MIR191
chr4 103715539 103749105 UBE2D3
chr6 39282473 39290330 KCNK16
chr6 161581163 161583014 AGPAT4-IT1
chr8 109455852 109499136 EMC2
chr9 100818958 100845365 NANS
chr10 63239798 63253189 TMEM26-AS1
chr12 49250915 49259653 RND1
chr15 41099273 41106767 ZFYVE19
chr15 72879557 72879654 MIR630
chr18 30252633 30352974 KLHL14
chr20 37377096 37401089 ACTR5
chr21 31581468 31584101 LINC00307
chrX 36011397 36019767 LOC101928564
chrX 55477927 55478015 MIR4536-2
うーん、これは相当素敵です。相当ファイン、sortFineですね。
| kubor | バイオインフォマティクス | 14:24 | comments(0) | - |
   1234
567891011
12131415161718
19202122232425
262728293031 
<< July 2015 >>

このページの先頭へ