Biopythonを使ってみる

俺は都会のマンションのベランダで、自分勝手に植物を育てる、ベランダーだ。
季節は春。俺は久しぶりに、ブログを更新することにした。
(植物男子ベランダー・一部改)

はてダがサービス終了ということで、はてブに移ってきました。

そもそもここにいろいろ書くのも最近はサボっていたので、まぁ、これを機会にまたちょこちょこ書くようにしていければとは思うんですけどねー(言うのは簡単)。

さてさて、ふとしたことでライフサイエンス系公共データをいじることとなり(ってまぁ、それが本業といえば本業なんですが)、いじるというか知識抽出なんですが、今までBioPerlでちゃらちゃらっとやっていたのを、せっかくだからBiopythonを使おうかと思い立ったわけです。

もともとpythonは10数年前から知ってはいてread-onlyだったです。今まで何回か、書きかけたこともあったですけどね。まぁ自分のころはperlが盛んにこの分野で使われていましたからねぇ。流行の最末端の自分ですが、世の中、だいぶpythonが隆盛してきたので(ってこともないですが)こういう機会でもないとpython使わねーしなというのもあり。

とまぁ、前置きが長いですが、
qiita.com
とか
biopython.org
を参考に書いてみたりして。

やったことはGenBankファイルのparse

元ファイル

LOCUS       EU919294                1220 bp    DNA     linear   INV 26-JUL-2016
DEFINITION  Lampides boeticus voucher RMBR:102827 cytochrome c oxidase subunit
            I (COI) gene, partial cds; mitochondrial.
ACCESSION   EU919294
VERSION     EU919294.1
KEYWORDS    .
SOURCE      mitochondrion Lampides boeticus (pea blue)
  ORGANISM  Lampides boeticus
            Eukaryota; Metazoa; Ecdysozoa; Arthropoda; Hexapoda; Insecta;
            Pterygota; Neoptera; Holometabola; Lepidoptera; Glossata; Ditrysia;
            Papilionoidea; Lycaenidae; Polyommatinae; Lampides.
REFERENCE   1  (bases 1 to 1220)
  AUTHORS   Lohman,D.J., Peggie,D., Pierce,N.E. and Meier,R.
  TITLE     Phylogeography and genetic diversity of a widespread Old World
            butterfly, Lampides boeticus (Lepidoptera: Lycaenidae)
  JOURNAL   BMC Evol. Biol. 8, 301 (2008)
   PUBMED   18973689
  REMARK    Publication Status: Online-Only
REFERENCE   2  (bases 1 to 1220)
  AUTHORS   Lohman,D.J., Peggie,D., Pierce,N.E. and Meier,R.
  TITLE     Direct Submission
  JOURNAL   Submitted (23-JUL-2008) Department of Biological Sciences, National
            University of Singapore, 14 Science Drive 4, Singapore 117543,
            Republic of Singapore
FEATURES             Location/Qualifiers
     source          1..1220
                     /organism="Lampides boeticus"
                     /organelle="mitochondrion"
                     /mol_type="genomic DNA"
...

書いたもの

#!/usr/env/python

import sys
import os
import inspect
from Bio import SeqIO


gb_file = sys.argv[1]     # input gb file path
for record in SeqIO.parse(gb_file, 'genbank'):

    id = record.id
    name = record.name
    desc = record.description
    xrefs = record.dbxrefs     # list
    annots = record.annotations     # dic
    features = record.features # list
    seq  = record.seq

    print('ID:\t' + id)
    print('name:\t' + name)
    print('desc:\t' + desc)
    print('seq:\t' + seq)

    # annotation
    ## source
    source = annots['source']
    organism = annots['organism']
    taxonomy = annots['taxonomy']   # list

    print('source:\t' + source)
    print('organism:\t' + organism)

    print('taxonomy:')
    for t in taxonomy:
        print('\t' + t)

    ## reference
    for ref in annots['references']:
        pmid = ref.pubmed_id
        medline = ref.medline_id
        journal = ref.journal

        print('pmid:\t' + pmid)
        print('medline:\t' + medline)
        print('journal:\t' + journal)

    # features
    for feat in features:
        type = feat.type
        print('type:\t' + type)

        for k in feat.qualifiers.keys():
            print('subtype:\t' + k)
            vals = feat.qualifiers[k]
            for v in vals:
                print('\t' + v)

 
