MotDB


AJACS32/bono

講習会のページに戻る

コマンドラインで遺伝子配列を解析する


目次


事前調査: Mac: 4, Windows: 16, N/A: 1

_ はじめに: 参加者アンケート

  • おもにMacの人:5
  • おもにWindowsの人:9
  • iPhoneの人:4
    • そのうちiPhoneは体の一部の人:2
  • iPadの人:4
  • バイオインフォマティクスを使っている人:0
  • プログラミングをたしなむ人:0
  • データベースを構築している人:0
  • データベースを構築したい人:0
  • ライフサイエンス統合データベースセンターの存在を知っていた人:3

_ 必要なコンピュータリテラシー

_ Google

http://www.google.co.jp/

Google(「グーグル」と読みます)検索することを「ググる」といいます。そこでインターネット上では、自分でインターネット検索もせずにあれこれと質問をしてくるユーザーに対して以下のようにいうことがあります。

ググれカス

  • 【実習B1】DBCLSでググりなさい。何件ヒットがありますか?

    ←こたえは左の+マークをクリックすると出てきます

    2012年7月26日現在、約 1,040,000 件

でもこれは「ウェブ全体から検索」した結果なのです。

  • 【実習B2】ググった結果のページ中で「日本語のページを検索」のリンクを探してクリックしなさい。

    ←こたえ。上記の結果と比較してどう変化するでしょうか?

    2012年7月26日現在、約 593,000 件

  • 【実習B3】さらに画面最下部の「検索オプション」をクリックして絞り込みをかけてみなさい。ドメイン.ac.jpを指定すると何件ぐらいに絞り込まれますか?

    ←こたえ。得られる結果にはどういった特徴があるだろうか?

    2012年7月26日現在、約15,700 件で、ドメインがdbcls.rois.ac.jpのサイトばかりがヒットしてくるといった特徴がある

  • 【応用B4】さらに見るべきヒットを絞り込むにはどういうオプションを指定すればいいだろうか?

    ←こたえ。

    例えば、「キーワードを含めない」オプションでroisを指定してみる

  • 【応用B5】'DBCLS' は「ライフサイエンス統合データベースセンター」の略号であるが、たまに「ライフサイエンス総合データベースセンター」と間違えられる。そう間違えられている例を"で囲うことでインターネット検索エンジンを用いて抽出しなさい。

    ←こたえ。

    "ライフサイエンス総合データベースセンター”でググって、「キーワードを含めない」オプションでライフサイエンス統合データーベースセンターを指定してみる

  • 【実習B6】'SPF'でググリなさい。どういったことが起こるか?

    ←こたえ。

    例えば、'DBCLS'でGoogle検索しても「ライフサイエンス統合データベースセンター」以外の'DBCLS'はインターネット上にそれほどないため困らないのであるが、短い略語の場合は同義語がインターネット上に多数存在して調べたい情報に行き着くまでに非常に苦労することになります。

つまり、こういうことです。

ググるなあぶない

_ 脱線: GGRNA

http://ggrna.dbcls.jp/

遺伝子名や断片配列(塩基配列やアミノ酸配列)でRefSeq?を検索できるサイト。DBCLS謹製。

http://lifesciencedb.jp/image/small_video_icon.png GGRNAで遺伝子をGoogleのように検索する

_ cygwin(windowsの人のみ)

cygwinをまだインストールしていない人は以下の動画を参考にインストールしてください。

http://lifesciencedb.jp/image/small_video_icon.png WindowsでUNIX! 1. Cygwin インストール編

_ UNIX コマンド

http://lifesciencedb.jp/image/small_video_icon.png WindowsでUNIX! 2. ファイル操作編

  • cd
  • pwd
  • ls
  • cp
  • mv
  • rm

http://lifesciencedb.jp/image/small_video_icon.png WindowsでUNIX! 3. ファイル操作応用編

  • wc
  • less
  • gzip
  • tar
  • mkdir
  • find

ここで紹介した以外に大きなファイルを扱う際に先頭の数行だけを表示する head、末尾の数行だけを表示する tailなどがあります。

% head sample1.fa

% tail sample1.fa

細かいオプションなどはググるまでもなく

% man head

などすることでマニュアルが参照できます。

_ Local BLAST

できる人は統合TVみてどんどん進めてってください。

ひょっとすると今使っているLANの制限で、以下動画の中で指示のあるFTPサイトにアクセス出来ないかもしれません。 BLASTのダウンロードサイト(臨時)として

http://dbcls.rois.ac.jp/~bono/docs/120728AJACS/LocalBLAST/

DBのダウンロードサイト(臨時)として

http://dbcls.rois.ac.jp/~bono/docs/120728AJACS/db/

を用意したので、こちらから必要なファイルをダウンロードしてください。

Windowsのcygwin上でのLocal BLASTを試みる人は、BLASTのプログラムのFTPサイト(のミラーサイト)

http://dbcls.rois.ac.jp/~bono/docs/120728AJACS/LocalBLAST/ から

ncbi-blast-2.2.26+-ia32-win32.tar.gz (32ビットのwindowsの方)

ncbi-blast-2.2.26+-x64-win64.tar.gz (64ビットのwindowsの方)

をダウンロードし、以下の統合TVのMacOSX版の方でやってみてください(中上級者向け)。その際、展開したプログラムファイルは

