ついにubuntuに移行した
これまでjupyternotebook便利すぎてwindows環境にこだわってたけど明らかにLinuxのが使いやすいしwindows対応してないパッケージ多すぎるしで結局RNAの発現解析調べるのはubuntuでやることにした。
まずBiocondaをインストールなどもろもろ環境整備してからfastqcをダウンロード。これで自由にfastqcファイルを見れるぞ！（これまではテキスト形式でみてた）
さて、バイオインフォマティクスを学んできて一か月、やっとお目当てのRNA発現量解析のところまできた。（長かった。。。）
今回はマッピングから。
ubuntuでBiocondaをインストールして、そこからhisat2を使って行う。
マッピングにはまず参照配列をインデックス化する必要がある。インデックスとはソフトウェアのマッピングプログラムから参照しやすくなるように整理・加工したもの。参照配列は必ずインデックス化する必要がある。
hisat-2を用いてインデックス化する際にはhisat2-build (ファイル名).fna (新しいファイル名).fnaでよい。
次にhisat2を用いてマッピングをする。
hisat2 (パラメータ指定) -x（参照配列） -1（入力配列） -2 （入力配列その２） -s（出力配列）
のようにマッピングする。マッピングのパラメーターには様々なものがある。パラメータ指定については
Manual | HISAT2 (daehwankimlab.github.io)
を参考にすればよい。
んでいざマッピング、と思ったら
Error, fewer reads in file specified with -2 than in file specified with -1
terminate called after throwing an instance of 'int'
Aborted (core dumped)
(ERR): hisat2-align exited with value 134
というエラー。どうやら2つ目のファイルがうまくトリミングできていないっぽい？
のでもう一度trimmomaticでトリミングしてから同じコードでマッピング。したらうまくいった。どうやらさっきやったトリミングはうまくできていないらしい。
ちなみにこのサイト覗いてみると、どうしてもread数がそろわないけどマッピングしたいときにはhisat2ではエラーが出てしまうのでSTARを使えばいいらしい。
Paired-endのデータが不揃いの場合に起こるエラーについて - Qiita
無事samファイルを出力しようと思ったのでbamファイルを作ろうと思い、
samtools view -@ 32 -b (ファイル名sam) -o (ファイル名bam)
で実行したら事件。
samtools: error while loading shared libraries: libcrypto.so.1.0.0: cannot open shared object file: No such file or directory
とまたもやエラー。Biocondaでsamtoolsをインストールするとlibcrypto.so.1.0.0を認識できずにこうしたエラーが起こるらしい。
色々ネットに解決策が転がっていたが、Ubuntuのアップデートによりopensslがアップグレードされてしまい、かなりめんどくさい状況。opensslをダウングレードすればいいっぽいが、すると他の色々もダウングレードされるので
conda uninstall samtools
でsamtoolsをアンインストールして、
sudo apt install samtools
でsamtoolsを再インストール。逃げの一手。笑。で、結局うまくbamファイルを作れました。めでたし。
samtools sort -@ 32 （もとファイルbam） > （ソートするファイル.sorted.bam）
以上のようにsortしてbamファイルの準備は完了。このソート作業でリードはマップされたリード開始位置の順になる。
またsortの際にmerging from 0 files and 32 in-memory blocks...と出るがどうやらエラーではないっぽい。sortするときにreadを一時的に保存した場所を知らせているだけらしい。
numpy np.r_/np.c_の使い方
numpyのnp.r_とnp.c_が便利そうだったのでまとめます。
下の記事の内容まんまです。
多次元配列の結合を行うオブジェクトnp.c_とnp.r_の使い方 - DeepAge
配列を結合したい際、次元数が違う時がある。そこでnp.r_を使えばちがう次元でも結合可能。
数値の文字列で表すことができます(a, b, cはそれぞれ整数)。デフォルトでは'0, 0, -1'となっています。
aは、どの軸(axis)の方向に沿って配列を結合するのかを指定します。bは、できあがる配列の次元数の最小値を指定します。cは、次元数の少ない配列の次元数の拡張を行った際、形状(shape)の表記として、どこに配列の最後の次元が置かれるべきかを指定します。
cについては、例えば (2,3) の形状(shape)を持つ二次元配列があったとします。この配列を他の(2,2,3)の形状(shape)を持つ三次元配列と結合したいとします。このとき、二次元配列は三次元に拡張しなければ結合できません。
c=-1またはc=1であれば、先頭に1が追加され、(1, 2, 3)の形状で結合することができます。c=0であれば、元の配列の形状が先頭に詰められるので、(2,3,1)となります。
上のリンクに書いてあったもの。np.r_に上の三種類の条件値を指定すれば結合の方法を指定できる。
それでnp.c_はこのうちa=-1,b=2,c=0を指定して入れる。
a=-1というのは最も低い次元のもとで結合していく。c=0は次元の違うものを結合させる際、末尾に1を追加するというもの。
ちなみにaxisとは座標軸のもの。行と列のみの行列なら列がaxis=0,行がaxis=1となる。
つまり、任意のx,yについてnp.gridmeshで行列をそろえ、np.ravel( )でリストにし、np.c_[ ]で結合すれば[x,y]の要素が入ったサイズが列数の行列を手に入れることができる。
WindowsからGraphviz利用
Graphvizを使うのに手間取ったので覚え書き。
決定木を書こうと思い、pydotplusを使いデータの読み込み、pydotplus内のgraphvizでデータの書き出しをしようとしたがここでエラー。
graphに決定木の学習データclfを入れてgraph.write_png('DecisionTree.png')でグラフを書こうと思ったがエラーが。
pydotplus.gaphviz.InvocationExceptionと表示され、pngファイルを認識できません、Use one of : でそのまま何も出ず。
graphvizはpip install graphvizでインストール済み。
これはgraphvizがただインストールしただけじゃダメな奴かな、と色々調べたところどうやらWindowsではコンパイラに何か問題があるぽい。やっぱMacほしいなあ。。。
しょうがないのでpip uninstall graphvizでgraphvizをuninstallしてgraphvizの公式サイトから直接ダウンロード。zipファイルを解凍してPathを通してリトライ。
コマンドプロンプトで dot -vと打ち込んでバージョン情報が出てくればgraphvizがインストールされている状態。無事バージョン情報が出てきたので、よっしゃこれでいけるわ！と思い再びpngファイルを書き出してみる。
エラー。
また同じエラー。嘘じゃん。
エラー内容にdot -cを入れてみたらどうですか、とあったのでコマンドプロンプトに戻りdot -cを入れる。
これが入らん。問題はdot -cが無事実行されれば解決されるだろうと思い、公式サイトのWindowsのインストール方法に飛ぶ。
New simplified installation procedure on Windows - Announcements - Graphviz
このサイトを読んでみると「管理者として実行せよ」との文字が。なるほど。コマンドプロンプトを管理者として実行してもう一度dot -cを打ち込んでみる。
Trueが返ってきた。どうやら管理者からしかアクセスできないようになっているらしい。
これでいけるやろと思いもう一度pythonからpngファイルの書き出し。ドンピシャ。無事決定木の画像データが表示された。
Graphvizは系統樹作成にも使えるから、graphvizをインストールできたのはでかい。
ちなみにdot言語を学ぶときは下のウェブサイトを参考にしよう
IGVの使い方
IGVの使い方について少し。
pythonのigv Jupyter Extensionを使えばどうやらpythonでも使えるらしいけど、リファレンスゲノムが多分サーバー上にあるやつしか使えないんじゃないかな。ということでサーバー上にない生物種を使うために結局アプリをインストール。
ちなみにigvでbamファイルを扱うためにはインデックス化してbam.baiの形にすることが必要です。
ここからはアプリでリファレンスゲノムを設定するための覚え書き。以下のサイト参照。リファレンスゲノムの読み込み | IGV ゲノムブラウザ (biopapyrus.jp)
要はリファレンスゲノムのFastaファイルを用意し、Genomeの欄からCreate.fileを選択してリファレンスゲノムをインポートすれば用意すればいいらしい。
OptionでGFFファイルを打ち込めるらしいが、必須ではないだろう。これは遺伝子名を入れたいときに必要なのか。
BioPythonまとめ その2
その1が長くなったのでその２です。
～NGS処理ツール～
まずはおなじみSAM・BAMファイルからの抽出をpythonでもやってしまおうというもの。できるんですね。Linuxでやると調べるのに苦労するからpythonで実装してくれていてとてもありがたい。ちなみにこいつらもともとC⁺⁺言語用だったようで、これらを覚えているとサクッといけたらしい。pythonでも実装しようということで、python処理用の薄い皮で包み込んだpysamを誰かが作ってくれたのです。ありがたい。
と思ったら事件。Windows10ではpysam使えないらしい。インストールできない。はあ。。。
仮想Linux環境作っておいてほんと良かった。。。
これからSAMファイルとかBAMファイル使うときにpythonいじるときはこっちからになりそうだな。。
一応BAMファイルのパーサーとしてbamnosticがWindowsに対応しているらしいが。。。本当かな？
bamnosticはpipから無事インストールできました！使い方は下参照。
Welcome to bamnostic’s documentation! — bamnostic 1.1.4 documentation
これでSAMファイルもpythonからいじれるようになる。。。かな？
一応ubuntuでpysamもやっとこう。ダウンロードの仕方などはこちら。
pysam does not install on windows · Issue #575 · pysam-developers/pysam · GitHub
ubutuでpysam自体はインストールできたが、bamfile=pysam.AlignmentFile('bamファイル名’,'rb')とやっても、could not retrieve index file for 'ファイル名'と出てきてファイルの読み出しができない。普通こういう時はパスが通ってない時だが、そもそもbioinfoのディレクトリにいるのだからパスもくそもない。。。
むむむ。。。ファイルが大きすぎるのだろうか。import osからos.listdir()で現在のディレクトリ内のファイルを確認しても目的のファイルは列挙されるし。。。
わからん！いかんせん仮想環境内だから対策も練りようがない。
ここは大人しくbamnosticにいこう！さらばpysam！君には期待していたのだが！
ちなみに、覚書なのだが、samtoolsはもともとC/C⁺⁺ベースで作られたツールのため、samtoolsをいじる際に出てくる-oとか-uとか、こういうのは全部C/C⁺⁺にしたがっているぽい。今後調べるときの参考に。
お次はbamnostic。コマンドプロンプトでうまくできなかったので、一度諦めてjupyternotebookにアップロードしてから動作を実行した。パスは一応デスクトップに通してあるはずなのだが、なぜだろうか。とりあえず今沼にはまりかけているので、できるところからやる。
bamnosticの使い方はpysamと大体同じっぽい？データの読み出しの時に、ランダムアクセスはできませんうんちゃらかんちゃらとの警告が出たが、無視してそのままヘッダーを表示。(.header)よし、うまくいった。
ubuntuでできなかったのも、ランダムアクセスが制限されているからなのでは？と考えて、ubuntu でもやってみた。一応。
すると。。。なんだ、アクセスできるじゃん！
流れとしては、pysamを使ってAlignmentFileからbamファイルを開いてみたが、could not retrieveとかでファイルにはアクセスできない風だった。
bamfile = pysam.AlignmentFile('unmapped.bam','rb')
が、次の行で
print(bamfile.header)
を使うと下の行にズラーッとヘッダーが表示されてきた。つまり仮想Linux環境下で、カレントディレクトリのbamファイルにアクセスできた、というわけである！
ここからは推測だが、AlignmentFileからアクセスできなかったのは、たぶんbamnositcと同様にランダムアクセスをしてしまっていたから。headerなど指定してアクセスする分にはうまくいく。ということで今後はAlignmentFileを打ち込んだのちに起こるエラーは無視することにしよう。
またここでヘッダーとは、リファレンスとは、などと疑問に思ったので、SAMファイルの構造から理解することにした。
参考となるURLは↓
crusade1096.web.fc2.com/sam.html
要はSAMファイルには下のように11列の配列情報が入っているということ。
1列目：リードの名前
2列目：リードの状況（どんな風にマッピングされてるか判断できる，以下参照）
3列目：張り付いた染色体，コンティグの名前
4列目：張り付いた場所
5列目：マッピングスコア
6列目：マッピング状況（indel入ってるとか，どのくらいマッチしてるとか判断できる，以下参照）
7列目：paired-endの時の相方の名前
8列目：paired-endの時の相方の場所
9列目：paired-endの時のインサートの長さ
10列目：リードのシークエンス配列
11列目：リードのクオリティ
さて、構わずbamnosticやpysamで処理を続けていたが、ここで問題が生じてきた。ヘッダー情報やレングス情報など、ただファイルから情報を取り出すだけの処理ならば問題なくできるのだが、例えばpileupやsortなど、がうまくできないなあ。
調べてみるとbamnositicにはそもそもsort機能がついてない。sort以外の機能、読み出しなどについてはpysamと同じ感じなのか。ただbamnositicはランダムアクセスをするとエラーを返してくるようなので、for構文を使ってうまく順々にファイルアクセスするようにしてあげなければならないぽいなあ。。。bamnosticはNGS処理では使い勝手が悪いか。。。
お次pysam。これはまだ未知数だが、windowsでは仮想環境上に構築しているため、いかんせん使い勝手があまりよくないのでござる。。。
ただし、sortとindexは成功したようだ。
unmapped.bamに対して、pysam.sort('-o','unmapped.sorted.bam','unmapped.bam')を実行したらsortファイルであるunmapped.sorted.bamがカレントディレクトリ上に出現した。すなわち成功。さらにpysam.index(sortファイル)を行ってもエラーが出現しなかったため、これも成功ということでよかろう。
ちなみにindexが無事終わると、.bam.baiファイルが出現する。.baiファイルを作ることをbamファイルにindexを付与すると言うのだろうか。このbaiファイルはマッピングデータを視覚的に表示できるIGVファイルで用いるときに必要なようである。他の用途はわからん。。。下のサイトに色々とIGVの使い方が書かれている。
ただしpysamで問題なのはfetch(出現させる)である。一度やってダメだったが、indexを付与したファイルを用いて、
samfile = pysam.AlignmentFile('unmapped.sorted.bam','rb')
を打ち込んでみる。すると。。。
なんと何のエラーもなく実行された！！！
これらのことから、普通にbamファイルを上のコードを起動するとエラーが起こるのは、indexが付与されていなかったからであり、index付与後には問題なく扱えることがわかった。
bamnosticよりもpysamの方がよさそうである。これから実際の解析にはこれを使おう。（先ほどbamnosticにはソート機能がないと書いたが、bam.baiにできるところがありそうなので、詳しく調べればソートできるところもあるのだろうか？）
さて、pysamをだいぶ使いこなしてきました。ソート後インデックスを付与したファイルを対象（今回はunmapped.sorted.bam）に、
countread = samfile.count()
print(countread)
と打ち込むとカウント数がでた！これでBAMファイルにいくらリード数が入っているかわかるね！ついでにcountの()のなかにcontigや範囲を指定すればそこにおいてマッピングされたリードの数を数えられます。
とりあえずこれでpysamは終わり。samtoolsのviewは使えるぽいけど、bamファイル同士を結合させるmergeはないのか？そしたらsamtools使用に逆戻りだなあ。。。
と思ったらpysam.merge(""all.bam"",''入れたいファイル１","入れたいファイル２",,,,,,)とやれば入るっぽい。
unmappedも取り出せるらしいし、先生に言われた処理全部パイソンでできそう。多分表記方法さえ変えれば全部pysamで代用できるっぽいな。すごいぞpysam！今日一日かけた甲斐があった。pileupは必要なときにまた。
GO解析に移ったが、正直何がなにやら
とりあえずgoatoolsをインストールすればよいのかしら。GO解析の入出力は主にエクセル形式で行われるぽい？
また、jupyternotebookにサンプルコードを残しているわけだが、まあ必要になったら参照する程度でよいかと、、、かなりムズイ。
一つ収穫は今回のように、pydotなど、graphvizを使ったものはdot.exeという実行ファイルにパスが通ってないと使えないということ。今回もjupyternotebookで実行するとうまくパスが通らずに'File Not Found Error'、'dot' not foundと表示された。
jupyternotebookでは結局うまくできなかった。コマンドプロンプトでは実行してないが、dot.exeにもパスは通してあるのでうまくできる。はず。これにこだわってもしょうがないので一回スルーしよう。
getattr使い方
もう完全にpythonのことしか書いてない・・・（笑）
getattrで苦戦したので書き残しておく。
問題のコードは下。recordはオブジェクト。
for u in dir (record):
if u[0] != '__' :
print(u,getattr(record,u))
これが意味わからんかった。record内の属性をdir及びforで一つずつ書き出し、__()__の形ではないものを（!=はノットイコールの意味）プリントする。プリントの際に問題となるのがgetattr。教科書的な説明では
object の指名された属性の値を返します。 name は文字列でなくてはなりません。文字列がオブジェクトの属性の一つの名前であった場合、戻り値はその属性の値になります。例えば、 getattr(x, 'foobar') は x.foobar と等価です。指名された属性が存在しない場合、 default が与えられていればそれが返され、そうでない場合には AttributeError が送出されます。
というもの。いろんなサイトを調べてよくわからなかったが、結局公式サイトの説明が一番わかりやすかった。getattrは指定されたオブジェクト及び属性の値を返すようだ。問題となったサンプルコードでは、
① for構文で属性を一つ取り出す
②printとgetattrにより、取り出した属性の値を文字列で返す。このときにu,が入っているため、属性の値だけでなく属性も返して見やすくしているようである。
とまあ、こんな感じ。理解すると簡単だけど理解するまでが難しいっす
Biopythonまとめ その１
日記風に書いていたが、BioPythonの難易度がpython初心者にとって高いので、これはまとめ用の記事。備忘録。じぶんのレベルが上がるにつれて適宜更新していきます。
ちなみに参考した本はオーム社の「Pythonによるバイオデータ・解析入門」山内長承著である。
step1.
重要なことはファイルの読み出し及びファイル作成であろうため、そこだけまとめる。用いるクラスはSeqRecordクラス。これでNCBIで使えるGenBank形式、ENAで使えるEMBL形式のファイルから読み書きできる。
SeqRecordクラスのもつアトリビュート（SeqRecordが持つ情報、SeqRecordオブジェクトにドットをつけてこの言葉を書けばオッケー、というもの）は配列自体（.seq）、内容の記述（.description）、アノテーション情報（.annotations→辞書型で保管されている）、フィーチャー情報（.features、SeqFeatureクラスで書かれる）、データベースへのクロスリファレンス（.dbxrefs）などがある。これらはGenBank形式やEMBL形式と一致している。
まずファイルの読み出しについて。
ファイルの読み込みはSeqIOクラス（データの入出力担当！）のSeqIO.read()によって行う。パラメータにファイル形式を指定することで、その形式に沿って情報を分解して読み取り、SeqRecordオブジェクトとして返す。
使用の際はrecord=SeqIO.read("ファイル名”,"ファイル形式（genbankなど）")のように用いる。画面に表示したいならprint(record)でなんとかなるが、print(repr(record))にしたほうが、より見やすいだろう（reprとは私たち操作する人間にみやすいようにわざわざ変換せず、処理の途中、ありのままの形式で表示する）
ちなみSeqIO.readオブジェクトの後ろにさらにSeqFeatureクラスのオブジェクトのようなもの（この書き方はたぶんよくない）が続いている。SeqFeatureクラスの中にはさらにtype,location,qualifierなど、配列データの細かい情報が続いている。これはリストの形式で保管されている。またmypositionに参照したい塩基の番号を入れ、if myposition in feature:の形にすれば比較したい位置がfeatureの範囲内にあるかどうかを判定してくれる。
次にファイルの書き出しについて。特にFASTAで書き出すときについて。SeqRecordクラスのコンストラクタを使う、らしいが、いまいちこのコンストラクタを使うというのがピンとこない。。。まあ雰囲気はわかるので、やり方だけ。
rec=SeqRecord(Seq("配列"),id="id",description="description")というように書いていけば何とかなる。うん。これでできた複数のrecをリストに入れれば、複数の配列を持った一つのFASTA配列を作ることが可能。イメージとしては、FASTAファイル作成に必要な情報をSecRecordの形で入れたものがrec、ということで、これをSeqIO.write("作ったオブジェクト","新しく作りたいファイル名","fasta")とすれば完成。
Step2
コマンドプロンプトからpythonを立ち上げ、clustalWを呼び出してアラインメントさせる及び樹形図の作成というもの。
まずClustalWをダウンロードしてPathをとおさねばならなかったが、これでだいぶ手こずった。。。ClustalWのダウンロードされた場所が不明。パスの通しかたも不明。
最終的にはパソコンのハードディスク上のフォルダの一つにあったので、（Program Files（×84）みたいなやつ）そのパスをコピーしてパス通してもよかったが、初心者のうちからパスの変更はあまりしたくないので、clustalw2.exeをデスクトップに移動してデスクトップにパスを通して操作完了。完全にやり方初心者。
上手くできたと思ったが樹形図が書けない。FASTAファイルに入っていた配列が一つだったからか？今更FASTAファイル作るのも面倒だし、アラインメントなら普通にウェブ上からできるので、ここは一旦飛ばして次に行きます。（英断）
step3
NCBIデータベースへのアクセス方法。
まずはEinfoについて。Einfoはまあ大体わかったが（そして多分つかう機会はないが）、ひとつprint関数についてよくわからないことがあったのでここではそれをまとめる。
普段よく使うprint関数。普通に使う分にはいいが、たまに中に％が入る形がでてくる。意味わからん。
調べたら色々でてきたのでまとめる。
まずsep=""。これはprintの()のなかに入れると区切り文字を変更することができる。例えば
# デフォルトの設定
print("Blue", "Red", "Green")
>> Blue Red Green
# 区切り文字を "+" に変更
print("Blue", "Red", "Green", sep="+")
>> Blue+Red+Green
# 区切り文字を無くして続けて出力
print("Blue", "Red", "Green", sep="")
BlueRedGreen
て感じ。
またend=''をprintの()の最後に加えると、出力後に自動で改行されていたのが、されなくなる。これはデフォルトではend='/n'になっているため、自動で改行される仕組み。だからend='/n'の/nの前に文字を入れれば、最後に文字を追加することも可能である。
またファイルをもともと書き出しモードで開いておき、printの最後にfile=ファイル名を入れればテキストファイルの書き出しが可能である。以下のようになる。
myfile = open("output.txt", "w")
print("Hello", file=myfile)
myfile.close()
さてこの後に本題の％について。
Python | 書式化演算子%を使った文字列の書式設定(printf形式の書式化) (javadrive.jp)
を参考にしました。てかこれ見れば解決。％は書式化演算子と呼ばれるもので、ひと昔前には主流の演算子だったよう。最近はformatメソッドが一般的になり、あまり使われていない模様。説明は参照先を読めば十分なので、今回は問題のコードについて。
for field in result['DbInfo']['FieldList']:
print('%(Name)s,%(FullName)s,%(Description)s' % field)
これを説明。％(Name)s、など%（）は、ここに何かを代入しますよ、ということ。（）の中身は引数。今回printの最後に%fieldがある。fieldは数あるFieldListの値の中の一つであり、辞書型になっている。辞書型のうちの一つを取り出してprintにしたいときに％は有用なよう。Name,FullName,Descriptionなどのkeyの値を取り出したいときに％の横に書けばよい、ということである。
step4
～NGS処理に関わるpythonライブラリ～
アノテーションなどのデータへのアクセスサポートから。まずEnsenmlがでてきたが、これはそもそも脊椎動物に特化したアノテーションのデータのサイトなうえに、現在もpythoでの使い勝手がよくないようなのでとばします。。。一応pyensemblrestをpipからインストールすれば使えるようになる、ぽいけど使う機会も来ないだろう。
そこでBioServicesから学ぶこととする。
BIOSERVICES: access to biological web services programmatically — bioservices 1.7.10 documentation
このURLを参照のこと。
まずはUniprotから。
Uniprotは主にタンパク質の配列やアノテーションデータをまとめたもの。
UniProtの初期設定としてu=UniProt(verbose=False)としてから色々はじめたほうがいいぽいね。
最初にUniProt内のIDがわかってる場合はu.retrieveを使って検索をかける。frmt="xml"
とすればxml形式で返してくれる。辞書型でaccession(IDのこと)、protein、gene、sequenceなどをキーとして返してくれるので値を取り出せばよい。sequenceの値の4番目に配列データが入っているので、指定の際にprint((オブジェクト['sequence'])[3].text)とすれば配列データをテキスト形式で返してくれる。
またu.retrieveにはXML形式を解析するために（XMLパーサー）にElementTreeライブラリ、HTMLを解析するために（HTMLパーサー）beautifusoup4が組み込まれている。そのため、それらのメソッドやアトリビュート、.findAll()や.root、.root.tag、.root.attribを用いることができる。
※Webスクレイピングについて
主にrequestsとbeautifusoup4を用いる。requests.get(url)でデータの獲得、beautifulsoup4でデータの解析を行う。その後、findallを用いて欲しいタグをすべて検索できる。
更にbioservices.KEGGモジュールを使えばKEGGサービスへのアクセスが可能である。KEGGとはKyoto Encyclopedia of Genes and Genomesのことである。Wikipediaによると、
KEGG（Kyoto Encyclopedia of Genes and Genomes：日本語では「京都遺伝子ゲノム百科事典」の意味）は、遺伝子やタンパク質、代謝、シグナル伝達などの分子間ネットワークに関する情報を統合したデータベースであり、バイオインフォマティクス研究に利用される。このデータベースは、細胞レベルでの生命システムの機能に関する知識を、分子間相互作用ネットワーク（代謝、シグナル伝達、遺伝情報等）の二項関係に基づいた情報としてデータベース化し (PATHWAY)、これを中心に据えているのが特徴である。さらに遺伝子カタログ情報 (GENES)、既知のタンパク質間の配列相同性情報 (SSDB)、機能的類似性情報 (KO)、生体関連化学物質に関する情報 (LIGAND) などに関する各データベースを統合している。
生物種としては、ヒトのほか動物・微生物を中心とした各種モデル生物が対象とされており、種による異同を調べたり、ネットワークの2要素の間で可能な経路を計算するなどが可能となっている。
ということである。bioservicesからKEGGをimportし、そののちKEGGからorganismIdを呼び出す。モデル生物一つにIdが割り当てられており、例えばヒトはhsaである。あとは.getをつなげればgeneの情報やpathwayの情報も得られる。pathwayは画像保存されており、get("pathway番号","image")としてダウンロードし、下のようにすれば画面に表示できる。
f = open("hsa05415.png","wb")
f.write(res)
f.close()
import cv2
img = cv2.imread("hsa05415.png")
cv2.imshow("hsa05415",img)
cv2.waitKey(0)
cv2.destroyAllWindows()
これもbioservicesからアクセスできる。bioservicesすごい。。。
from bioservices import QuickGO
でQuickGOを捕まえ、g=QuickGO(verbose=True)にする。なんでここでTrueにするかはよくわからんが、とりあえずTrueにしといたほうが良いっぽい。
参考文献の情報とは違うのだが、最新のコードは10. Utilities — bioservices 1.7.10 documentationの通り、g.get_go_terms("id")でデータを出力する。
Annotation部分は参考文献と変わらず。g.Annotationの後にtaxonIdとwithFromで情報を検索する。
ちなみにこうした辞書型で出てきたものの清書の仕方は、
for 変数 in (出てきたオブジェクト).keys():
print(変数,(オブジェクト)[変数])
のようにすれば可能である。