で、出てきたもの

ID:     EU919294.1
name:   EU919294
desc:   Lampides boeticus voucher RMBR:102827 cytochrome c oxidase subunit I (COI) gene, partial cds; mitochondrial
seq:    AACATTATATTTTATTTTTGGAATTTGAGCAGGAATATTGGGAACATCATTAAGAATTTTAATTCGAATAGAATTAGGAACTCCAGGATCTTTAATTGGTGATGATCAAATTTATAATACTATTGTAACAGCTCATGCTTTTATTATAATTTTTTTTATAGTTATACCTATTATAATTGGTGGATTTGGAAACTGATTAATCC
...
TTTCAGTTGATTAGCTACTATTTATGGAACTCAAATTAACTATAGACCTTCAATATTATGAAGTTTAGGTTTCATTTTTTTATTTACAGTAGGCGGATTAACAGGAGTAATTTTAGCTAACTCTTCAATTGATATTACTTTACATGATACTTATTATGTTGTAGCACATTTCCATTATGTTTTATCTATAGGAGCTGTATTTGCAATTTTTGGAGGATTTATTCATTGATATCCTTTATTTACTGGATTAATAATAAATCCTTATCTATTAAAAATTCAATTTATTATTATATTTATTGGAGTTAATT
source: mitochondrion Lampides boeticus (pea blue)
organism:       Lampides boeticus
taxonomy:
        Eukaryota
        Metazoa
        Ecdysozoa
        Arthropoda
        Hexapoda
        Insecta
        Pterygota
        Neoptera
        Holometabola
        Lepidoptera
        Glossata
        Ditrysia
        Papilionoidea
        Lycaenidae
        Polyommatinae
        Lampides
pmid:   18973689
medline:        
journal:        BMC Evol. Biol. 8, 301 (2008)
...
type:   source
subtype:        organism
        Lampides boeticus
subtype:        organelle
        mitochondrion
subtype:        mol_type
        genomic DNA
...

 
きったねーコードというのは承知しているんですけどね。
まぁ、ボチボチ手直ししつつ、さらにこれを使って先に進めれれば。

PlatanusをMacで使う