cp ~/ncbi-blast-2.2.26+/bin/* /usr/local/bin

などとしてコマンドサーチパスが通ったところにコピーしておいて下さい。

_ Windows

http://lifesciencedb.jp/image/small_video_icon.png Local BLAST の使い方〜導入・準備編〜 2011

http://lifesciencedb.jp/image/small_video_icon.png Local BLAST の使い方〜検索実行・オプション編〜 2011

使用するテスト配列のファイルはこちらです。

_ MacOSX

http://lifesciencedb.jp/image/small_video_icon.png Local BLAST の使い方〜導入・準備編(MacOSX版)〜 2011

http://lifesciencedb.jp/image/small_video_icon.png Local BLAST の使い方〜検索実行・オプション編(MacOSX版)〜 2011

使用するテスト配列のファイルはこちらです。

_ データ後処理

_ 数える

大規模にBLASTを計算するといくつヒットがあったかなどが一目ではわからなくなってきます。そういった「数える」必要があるときにコマンドラインでの処理が役に立ちます。

  • FASTAファイルの配列のエントリ数を。例えば、上で用いたyeastのアミノ酸配列のDB(FASTA形式)の配列エントリ数を数える場合。

    grep -c ^\> yeast.aa

6298と出れば正解です。これが(データベース中の)yeastのアミノ酸配列数です。

  • BLASTの結果ファイル(-outfmt 6オプションでタブ区切り出力したもの)で、ヒットのあったエントリ数を。自分で結果ファイルを用意しても構いませんが、ここでは以下の結果ファイルのサンプル(ファイル名:testfmt6.out)で説明します。

処理する前にこれがどんなファイルか見ましょう。

less testfmt6.out

一行が長くて見づらいですね。lessのオプションに-S (Sは大文字のS)というものがあります。すると…。

less -S testfmt6.out

見やすくなりましたね。

  • 1カラム目(一番左のカラム)を抜き出し(cut -f1)、ソートして重なりのない状態にして(sort -u)出力する。

cut -f1 testfmt6.out | sort -u | less

さらに、数える(wc)。

cut -f1 testfmt6.out | sort -u | wc

これによってヒットがあったクエリの数が数えられます。これぐらいの数なら人の目でも数えられますが、多くなると不可能になってくるので、こういうコマンドが大変重宝します。

  • ヒットのあったDB中のエントリが記載されている2カラム目を抜き出し(cut -f2)、
    • ソートして(sort)数える(wc)。

      cut -f2 testfmt6.out | sort | wc

  • ソートして重なりのない状態にして(sort -u)数える(wc)。

    cut -f2 testfmt6.out | sort -u | wc

この違いは重複したエントリによるものです。それが何回でてきたか知るには?

  • ヒットのあった数をそれぞれのエントリごとに。count.prlというPerlスクリプトを使って。

which perl

といれて

$ which perl

which: no perl in (/usr/local/bin:/usr/bin:/cygdrive/c/ProgramFiles?/Parallels/ParallelsTools?/Applications:/cygdrive/c/Windows/system32:/cygdrive/c/Windows:/cygdrive/c/Windows/System32/Wbem:/cygdrive/c/Windows/System32/WindowsPowerShell?/v1.0:/cygdrive/c/Program Files/MacType?)

のように no perl と出たらPerlがインストールされていません。cygwinのインストール画面でオプション設定でPerlをインストールしてください。ちなみにMacOSXの場合は買った状態でPerlがインストールされています。

while(<>) {

my($word) = split;

$num{$word}++;

}

foreach (sort keys %num) {

print "$_\t$num{$_}\n";

}

cut -f2 testfmt6.out | sort | perl count.prl | less

_ ベストヒットだけを抽出

BLASTの出力ファイル(-outfmt 6オプションでタブ区切り出力したもの)をテキスト処理するやり方もありますが、Local BLASTが実行可能となった今、普通にBLASTのオプションで -max_target_seqs 1 を指定するほうが早いです。例えば、

blastp -db yeast.aa -query mito.aa -out mito_vs_yeast.out -outfmt 6 -max_target_seqs 1

これでquery: mito.aa db: yeast.aaのベストヒットが抽出できます。

_ おまけ:改行コード問題

出てきたBLASTの結果。テキストファイルなのにある種のソフトウェアではうまく認識されないことがあります。その場合、この問題を疑ってください。

_ 改行コード問題とは?

使用プラットフォーム(OS)によって改行コードが異なるため、、

  • Macのテキストファイル

    % od -c mac.txt

    0000000 c a r r i a g e r e t u r n \r

    0000020 m a c

  • Unixのテキストファイル

    % od -c unix.txt

    0000000 l i n e f e e d \n u n i x

    0000016

  • Windowsのテキストファイル

    % od -c win.txt

    0000000 C R a n d L F \r \n w i n

    0000016

_ その対処法

いくつかやり方がありますが、一番汎用性の高いものを紹介します。

以下のperlワンライナーでmac改行形式のファイル(mac.txt)がUNIX改行形式のファイル(mac_conv_unix.txt)に変換されます。

perl -pe 's/\r/\n/g' mac.txt > mac_conv_unix.txt

つぶやけ、されば救われん

 
添付ファイル: filetestfmt6.out 623件 [詳細] filetestdouble.txt 588件 [詳細] filetestnt.txt 618件 [詳細] filetest1.txt 631件 [詳細] filetest.txt 617件 [詳細] filesample.zip 465件 [詳細] filewin.txt 581件 [詳細] fileunix.txt 629件 [詳細] filemac.txt 609件 [詳細]
 
Link: AJACS32(2488d)
Last-modified: 2012-07-28 (土) 11:38:33 (2488d)