Macintoshが発売されてから30年。自分とMacintoshは…と考えると、分子生物学の研究室に配属された1994年から、なのでちょうど20年。研究室にあって便利な文房具としてすぐに中古でColor Classicを買い、それで卒論も書き上げた。そう考えると付き合いはとても長い。 それから何台買ったのか、もう記憶に無い。現在の仕事道具(MacBook Pro)としても、自宅のオブジェ(Mac mini)としても、さまざまな役割をMacは担っている。例えば、今このブログを書いているのもMacだし、インターネット接続のためのWiFiルーターの役割を担っているのもMacだし、そのブログのサーバーとなっているのもMacである。そう考えると
No Macintosh, no life.
なのだなあとこのタイミングで再認識させられた。ありがとう、Macintosh、これからもよろしく!
一対一対応のデータなら対応付けは簡単。しかし、実際の研究活動においてはそんなことばかりでなく、特定のIDに対応づくデータが複数あり、それをまとめて表示することが必要な場合はままある。例えば、特定の遺伝子にアノテートされた機能情報などがそれである。1つの遺伝子について機能情報が綺麗にまとまっていれば対応付けも楽なのだが、そんなことはない。
そこで以下のような捨てコードで、同一のIDにアノテートされた機能情報をまとめて人間が見やすい**かたちに整形する必要があったりする。入力とするデータはID順にsortされている必要はないが、重複があるとその分、出てきた結果が見づらくなるので [shell] sort -u go_annotation.txt > go_annotation-unique.txt [/shell] などとして重複を除いておくべき。
1007_s_at ATP binding 1007_s_at collagen binding 1007_s_at metal ion binding 1007_s_at protein binding
のような入力データを処理して
1007_s_at ATP binding / collagen binding / metal ion binding / protein binding
のように同じIDのものは一行にまとめてくれるようにするには以下のようなPerlスクリプトで実現可能である …
タブ区切りテキストのカラム切り出しはUNIXコマンドラインでは、これまでも紹介したように [shell] cut -f2 arraytable.txt > column2.txt [/shell] みたいな感じで、arraytable.txtファイルの左から2番目のカラムだけを切り出して、column2.txtファイルに書き込むみたいなことを日常茶飯事的に実行する。 同じことをawkでもワンライナーで実行できて、 [shell] awk '{ print $2 }' arraytable.txt > column2.txt [/shell] とやればいい。
しかしながら、カラムの値にスペースが入っていたりするとこれがうまくいかなくなる。それはカラムをわけるデリミタ(とかファイルセパレータ)がdefaultではスペースもしくはタブという具合になっているからだ。そうなるとちょっと工夫が必要となってくる。 もちろんcutコマンドの-dオプションで指定してもいいのだが、敢えてawkでそれをやるスクリプトを紹介する。 [shell] BEGIN { FS="t" } { if( $16 > 1 ) { print $11 "t" $16 }} [/shell] この例では11カラム目がIDで、16カラム目が発現変動の比で、それが1以上になったものだけをIDとともにタブ区切りテキストで出力するというものである。この内容をcut …
サイト投稿によってコンテンツの増えた「ぼうのブログ」の見てくれ変更。ブログの右カラムにInstagramで直近にアップした画像が出るようにして、冗長になっていた「アーカイブ」や「直近の投稿」を「お休み」状態にしたり。Alpine PhotoTile for Instagram プラグインの設定がちょっと手間がいったが、いろいろとカスタマイズできそうで、これからもいぢって行くのではないかと。やっぱりブログ、楽しい。
今年も職場の同僚たちに明日に控えた誕生日を祝ってもらえた。それぞれ忙しい中それを忘れずに覚えていてくれたり、寒い中ケーキを買いに行ってくれたり…。大変嬉しいかったです。ありがとう。 前夜祭は独り、お気に入りの GERARD BERTRAND DOMAINES DE VILLEMAJOU 2011。 [amazon template=image&chan=default&asin=B00FE8BWYG]
まだまだ、我々のやっていること、広まっていない。個別に、じみじみと広めていかねば。 友人の学問に対する前向きな姿勢、本当敬服。息切れしかかっていた私にとって、なによりの「チャージ」。いろいろあるけど、「夢の実現」から「継続」へとコマを進めていかねば。
年末の出張の事務書類を出すところからスタート。その後、大ボスに呼び出されるなど。某仕様書、急ぎます。遅れていてすみません…。 昨晩、突貫工事で「統合」したブログの微調整も。日付がキッチリと入らないエントリがあるのはしょうがないのか?
昨年2013年は、以前より進めてきた
次世代シーケンサからのデータアーカイブSequence Read Archive(SRA)を使い倒すための仕組みDBCLS SRAの論文
の2つを世に送ることができた。2013年年頭に掲げた目標
「実際に、現場で、使った論文等を世に送り出す」
はなんとか達成できたのではないかと。それに慢心せず、本年2014年は、これまではやってこなかった新しいことに挑戦していきたい。もちろん、これまでどおり、DBを使って研究をされている皆さんとの出会いを大切に、「知のめぐり」の向上にDBを使いやすくする仕組みやコンテンツの作成維持管理に努めていくこともやった上で。
来年度頭には、DBCLSの柏の葉キャンパス駅前と国立遺伝学研究所(三島)内への移転が予定されており、変化の年になるはず。遺伝研の中に仕事場を頂くからには、現場で自然に使われる「統合DB」を目指して、他の遺伝研の部門と融合的な活動を推進していきたい。
何エントリか前に紹介した遺伝子対応表を作る際に、一般的には双方向ベストヒットということをよくやります。これは2種類の生物種にコードされた遺伝子セット全体を比較する際、クエリとした生物種1とデータベースとした生物種2でBLASTをかけた後にクエリに生物種2、データベースに生物種1を割り当てて逆に計算し、その対応関係において双方向で相思相愛になるペアだけを選ぶというものです。詳しくは坊農の博士論文[PDF] 2.3オーソロガス遺伝子の推定(p.12-14)を見てください。それをやる具体的な手順を今回紹介します。
BLASTをかけて-outfmt 6で得られた結果に対して、cutコマンドを使って1カラム目と2カラム目だけを週出してuniqコマンドで重複を除きます。 [shell] cut -f1,2 kaiko_human.txt |uniq > kaiko-human1.txt [/shell] また逆のBLAST検索の結果も同様に処理します。 [shell] cut -f2,1 human_kaiko.txt |uniq > kaiko-human2.txt [/shell] で、得られた結果の2つのファイルの同じ部分が双方向ベストヒットのペアということになります。 それを抽出するにはちょっとトリッキーですが、以下のようにしています。 [shell] cat kaiko-human1.txt kaiko-human2.txt | sort …
前回の配列セットを取得するはTogoWSを使うものでした。TogoWSは汎用で前回例で示した塩基配列データベース以外のデータベースにも使えるので便利なのですが、取得が1エントリづつで遅いという欠点がありました。
塩基配列に限っては、DDBJのgetentryを使うという方法が数万のエントリ取得において劇的に早いことが関係者の方の手助けを得て分かりました。 前回紹介したAK377185-AK388575の連番エントリを一気に取るには、以下のコマンド(curl)でいいようです。
[shell] curl -o fasta.txt.gz 'http://getentry.ddbj.nig.ac.jp/getentry/na/AK377185-AK388575/?format=fasta&filetype=gz&trace=false&show_suppressed=false&limit=0' [/shell]
これでfasta.txt.gzに一万超エントリのFASTA形式ファイル(+gzip圧縮)が得られます。もちろん、シングルクオートの中身(URL) [shell] http://getentry.ddbj.nig.ac.jp …
cDNA配列のセット(多くの場合数万のオーダーの配列群)が論文発表されたのに、それがウェブサイト上にないとかいう話をよく聞きます。その場合、大変申し訳無いのですが論文の読み込みが足りないと言わざるを得ません。論文発表されたからには必ずデータは公共データベースに公開されているはずです。していないのは「公開した」とは言えないからです。 ＃まれにNGS黎明期のゴタゴタで配列を公開していないのに論文が通ってしまっている例がありますが、それは例外として。
論文にはデータベースのアクセッション番号がどこかに書いてあるはずです。まずそれを探しましょう。例えば、
All full-length cDNA sequences obtained in this work appear in public databases under accession nos. [DDBJ:AK377185-AK388575].
とあると、アクセッション番号AK377185-AK388575で公開しているわけですから、以下の様なPerlのコード(get.prl)で配列が取得できます。 [perl] my $togows = "http://togows.dbcls.jp/entry/nucest"; #baseURL my $start …
本日(2013年7月25日)、PLOS ONEから論文を発表しました。お世辞にもいい機能アノテーションがなされているとは言えない「非モデル」生物種カイコにおいてマイクロアレイデータ解析を可能とし、ヒトで使われているパスウェイ解析ソフトが利用可能となったことが解析の「突破口」となりました。それを可能としたのはカイコとヒトの「遺伝子対応表」です。
もちろん、あらかじめ個別のモデル生物データベースセンター的なサイトで用意されてあればそれを使うのがいいのですが、必ずしもそれがあれば利用者側で何もやらなくていいというわけではありません。今回の場合もそうだったのですが、そういった提供されているデータセットと、マイクロアレイ解析の場合など調べたいデータセットとのIDの対応関係が簡単にとれない状況であることもあります。そういった場合、やりたい解析(今回の場合、カイコをヒトと比べてこの二種類の生物種間で対応表が作成したかった)に合わせて「対応表」を個別に作成する必要があります。それをどうやったかを手順を追って解説したいと思います。
まず、どこのデータセットを比較対象とするかの取捨選択から解析は始まります。世の中にはすでに利用可能な、たくさんのデータセットがあります。そのうち、どれを使うべきか?これには実は答えがなくて、それぞれの状況や好みに応じて選ぶしかないと思います。UCSCのサイトのものを使うやり方もありますが、Ensemblの方が整然としていて使いやすいので、今回もEnsemblの方を使いました。ここにあるデータダウンロードindexページから必要なデータをダウンロードしてきます。今回の場合、カイコアレイの元のEST(cDNA)配列に対して、ヒトのcDNA配列の対応を付けたかったので、それに相当するHuman cDNA(FASTA)をクリックします …
マイクロアレイデータを再利用する際、遺伝子方向のみならずサンプル方向のデータの相関が気になってくると思います。justRMAでnormalizeしたデータで次に行うのは、おそらくサンプル方向のPCA(Principal Component Analysis:主成分分析)でしょう。ここではRでPCAをして第1,2主成分でグラフを一気に作成するコマンドを紹介します。
まず以下の様な内容でPCA.rという名前のファイルを作成しておきます。もちろん、Rを起動して一行ずつコピペしていっても構わないのですが、最後の行にあるパラメータをいじったりするなど何回もこれを実行することになると思うのでそういう場合はファイルに書き下して「スクリプト化」しておいたほうが楽です。 [shell] data <- read.table("RMA.txt", header=TRUE, row.names=1, sep="t", quote="") data.pca <- prcomp(t(data)) names(data.pca) plot(data.pca$sdev, type="h", main="PCA s.d …
ファイルが何行あるか知りたいだけなら、いちいちテキストファイルをWordで開いて文字カウントツールを使うまでもありません。そのファイルをfile.txtとするとコマンドラインでは [shell] wc -l file.txt [/shell] で良くて、文字数とかバイト数も合わせて知りたい時には [shell] wc file.txt [/shell] で知ることができます。
また、出現する文字の種類が何種類あるか知りたいだけならsortとuniqコマンドを使って [shell] sort file.txt | uniq | wc -l [/shell] すれば良い(ただし、この場合は一行に一単語というファイルの中身の前提条件があって、だが)。バイオインフォマティクスな業界ではそれが何回出てきたかが気になる状況が多々あります。そういう需要に対応するために以下の様な簡単なPerlスクリプト(count.prl)を2010年代の今になっても頻繁に利用しています。 [perl] while(<>) { my($word) = split; $num{$word}++; } foreach (sort keys %num) { print "$t …
最近のバイオなデータ解析では、ファイルサイズが大きすぎたり、または多すぎたりで、GUIな環境では快適にファイル操作できないこともあって、コマンドラインでデータ解析をする人が増えてきたのはいい流れと思っています。基本的なシェルの操作は私がまた一から書くようなものではないと思うので、そのレベルの情報が必要な場合は、若葉本(初心者でもわかるバイオインフォマティクス入門)の「第2章必要なコンピュータ環境を作ろう」p38-92とか、Mount本(バイオインフォマティクスーゲノム配列から機能解析へ)第2版の「付録: UNIXをはじめよう」p565-570とか参照してください。
割と当たり前すぎて紹介されにくい、普段から私などがよく使う(技というにはおこがましいという人がいるかもしれないレベルの)tipsを何点か紹介します。
コマンド入力は、誰しもほぼ間違いなくタイプミスします。それを防ぐためにどうしているかというと、タブ補完を利用します。コマンドを途中まで入力して(たとえばcuff)、そこでタブ(tab)キーを入力するとコマンドサーチパスにその単語で始まるコマンドがsuggestされます。 [shell] % cuff cuffcompare cuffdiff cufflinks cuffmerge [/shell] さらにlと打って、タブキーを押すとコマンドすべてを入力しなくてもコマンド名が「補完」されてcufflinksとコマンドライン上に表示されます。 cufflinkだったか、cufflinksだったか、うろ覚えで、それに続くさまざまなパラメータを書いてしまってcommand not foundと言われてがっかりしなくてすみます。
直前に使ったコマンドを一部オプションを変えて再実行したいということが頻繁にあります。その際 …
マイクロアレイデータを再利用する場合、それらのデータのnormalizeは必須です。 GeneChipの場合、normalizeしたいCEL形式のファイル群を同じディレクトリに入れておいて、一気にやってしまうことがjustRMAで可能です。
まずRが入っているかどうか。入っていない場合、Mountain lion(10.8.4)の場合、homebrewで [shell] brew install -v R [/shell] とするだけでRの最新版(2013/6/17現在、3.0.1)をインストールできます。 WindowsだとCRAN、たとえば筑波大のCRAN mirrorからダウンロードするとよいでしょう。
[shell] source("http://bioconductor.org/biocLite.R") biocLite("affy") library(affy) write.exprs(justRMA(), file="RMA.txt …
2010年代になった今でも相変わらず使用しているプラットフォーム(主にOS)の違いによって受渡したデータが閲覧できない等の問題は多く発生しているようです。
一番顕著なのはAJACS名古屋のテキストで詳しく書いたUNIX, MacOS, Windowsのテキストファイルでの改行コードの違いの問題でしょう。これらは歴史的な経緯があってこうなってしまっていて、普段はアプリケーションレベルで回避しているので問題にならないものの、自らが生データを扱うようになるととたんに問題の原因となっているようです。
OSの選択肢の傾向(データを渡す側はWindowsで受け取る側はUNIX)から一番良く起こると考えられる、WindowsのテキストファイルをUNIXのそれに変換する場合、以下のperlのワンライナーで変換できます。
perl -pe 's/r//g' win.txt > unix.txt
変換したいファイル名を変えずに上述の改行コード変換だけをしたいのなら以下のワンライナーがおすすめです。
perl -i~ -pe 's/r//g' data.txt
data.txtの改行コードが変更され、元ファイルはdata.txt~というファイル名に変更されます。
また、改行コード問題以前に、半角の￥(円記号)と\(バックスラッシュ)が違いない(同じ意味である)ということについて …
遺伝子上流配列など、ゲノム配列の任意の場所を切り出すのには、BLATパッケージに含まれている nibFrag をこれまで勧めてきました。FASTAフォーマットの配列ファイルに対して別個にnibインデックスを作らないといけなくて管理が面倒なのと、やはりBLATはfor profitには有償という大きな難点がありました(ライセンス料は決してお安くない)。
そこでpublic domainなソフトウェアのNCBI BLASTパッケージにも似たプログラムがきっとあるはずだろうということで探してみると… fastacmd コマンドがそれだったようです。しかしながらBLASTパッケージのバージョンアップでBLAST+になってからは fastacmd が含まれなくなっているようです。更に調べると blastdbcmd というプログラムがその後継という位置づけのようで、BLAST+をインストールする時に blastdbcmd は同時にインストールされるプログラムです。その使い方を調べてみました。
まず、makeblastdb コマンドでBLAST用のindexを作成します(かつての formatdb コマンド相当だと思います)。-dbtype nucl で核酸配列であることを指定し、-hash_index でインデックスの作成も指定しておきます(追記: -parse_seqids を足さないと個別エントリの切り出しには対応しない模様。別のblastdbcmd関係エントリ参照。ここで -in で指定している 1.fa はゲノム配列(1番染色体)が1本だけ入ったFASTAフォーマットのファイルです。
1 … |
MountainLion(MacOSX10.8.4)のhomebrewでBLASTがとくにいぢらなくても入るようになった模様です。バージョンの確認は、左上のリンゴマークメニューから「この Mac について」を選んで出てくるウインドウで確認できます。 homebrewのインストールは「homebrew インストール」でググって入れてもらうとして、ここではlocal BLAST(以下、単にBLASTと書きます)を入れる際につまづきそうなことだけ。 普通に [shell] brew install blast [/shell] とやってもBLASTなんて知らんといわれてインストールされません。 [shell] brew tap homebrew/science [/shell] としてレポジトリを追加する必要があります(この場合、homebrew-science)。そうしてから [shell] brew install -v blast [/shell] とすると…2013年6月6日現在だと2.2.28のBLASTのコンパイルが始まります。コンパイルには結構時間がかかると思いますが、終わったら。 [shell …
中5週間。今週外部での講演もあるので。
GUIとはGraphical User Interfaceの略で、今使っているInternet ExplorerやGoogle Chromeなどのウェブブラウザのように主にマウス操作と最小限のキーボード入力でコンピュータが動いてくれるインターフェースのことです。ウェブブラウザ以外だと、表計算以外の目的でも頻繁に使われているExcelとか遺伝子発現解析ソフトウェアの定番GeneSpringなどがGUIなソフトウェアということになります。それらの操作法については、統合TVで動画としてプロトコールの行間を埋める形で多数紹介してありますので、そちらを御覧ください。 2010年代の最近ではGUIでも高度な操作ができるようなものも使われるようになってきました。例えば、Biomartを使うと遺伝子の上流配列を取得することが簡単にGUIでできますし、Galaxyを使うと特定の転写因子予測結合領域と遺伝子上流領域の「交差点」をリストアップすることが可能です。しかしながら、一般的にはGUIでは研究者が論文の図や表として求めるレベルのデータを作るには機能不足という状況がよくあります。 それに対してCUIとはCharacter User Interfaceのことで、コマンドをキーボードから直接入力してコンピュータにジョブを投入する方式のことです。CUIなツールに関しても一部、LocalBLAST、cygwinやBioperlの使い方を通して統合TVでも紹介してきましたが、そのバリエーションの多さ故にそのありがたさの多くを伝えきれていません。 そこで、ここではその「バリエーション」に敢えて向き合い、いろんな使い方を「プロトコール」として紹介していきたいと思っています。
電子化されコンピュータ上で扱えるデータが主に核酸配列とそれを翻訳して得られるアミノ酸配列だったこともあり、それらのデータ解析（配列解析）がバイオインフォマティクスと呼ばれてきました。その流れで書かれた教科書が我々が翻訳したMountのBioinformaticsです。 マイクロアレイの発明により遺伝子発現が網羅的に測定できるようになるとそこから得られた数値のデータもバイオインフォマティクスで扱うデータの対象となってきました。それを考慮して第二版のMountのBioinformaticsはマイクロアレイデータ解析の部分が加筆されています。 [amazon asin=4895924262&template=thumbnail] また、知識として蓄えられた代謝経路などのパスウェイのデータもバイオインフォマティクスとして認知されるようになってきました。現在ではそれらのデータは、人類共通の知的な財産としてWikipediaなどにまとめられ(Gene WikiやWikiPathwaysなど)、誰でも自由に使えるようにする流れになっております。 それらのデータを使ったシミュレーションもバイオインフォマティクス、さらには測定機器から得られるイメージデータのデータ処理もバイオインフォマティクスと呼ばれているのが現状だと思います。さらには、今後得られるであろう臨床データもバイオインフォマティクスとされていくのではないかと思います。つまり、ライフサイエンスを研究すること＝バイオインフォマティクスになりつつあるのです。
[caption id="attachment_35" align="alignnone" width="300"]
バイオインフォマティクスの範疇[/caption]
中7週間。大分暖かくなってきたので。
私は、約十年前に「初心者でもわかる!バイオインフォマティクス入門―やさしいUNIX操作から遺伝子・タンパク質解析まで」と題した本を書きました。表紙などにある若葉マークが印象的なため、「若葉本」とその筋の人たちには呼ばれています。2001年の秋から1年間、羊土社刊の「Bioベンチャー」誌に創刊号から連載させていただいた、ぼうのうによる「バイオインフォマティクス入門」の連載記事を纏めて加筆した本です。とくに第2章「必要なコンピュータ環境を作ろう」、第4章「最新の研究をのぞいてみよう」は完全に書き下ろしです。情報科学方面から参入しようとしている方には生物学者に対する間違った認識による行き違いのないように、とくに第1章4節「バイオインフォマティクスと情報科学」の一読をお勧めします。 ぼうのうがよく周りから訊かれること、普段から考えていることを文章にした本で、また本にすることによって考えがまとまり非常にすっきりしました。本当にいい経験をさせてもらい、また周りの人にも自信を持って勧めれる本になったのでいうことなしでした。[amazon asin=489706290X&template=thumbnail]
それと同時に海外でも使われているBioinformaticsの教科書の翻訳、監訳にもたずさわり、バイオインフォマティクスが日本で普及することに加担してきたつもりです。とくに、Cold Spring Harbor Pressから出ているBioinformaticsは、Mountさんが著された本なので、「Mount本」と呼ばれています。日本語訳は 「バイオインフォマティクス ゲノム配列から機能解析へ …
twitterを使うようになって、ブログをあまり書かなくなった。しかし、やはり、まとめた考えを発信したいと思い、「某のウェブログ」改め「ぼうのブログ」として再開します。
「2010年代のバイオインフォマティクス」というカテゴリで新しい企てを始めています。詳しくは、「はじめに」を見てください。
中4週間。今回はうちで。
なる会議を100kmほど離れた別部署にて。相互理解が深まったのではないかと勝手に確信。
学会(の年会)じゃないタイミングで参加するのはおそらく初めて。主賓の一人も神戸でのオフ会で知り合った人だったが、もう一人も実は知られていたというオチ。世間は狭い。
復帰2回目。有酸素運動にウォーキング25分、筋トレ、水泳25m往復12回。利用券を使い果たしたので、早速次を頼んでおく。
昨年と打って変わって、不気味なほどにdisられず、むしろ恐ろしい。総評で言われた「対等なコラボをして、biologistだけじゃできない成果を出す」ことの実現に向けて、さらに頑張って行きたい。