縁あって、ゲノム解析とか発現解析をしています。
de novoでゲノムデータをassembleするというので、Platanus ( http://platanus.bio.titech.ac.jp/ ) というのを使ってみることとしました。現在のバージョンは1.2.4のようで(去年に出たようだ)。サイトでは1.2.1ってのも配られていて、こちらはMacのバイナリもあるようですが、最新版の1.2.4はLinux 64 bit binary (precompiled)かソースで、という状況。
そこいらのマシンで自分のMacProが一番性能がよさげなので、これで動かすべくソースを落としてみました。
一応、trimmerも使おうかと思うので、まずはこちらから。ほどきまして、makeしますと、

nakazato@grouper:~/Downloads/Platanus_trim_v1.0.7$ make
g++ -o common.o -c common.cpp -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -DRUN_MODE=\"\"
clang: warning: argument unused during compilation: '-s'
In file included from common.cpp:2:
./common.h:13:10: fatal error: 'omp.h' file not found
#include
^
1 error generated.
make: *** [common.o] Error 1

はい、怒られました。もろもろ調べたところ、Macのg++はclangの方で(イミフ。でもよし)OpenMPは使わないようになっているとか。そこで(brewbrew install gcc49してからか?)Makefile

CXX = g++
 ↓
CXX = g++-4.9

しますと

nakazato@grouper:~/Downloads/Platanus_trim_v1.0.7$ make
g++-4.9 -o common.o -c common.cpp -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -DRUN_MODE=\"\"
g++-4.9 -o trim.o -c trim.cpp -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -DRUN_MODE=\"\"
g++-4.9 -o main.o -c main.cpp -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -DRUN_MODE=\"\"
g++-4.9 -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -o platanus_trim common.o trim.o main.o
ld: warning: option -s is obsolete and being ignored
g++-4.9 -o main.oo -c main.cpp -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -DRUN_MODE=\"internal\"
g++-4.9 -s -std=c++0x -O3 -funroll-loops -fomit-frame-pointer -fopenmp -o platanus_internal_trim common.o trim.o main.oo
ld: warning: option -s is obsolete and being ignored

なにやら怒られましたが、platanus_internal_trimとplatanus_trimができたのでよしとしましょう。(ちゃんと動きました)。
はてさて、今度はPlatanus本体ですが

nakazato@grouper:~/Downloads/Platanus_v1.2.4$ make
g++ -o main.o -c main.cpp -std=c++0x -O3 -funroll-loops -Wall -fopenmp -finline-limit-50000 -lm -Dnullptr=0
clang: error: unknown argument: '-finline-limit-50000'
clang: warning: -lm: 'linker' input unused
make: *** [main.o] Error 1

はい、やっぱり怒られました。で、trimmerのと同様に直したわけですが

nakazato@grouper:~/Downloads/Platanus_v1.2.4$ make
g++ -o main.o -c main.cpp -std=c++0x -O3 -funroll-loops -Wall -fopenmp -finline-limit-50000 -lm -Dnullptr=0
clang: error: unknown argument: '-finline-limit-50000'
clang: warning: -lm: 'linker' input unused
make: *** [main.o] Error 1
nakazato@grouper:~/Downloads/Platanus_v1.2.4$ cp Makefile.org Makefile.arrange
nakazato@grouper:~/Downloads/Platanus_v1.2.4$ vi Makefile.arrange
nakazato@grouper:~/Downloads/Platanus_v1.2.4$ make clean
rm -f platanus main.o assemble.o scaffold.o scaffoldGraph.o gapClose.o common.o baseCommand.o seqlib.o mapper.o gapCloseOLC.o
nakazato@grouper:~/Downloads/Platanus_v1.2.4$ make -f Makefile.arrange
g++-4.9 -o main.o -c main.cpp -std=c++0x -O3 -funroll-loops -Wall -fopenmp -finline-limit-50000 -lm -Dnullptr=0
In file included from assemble.h:25:0,
from main.cpp:22:
counter.h:455:263: error: redeclaration of 'void Counter::countKmerOrWriteTemporary(bool&, const typename KMER::keyType&, DoubleHash*, FILE*, omp_lock_t*, const KMER&, unsigned int)' may not have default arguments [-fpermissive]
inline void Counter::countKmerOrWriteTemporary(bool &loopFlag, const typename KMER::keyType &key, DoubleHash tmpOccurrenceTable, FILE *unmappedFP, omp_lock_t lock, const KMER &kmer, const unsigned iterateTimes=32)
^
counter.h:684:190: error: redeclaration of 'void Counter::countKmerOrWriteTemporary(bool&, const typename KMER::keyType&, FILE*, omp_lock_t*, const KMER&, unsigned int)' may not have default arguments [-fpermissive]
inline void Counter::countKmerOrWriteTemporary(bool &loopFlag, const typename KMER::keyType &key, FILE *unmappedFP, omp_lock_t lock[], const KMER &kmer, const unsigned iterateTimes=32)
^
In file included from assemble.h:26:0,
from main.cpp:22:
graph.h:1764:80: error: redeclaration of 'void BruijnGraph::cutBranchIterative(long long unsigned int)' may not have default arguments [-fpermissive]
void BruijnGraph::cutBranchIterative(const unsigned long long numThread=1)
^
make: *** [main.o] Error 1

やっぱり怒られました。で、うねうね調べまして、むりくり先に進めるみたいなオプションをつけるとうまくいくとかで、上記のに加えて

CXXFLAGS = -std=c++0x -O3 -funroll-loops -Wall -fopenmp -finline-limit-50000 -lm -Dnullptr=0

CXXFLAGS = -std=c++0x -O3 -funroll-loops -Wall -fopenmp -finline-limit-50000 -lm -Dnullptr=0 -fpermissive

としますと、まぁやっぱりくどくど怒られるんですが、無事にplatanusの実行ファイルができまして、ちゃんととりあえず、動き始めたんで一安心みたいなところですかね。
週末を越えて動かしているんですが、ちょっとファイルが大きかったのか、スレッド数とかメモリが小さかったのか、"WARNING: Sorry, memory exceeds specified value!!"と怒られつつ、まだ計算しております。。。(そのあたりはまた後日?

MacでのJavaのバージョンアップ

Javaって、まぁ、セキュリティのごにょごにょもあって最新版にせーとか何とか言われたりするわけでありますが、今般、8でないと動かないソフトを入れるにあたり、6から8にバージョンアップしようとしたわけであります。
そのソフト、ターミナルで入れるので、ならJavaのバージョンもターミナルで確認しますわな。

nakazato@gardeneel:~$ java -version
java version "1.6.X_XXX"
Java(TM) SE Runtime Environment (build 1.6.X_XXX-bXX)
Java HotSpot(TM) 64-Bit Server VM (build XX.XXX-bXX, mixed mode)

今はもう8になっているのでちとごまかしてますが。で、ググると、システム環境設定の中にJavaの項目があって、そこからアップデートできるというわけです。確かにそうなっていて、アップデートを促され、ポチッとやることで6から8にアップデートされたかのようになってめでたしめでたし、とは行かなくてですね。上のコマンドをたたいてもバージョンは6のままでございまして。なんかインストールするところが変わったのか、一部しかインストールされていないのか。
で、またごにょごにょ調べまして、インストールされているところにシンボリックリンクをはるかどうかの寸前まで行ったわけですが、結局、システム環境でインストールされるのはJREの方でございまして、コマンドラインの方で見てるのはJDKである、ということがわかりまして、普通にJavaのダウンロードページJava SE - Downloads | Oracle Technology Network | Oracleに行って、JDKのを落としてきてインストールしたらば

nakazato@gardeneel:~$ java -version
java version "1.8.0_101"
Java(TM) SE Runtime Environment (build 1.8.0_101-b13)
Java HotSpot(TM) 64-Bit Server VM (build 25.101-b13, mixed mode)

無事にアップデートできましたとさ。

サンタ

毎度おなじみ BioHackathon で山形は あつみ温泉に来ている。
外国からも交えて100人くらいか。今年は大所帯。
夜にNCBIとうちのと話していたらば、そのうちのが今日は何をやっていたかと訊かれていて、リストを作っていたと返事をしたところ、チェックは2回したか、お前はサンタか、ということを真顔で言われたわけである。実はこれは以下をふまえたジョークである。今度どこかで使いたい。

Santa Claus is Coming to Town(サンタが町にやってきた)

He's making a list
And checking it twice
Gonna find out Who's naughty and nice
Santa Claus is coming to town

昆虫学会+応動昆@大阪

来てます。
今回はポスター発表にしてます。内容は、まぁ、いつもどおり。本当はもっとモノに根ざした発表をしたいんですが。
せっかくなので、公共NGS中の昆虫データについて調べてみました。

↑ ポスターの図より。数字はプロジェクト数。1つのプロジェクトで複数の生物種が出てきた時は複数にカウントされているから、厳密にはこの表現はおかしいのかとも思う。
ハエが多いのはDrosophilaとかその近縁とかカが含まれているから、と思われ。
チョウ目はカイコとタバコスズメガか?、セミ目にはウンカやアブラムシが入っているから。甲虫が少ないですね、という話もあったが、生活環が長いから研究には不向きなのだ、という別の指摘が入った。
生物種は何種類か、という質問を受けたのでその場で調べてみたところ、1256種類だった。
上位の方。

770 7227 Drosophila melanogaster
83 7165 Anopheles gambiae
51 7460 Apis mellifera
40 7091 Bombyx mori
37 7159 Aedes aegypti
36 7240 Drosophila simulans
24 76194 Papilio polytes
24 7160 Aedes albopictus
23 7237 Drosophila pseudoobscura
21 7245 Drosophila yakuba
21 7029 Acyrthosiphon pisum
19 7173 Anopheles arabiensis
19 7038 Bemisia tabaci
15 7244 Drosophila virilis
15 7070 Tribolium castaneum

なんか、上の図と合わないんだが、なんで合わないかは、そのtaxonomy ID決め打ちか、下も見ているかの違いだったかと記憶している。(グラフのDrosophila melanogasterはDrosophila melanogaster なんちゃら/なんちゃらグループみたいのも含んでいる、ということ)
だいぶdryをやっていると忘れそうな感覚だが、この学会に来ている人はこういう学名を見てだいたいあれね、というのがわかる。(まぁ、自分も有名どころとか日ごろ図鑑を見つめている蝶関係とかならだいたいわかるが)。それは、wetの研究者が遺伝子発現で発現量が変動した遺伝子上位20とか見せられて、あー、あの遺伝子ね、とわかるのと同じことである。
oral聞いていると、なんちゃらアザミウマとかにわかに姿が浮かばない昆虫が乱発してニッチェ感が満載である。それに対して、私のやっているなんちゃらハマキガでは、とか細かいどうしで質問していてむしろこの人たちの研究対象をこうやってリスト化したいw。

BioProjectとかBioSampleを眺めてみた(続き)

普段の職場を飛び出してなので、(いまだに書類をこさえたり調整してメールしたりもありますが)コード書いたりするのがなんか久しぶりな気がします。ま、こうやってわいわいやった方が、ひとりでこもってやるよりシナジー効果が出ますわな。
昨日の続きとして、接頭辞のをちゃんと見たりもしてみた。(この後ろに数字が来る)

27780 SAMD
17 SAME
480372 SAMEA
3628607 SAMN

ついでにBioProjectでも同じことをしてみた。

288 PRJDA
3167 PRJDB
499 PRJEA
9964 PRJEB
143180 PRJNA

それから、仕事柄、文献とのリンクも見たいのだが、実際に文献で引用されている正例がとれずに、まだそこは野ざらしである。
本業の方だが、懸案の自動でAmazonのデータを更新の部分で、データベースをこちらでこさえて、転送するところまでは済。きちんとそれがcronでまわるのか(自分で打つ分には実行できるけど)、向こうできちんとロードできるかはまだなので、それはある程度 時を待ちたい。
あ、あとは15日をすぎたので、定例の月イチの統計情報の更新もdone
なかなか目に見えるものができないってのはアレですなぁ。

BioSampleを改めて眺めてみた

このブログも昨年は書いていなかったので、これは再開するのはいかがなものかと思ったのだけれども、よく見たら(自宅でブログを書いていたので)3年くらい書いていなかった時期があって、まぁいいか、と思い、また書いてみました。
次世代シーケンサーの公共データベースSRA(Sequence Read Archive)の検索エンジンを仕事で作ったりなどしています。もともとはSRAだけで閉じていたのですが、発現情報はもともとマイクロアレイのデータのレポジトリであるGEOにも入るようになり、つまり、発現情報はGEOとSRAの両方に入ることとなって、両方見ないといけなくなって、共通であるプロジェクト部分はBioProjectに、サンプル情報はBioSampleに外だしされることとなりました。
で、BioSampleのFTPサイトを改めて眺めてみたのですが、
DDBJの方





SAMD00000001


Bradyrhizobium sp. DOA9
MIGS Cultured Bacterial/Archaeal sample from Bradyrhizobium sp. DOA9

Bradyrhizobium sp. DOA9

NCBIの方は





SAMN00000002
19655
SRS000002


Alistipes putredinis DSM 17216

全然 違うじゃないですか。。。(NCBIの方はDDBJのを含んでいるようなのでこっちを見るか)
とりあえず、IDのprefix(ようするに頭の方)を調べてみました。とりあえず5文字分

27780 SAMD0
17 SAME5
480372 SAMEA
3628607 SAMN0

SAMと来て、DとかEとかNとか来て(DDBJ、EBI、NCBIなのは想像に難くない)で、なんちゃらだが、その後ろが数字かと言うとそうでもない、ということ。