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(変数,(オブジェクト)[変数])

のようにすれば可能である。