2017.08.04版

スライド30までは自習。8/28の自習日
か、8/29の10:30までにスライド31から始められる状態にしておいてください

NGS解析（初～中級）
ゲノムアセンブリ後の各種解析

東京大学・大学院農学生命科学研究科
アグリバイオインフォマティクス教育研究プログラム
門田幸二（かどた こうじ）
kadota@iu.a.u-tokyo.ac.jp
http://www.iu.a.u-tokyo.ac.jp/~kadota/
Aug 29-30 20171

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017

2

予習事項


平成28年度NGSハンズオン講習会①乳酸菌ゲノム配列決定論文のPacBioデータを用いて、de novoアセンブリ後の各種解析を行う。②主目的は、LinuxおよびNGS解析周辺のスキルアップ。昨年度からの継続性も重視

 日本乳酸菌学会誌NGS連載第7回のウェブ資料W9-3までの内容を実施



平成29年度NGSハンズオン講習会
 8/29-30は、第7回の残りと第8回の内容が中心 ②①
Aug 29-30 2017

Tanizawa et al., BMC Genomics, 16: 240, 2015

3

予習事項

乳酸菌NGS連載第4回までは自習。第5-7回途中まで
は、赤枠内の平成28年度講習会の8/1-3で実施済み。特に、①2016年8月3日の講習会資料のスライド121
以降を眺めておき、どんなデータのアセンブリ後の解析をやろうとしているかは把握しておきましょう①

Aug 29-30 2017https://biosciencedbc.jp/human/human-resources/workshop/h28-2

4

（実質的に同じなので）もちろん①乳酸菌
NGS連載第7回の、②ウェブ資料のW9あたりまでを眺めておくのでも構いません

予習事項

①

②

Aug 29-30 2017

5

第7回W9-2で、DDBJ Pipeline上で実行した
HGAPアセンブリ結果のページを眺めている。①
Total contig sizeは、乳酸菌の一般的なゲノムサイズ（2-3Mb）と近く妥当。②Maximum contig
sizeのものが全体の9割以上を占めていることから、これが乳酸菌の染色体ゲノムなのだろうと妄想する。③コンティグ数は4つ。このデータの正解は3つ（1 chromosome and 2 plasmids）

W9-2：結果を眺める③①②

Aug 29-30 2017

6

この後の展開は…


2016年8月3日のNGS講習会資料の最終ス
ライド。平成29年度講習会では、この中の一部を講習会用にアレンジして解説します

W10：multi-FASTAファイルの分割
 プログラムによっては、single-FASTAのほうが取扱いやすい場合もある



W11：FASTQファイルの分割とクオリティスコア分布
 FASTAファイル用がうまく動かなくても、4行で1リードだということが分かっていれば、Linuxコ

マンドでどうにかなる、という話
 PacBioはFASTQファイルも出力する。Rで(single-)FASTQファイルを読み込んでクオリティス
コア分布を描画し、PacBioデータは両末端のクオリティが低い傾向にあることを確認

環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認
 W18：NCBI blastのやり方も示し、様々な解析手段を伝授

Aug 29-30 20177

参考Tips：名前の変更①のBioLinux8という名前をNGS20170829に変更する。②電源オフになっている状態で、③設定③①②

http://www.iu.a.u-tokyo.ac.jp/~kadota/bioinfo_ngs_sokushu_2016/20160719_tips_20160704.pdf
Aug 29-30 20178

①の部分を編集して、BioLinux8という名前を参考

Tips：名前の変更①

Aug 29-30 20179

①NGS20170829に変更して、②OK

参考

Tips：名前の変更①②

Aug 29-30 201710

参考 ①無事NGS20170829になりました。この部分はただの識別用。8/31、および9/1のovaファイルの名前と被るときや自分好みの名前に変えたいときにご利用ください

Tips：名前の変更①①

Aug 29-30 201711

①でBio-Linuxを起動。②パスワードは
pass1409です（連載第3回W2あたり）

起動
①

②

Aug 29-30 2017

Bio-Linux: Field et al., Nat Biotechnol., 24: 801-803, 2006

12

（もしこのようなポップアップが出たとしたら）①Don’t Upgrade。やってもいいとは思いますが、やったことがないのでその後の保証は（やらなくてもですがw）できません

起動後①

Aug 29-30 2017

13

①OK

起動後

①

Aug 29-30 201714

①Terminal

ターミナルを起動①①

Aug 29-30 201715

こんな感じになります

ターミナル起動後

Aug 29-30 201716

backupがないのはダメ①①lsで確認。②backupという名前のディレクト
リがないヒトは平成29年度講習会用として指定されたovaファイルを正しくインポートできていないので、指示通りにやり直してください②

Aug 29-30 201717

参考URL①ここにコピペ用コマンド集がありますので、できるかぎりこれを有効活用してください。次のスライドで見えているresult.zipファイルは、②からダウンロード可能です①②
Aug 29-30 201718

W9-3：ダウンロード1

PacBio用de novoアセンブリプログラムHGAP実行
結果である、①result.zipというzip圧縮ファイルを共有フォルダ（ホストOS側はDesktop/share）にダ
ウンロード。②MD5については、連載第3回W12で説明した。③リンク先のMD5チェックサム値

①

②

③

Aug 29-30 2017

HGAP: Chin et al., Nat Methods, 10: 563-569, 201319

①cdコマンドで共有フォルダに移動

W9-3：ダウンロード2①

Aug 29-30 201720

①pwdで作業ディレクトリをフルパスで表示

W9-3：ダウンロード3①

Aug 29-30 201721

W9-3：ダウンロード4①wgetコマンド（初出はたしか第3回W11-1）でアセンブリ結果ファイル(result.zip)をダ
ウンロード。ここでは、②オプションを-cのみとして、途中経過を表示させている②①

Aug 29-30 201722

①リターンキーを押した直後

W9-3：ダウンロード5①

Aug 29-30 201723

W9-3：ダウンロード6①ときどきネットワークの不調で、こんな感じで失敗することがあります。そんなときは…

①

Aug 29-30 201724

W9-3：ダウンロード7

数分後にwgetをリトライするのもアリですが、①のような感じで、②~/backupディレクトリ内にあるresult.zipを、③カレントディレクトリにコピーして、wgetしたつもりになってください①②

Aug 29-30 2017③25

①ls –lで確認。②約13MBの
result.zipが確かに存在します

W9-3：ダウンロード8①②

Aug 29-30 201726

W9-3：ダウンロード9①md5sumでresult.zipの②MD5チェック
サム値（第3回W12）を表示。この値と…

①②

Aug 29-30 201727

W9-3：ダウンロード9（今回はダウンロード元がDDBJ Pipelineではなかったが、そうだったとして）①をクリッ
クして見られる②MD5チェックサム値と比較し、同じであればファイルの中身は同じである（ファイルが壊れていない）と判断する①②

Aug 29-30 2017

DDBJ Pipeline: Chin et al., Nat Methods, 10: 563-569, 2013

28

ゲストOS（Bio-Linux）上の①のディレクトリは、ホストOSのデスクトップ上にあるshareフォルダとつながっています。そのような設定を行ったからです

共有フォルダ①

Aug 29-30 201729

共有フォルダ

ゲストOS（Bio-Linux）上の①のディレクトリは、②ホストOS（この場合Windows 10）のデスクトップ上にあるshareフォ
ルダとつながっています。つまり共有フォルダです。③確
かに同じresult.zipファイルが見えていますね。この（程度の）話についてこれるヒトが予習をやったと言えるヒト

②
②
③
①

③

Aug 29-30 201730

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 201731

W9-4：解凍して概観①result.zipを解凍。zip圧縮フ
ァイルの解凍コマンドはunzip①

Aug 29-30 201732

W9-4：解凍して概観①lsで中身を確認。②resultディレクトリの中には、③赤枠で示す計4つのファイルが存在する①②③

Aug 29-30 201733

W9-4：解凍して概観①がPacBio用のde novoアセンブリプロ
グラムHGAPの主な実行結果。拡張子が.fastaなので、multi-FASTA形式ファイル①

Aug 29-30 201734

Tips：ls -ld①ls –ldとすれば、ディレクトリの中身を非表示にできる(第6回W11-2)。②ls –lとの違いは③一目瞭然

②

③

①
③

Aug 29-30 201735

第7回原稿p104右下

今は①赤枠あたりです。このあと②赤
下線部分あたりの話をしていきます①②

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 201636

W9-5：コンティグ数①resultディレクトリに移動し、②コンティグ数を表示。FASTA形式ファイルなので”>”を含む行数がコンテ
ィグ数に相当する。③description行を表示。こんな感じの記述内容か～と思うだけ。④行数は、40,567行

①

②
③

④

Aug 29-30 201737

W10-2のイントロ1①cdでホームディレクトリに移動。②lsで確認。これから作成予定のbinというディレクトリがないことを確認しているだけ①②

Aug 29-30 201738

①binというディレクトリを作成して、②そこに移動

W10-2のイントロ2①②

Aug 29-30 201739

W10-2のイントロ3

multi-FASTA形式ファイルを入力として、指定した配列長未満の配列を除くPythonプログラム(fastaLengthFilter.py)を①wgetでダウンロード。このプログラムは配列長順にソートした結果を返す①

Aug 29-30 201740

①ダウンロードが無事成功したときの典型的な画面

W10-2のイントロ3①

Aug 29-30 201741

①lsで確認。②headで最初の10行分を表示

W10-2のイントロ4①②

Aug 29-30 201742

W10-2のイントロ5①chmod 755として、fastaLengthFilter.py
プログラムの実行権限を付与。②ここがXとなっていることを確認。第4回W3-1①②

Aug 29-30 201743

W10-2の本番1①取り扱いたいアセンブリ結果ファイル(polished_assembly.fasta；スライド34)のある②ディレクトリに移動して、③fastaLengthFilter.pyを実行

②

①

③

Aug 29-30 201744

W10-2の本番1①~/bin/fastaLengthFilter.pyのように書くことで、~/binに「パスを通す」作業をしていなくても実行できる。パスについては、第4回のW9-5やW15-5、第6回のW12-3など①

Aug 29-30 201745

W10-2の本番1

Aug 29-30 2017

fastaLengthFilter.pyの実行結果ファ
イルをLH_hgap.faという名前で保存46

W10-2の本番2

fastaLengthFilter.pyの①入力と②出力。③出力
ファイルのほうが若干サイズが小さくなっている②①③

Aug 29-30 201747

W10-2の本番3

multi-FASTA形式の①入力ファイルと②出力ファイ
ルのdescription行部分を表示。fastaLengthFilter.py
は②出力ファイルのdescription行部分をsequence…
のようにシンプルに記載する。この文字数の違いもファイルサイズ削減に寄与している①②

Aug 29-30 201748

W10-2の本番4①入力ファイル(polished_assembly.fasta)は40,567行。②出力ファイル(LH_hgap.fa)は8行。この改行コードのバイト数の違いがファイ
ルサイズの違いに大きく寄与している①②

Aug 29-30 201749

W10-2の本番4

赤下線部分のファイルサイズの違いは、主に①polished_assembly.fastaの塩基配列情報部分において、60塩基程度ごとに改行が入っているから。具体的には、改行コードが(40567 – 8) = 40559
個分、つまり40559 bytes分だけ大きくなっている①②

Aug 29-30 201750

W10-2の本番4①polished_assembly.fastaの2,474,245 bytesから、 (40,567 – 8) = 40,559 bytesを引くと、2,474,245– 40,559 = 2,433,686 bytesとなる。②LH_hgap.faの2,433,662 bytesと酷似しており妥当①②

Aug 29-30 201751

W10-2の本番4①polished_assembly.fastaの2,474,245 bytesから、 (40,567 – 8) = 40,559 bytesを引くと、2,474,245– 40,559 = 2,433,686 bytesとなる。②LH_hgap.faの2,433,662 bytesと酷似しており妥当。2,433,686– ②2,433,662 = 24 bytesの違いは、赤枠部分の文字数の違いに起因。③6文字×4で24 bytes

③

③
①
②

Aug 29-30 201752

あくまでもこの後の処理をやりやすくするための前処理の話。（何をもってやりやすいと思うかは
ヒトそれぞれだが）Linuxコマンドを駆使して解析を行う場合は、①LH_hgap.faのように「コンティグ
1個あたり2行」で表すほうが取り扱いやすい

ここまでの話は…

①

Aug 29-30 201753

①LH_hgap.faは全部で8行からなる。②最初の2行分がsequence1、③最後の2行分（7-8行目）がsequence4の情報となる。fastaLengthFilter.pyを実行したのは、このようにすっきりさせるため

LH_hgap.fa②③①

Aug 29-30 201754

第7回原稿p106左上

今は①赤枠のファイル分割の前処理部分の話。前処理が終わったので、②Linuxコマンドを組み合わせたファイル分割を行っていきます②①

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 201655

sequence1は、最初の2行分に相当する。それゆえ、①headコマンドで最初の2行のみ抽出した結果をsequence1.faというファイル名で保存している。それ以外の②grepや③wcはただの確認

W10-3：ファイル分割2①②③

Aug 29-30 201756

W10-4：ファイル分割2①sequence2と②sequence3は、headとtailを組み合わせて目的の配列のみ抽出。連載第3回のW19-3にもあり①②

Aug 29-30 201757

①sequence4は、tailだけでよい

W10-4：ファイル分割2①

Aug 29-30 201758

W10-5：lsで確認①ls –lで配列ごとのファイルサイズを確認。②のように書いてもよい。[0-9]は任意の数字1字という意味。③ここのファイルサイズは、配列長とほぼ同じ。④例えばsequence1の2,289,509 bytesは

①

②

④

③

Aug 29-30 201759

スライド6の①Maximum contig sizeの2,289,497 bpとほぼ同じということです。 sequence1の2,289,509 bytesは、「>sequence1」の計10文字と改行コード2文字分からなる、計12 bytes分の情報を余分に含むためと解釈

W9-2：結果を眺める①

Aug 29-30 201760

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 201761

第7回原稿p106左上①のあたりの話に移行。HGAP実行結果のFASTQファイルを眺めるところからです①

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 201662

W11-3のイントロ1①lsで、②FASTAとFASTQのサイズ比が1:2になっているので妥当①②

Aug 29-30 201763

W11-3のイントロ2

主に①FASTQファイルの行数を確認。②16行。4コンティグというのは既知。1コンティグあたり4行で表されているので、headとtailコマンドの組み合わせでどうにかなるのではと閃く①②

Aug 29-30 201764

W11-3：FASTQ分割1①W10-3やW10-4のファイル分割のやり方と若干違うのは、このような記述の仕方でもOKであることを示すとともに、シェルスクリプ
ト利用のイントロ（後述のスライド280～）①

Aug 29-30 201765

①ファイル分割前は、分割後のFASTQファイルが配列長順になっているかどうかは不明であった。②lsでsequence*のfaとfqのサイズが、③約1:2になっていることを確認して配列長順になっていると確信

W11-3：FASTQ分割1

①

②
③
③
③
③

Aug 29-30 201766

W11-3：FASTQ分割2

FASTQファイルのdescription部分の行
頭は①@および②+なので、念のため両方で調べている。③行数は全部4行

①

②

③

Aug 29-30 201767

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 201768

Tips：LinuxでR

Linux環境でのRの基本的な利用方法は、①H28年
度講習会の②8/1の③講義資料（スライド100～）①②

Aug 29-30 2017③

69

W11-4：スコア分布

赤枠部分がFASTQファイル中の1文字表記のクオリティスコアを数値化（PHREDスコアに変換）して保存するRコード。講習会用に①ShortReadパッケージをインストール済み①

Aug 29-30 201770

W11-5：入出力の関係

これは、①sequence4.fq（一番短い11,372
bpのコンティグ）を入力ファイルとして、②2つのファイルを出力するコード①②

Aug 29-30 201771

コードの枠内を全選択（Windowsのヒトはトリプル
クリックで全選択できます。CTRLキーを押しながらワンクリックでもOK）して、右クリックで①コピー

コピペ①

Aug 29-30 2017

72

（Bio-Linuxのターミナル上で）①ペースト

コピペ①

Aug 29-30 201773

W11-6：実行結果

コピペ実行結果後に①lsで確認（自動でここまでなっているはず）。②確かに指定した名前の2つのファイルが作成されていることがわかる①②

Aug 29-30 201774

Rも自動で終了し、pwdや①ls -l sequence4*まで実行されたのは、①コピペしたコードの最後のほ
うに、②のような感じで書きこんでいるからです

余談1

①

②

Aug 29-30 201775

もちろんペーストは、①ShortReadパッケージがインストールされていれば、ホストOS（この場合はWindows10）上のRでやっても構いません。これが共有フォルダの便利なところ

余談2①

Aug 29-30 201776

W11-7：pngファイル①pngファイルのほうは、②で横幅700 pixelと縦幅
350 pixelにサイズ指定しているので横長になっている。横軸はコンティグのposition、縦軸はPHRED
スコア。数値が大きいほどクオリティが高い。尚、Linux上で出力ファイル（sequence4.png）を開く場
合は、左側の引き出しアイコンから辿っていくのでもよいし、「eog sequence4.png」と打つのでもよい①

350 pixel②

Aug 29-30 2017

700 pixel

77

W11-8：テキストファイル

一番短い11,372 bpのコンティグなので、①sequence4.txtは11,372行×2列になる①

Aug 29-30 201778

①最初の5行分と②最後の4行分を表示。③1列目はposition番号、④2列目がPHREDスコア

W11-8：テキストファイル

①

②

③

Aug 29-30 2017④

79

W11-9：sequence3.fq①2番目に短い45,853 bpのファイル(sequence3.fq)に対しても同様な作業を実行。
②pngファイルを眺めると、確かに③コンティグ
両末端部分のクオリティが低いことがわかる①②③
Aug 29-30 2017③
80

クオリティスコアが0なのは、最初から①1,000塩基目
ちょっとと、②最後から1,500塩基ちょっと（44,000塩基目）くらいかな…とか妄想し、ある程度あたりをつけておいて、③sequence3.txtをlessで眺めたりしておくW11-9：sequence3.fq③①
Aug 29-30 2017②81

①sequence3.txtをlessで眺めようとしている

less sequence3.txt①

Aug 29-30 201782

less sequence3.txt①

Aug 29-30 2017②

less実行直後。①1列目がposition番号（塩基番号）なので、ここが1,000ちょっとくらいにな
るまで、ざばーっと「スペース」キーで進み、② PHREDスコアが0より大きくなるところを探索する。残りは上下矢印キーで微調整83less sequence3.txt①1,224番目の塩基以降から0以上のPHREDスコアになっているようですね①

Aug 29-30 201784

less sequence3.txt

Aug 29-30 2017

「G」と打って、最終行に移動したところ。この状態からキーボードのuを押して効率的にページ上部を探索する

85less sequence3.txt①44,518番目の塩基までが0以上のPHREDスコアになっているようですね①

Aug 29-30 2017

86

qと打って、lessから抜けたところ

less sequence3.txt

Aug 29-30 201787

①の話までが終了しました。赤下線部分が重要第7回原稿p106左中①

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 201688

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 201789

第7回原稿p106左下①コンティグが環状になっているかどうかを調べる基本的な戦略が②ドットプロット。③
単純な塩基配列のFASTAファイルをLinux上で作成するところ（W12-4）から解説します

①

②

③

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 201690

W12-4：hoge.fa①

seqinrパッケージのdotPlot関数実行時に入力として用いるファイルhoge.faを作成する。
①とりあえずここで作業。②と③は作成した
いhoge.faがないことを確認しているだけ②③

Aug 29-30 201791

W12-4：hoge.fa①echoコマンドは、第4回W9でも利用。ここでは任意の文字列「>test」を表示させているだけ。
②リダイレクト(>)でhoge.faファイルを新規作成
して書き込み。③moreでhoge.faの中身を表示①②③

Aug 29-30 201792

①ACTCGTCAGAという文字列をhoge.faに追加書き込み。②既存のファイルに追加書き込
みをするときは>>でした（第4回W12-2；第5回
W17-2）。③これでFASTA形式ファイルの完成

W12-4：hoge.fa②①③

Aug 29-30 201793

W12-5：dotPlot実行①赤枠内がドットプロット（300×300ピクセ
ルのhoge1.pngファイル）を作成するコード①

Aug 29-30 201794

W12-5：dotPlot実行

コピペ実行後に①lsで確認。②hoge1.pngが確かに作成されている。③中身はこれ。「eog hoge1.png&」
で開くことができます（スライド77）③①②

Aug 29-30 201795

ACTCGTCAGA

W12-6：解説

ドットプロットの解説。Rのseqinrパッケージ中のdotPlot
関数実行結果ファイルは、①左下を原点として比較す
る2つの配列を並べている。一致が黒、不一致が白A C TGGCA
A CT① ACTCGTCAGA

Aug 29-30 201796

同一の配列を比較するときは、①必ず対角線上の塩基が一致（つまり黒）する

W12-6：解説

ACTCGTCAGA①

A C TGGCA
A CT

ACTCGTCAGA

Aug 29-30 201797

W12-7：環状コンティグ例

アセンブリ結果として、①最初と②最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断

ACTCGTCAGAA C TGGCA
A CT

ACTCGTCAGAACTC①

Aug 29-30 2017②

98

W12-7：環状コンティグ作成

両末端の4塩基がACTCで同じ、
計14塩基からなるミニ環状コンティグファイル(hoge2.fa)を作成。最初の10塩基分は、W12-4で作成したhoge.faと同じ

ACTCGTCAGAACTCAug 29-30 201799

ACTCGTCAGAACTC

W12-8：ドットプロット①両末端の4塩基がACTCで同じ、計14塩基からなるミニ環状コンティグファイル(hoge2.fa)を入力として、再度ドットプロットを実行した結果。赤枠以外がACTC追加部分

A C TGGCA
A CT

ACTCGTCAGAACTC①

Aug 29-30 2017①

100

W12-8：環状の場合
ACTCGTCAGAACTC②

こんな感じに見えます。①対角線上にプロットされるのは同じですが、②対角線と平行に末端部分もプロットされるのが環状の特徴①

A C TG
G
C
A
A CT
②

ACTCGTCAGAACTCAug 29-30 2017101

ACTCGTCAGAACTC

W12-8：環状の場合

これは、コンティグの両末端が同じ配列であることを意味する。「重複する部分の除去」は、「complete
genomeにする操作(finishing)」の一部に相当する

A C TGGCA
A CT

ACTCGTCAGAACTCAug 29-30 2017102

ACTCGTCAGAACTC

W13-1：重複除去

重複除去（トリミング）の選択肢は、（この場合は結果的に同じになるが）①5通り存在する。通常（推奨）は、両末端はクオリティが低いので、②中央部分を残して両端をトリムする選択肢（*のついたもの）を採用する

A C TGGCA
A CT

ACTCGTCAGAACTCAug 29-30 2017

ACTCGTCAGAACTCACTCGTCAGA
CTCGTCAGAA ②
TCGTCAGAAC*CGTCAGAACTGTCAGAACTC①

103

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017104

第7回原稿p106右中①簡単な例でLinuxコマンドを用いた実際の重複除去の例を示します①

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 2016105

W13-2：cutコマンド①

特定の範囲の切り出しはcutコマンドを利用。①
hoge2.faはsingle-FASTA形式。②「tail –n 1」で、最後の1行分のみ取り出している。③パイプで流してcutコマ
ンドを実行し、3-12文字目を表示。これが④環状コンテ
ィグの重複除去後の塩基配列に相当する。これは1行1配列にしているからこそできるワザ！（60塩基程度ごとに改行が入っている通常のFASTA形式だとムリ！）②③

ACTCGTCAGAACTCACTCGTCAGA
CTCGTCAGAATCGTCAGAAC* ④
CGTCAGAACTGTCAGAACTCAug 29-30 2017106

W13-3：続cutコマンド①最初の10塩基分のみ取り出しても、②
最後の10塩基分のみ取り出しても、この
シンプルな環状ゲノムの場合は③アスタ
リスクのついたものと結果的には同じ

①
①
②
②

ACTCGTCAGAACTC① ACTCGTCAGA
CTCGTCAGAATCGTCAGAAC* ③
CGTCAGAACTGTCAGAACTC ②

Aug 29-30 2017107

理由はスタート地点をどこにするかという違いのみだから。本物の環状染色体の場合は、特定の遺伝子配列が先頭になるように回転させる慣例がある

W13-3：続cutコマンド

A C TG
G
C
A
A CT
①

③

ACTCGTCAGAACTC① ACTCGTCAGA
CTCGTCAGAATCGTCAGAAC* ③
CGTCAGAACTGTCAGAACTC ②

Aug 29-30 2017108

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017109

第7回原稿p106右下①このあたりの話に移行します。seqinrパ
ッケージのdotPlot関数で行うのが難しい
ところはすっ飛ばします。気になる人は
第7回W14-1を自分でやってみてください①

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 2016110

W14-2：dotter①sequence3.fa同士のドットプロットをdotterで実行。画面はリターンキーを押して数秒後の状態①

Aug 29-30 2017

Sonnhammer and Durbin, Gene, 167: GC1-10, 1995

111

W14-2：dotter①

計3つのウィンドウが立ち上がる。①Greyramp Tool
はよくわかりませんが、ドットプロットのコントラスト
調整用なのだろうと思います。②のあたりをクリック
して、アラインメントのウィンドウをアクティブにする②

Aug 29-30 2017112

W14-2：dotter①①Alignment Toolは、比較している2つの配列のアラインメント結果を表示。②今比較しているのは同じ配列なので完全一致。③RevCompと書いてあるので、④片方をReverse Complement（逆相補鎖）にした結果も表示されていることがわかる②③

Aug 29-30 2017④

113

W14-2：dotter

裏側に見えているのが主目的のドットプロット。①このあたりでクリックして、ドットプロットのウィンドウを手前に表示①

Aug 29-30 2017114

W14-2：dotter

赤下線部分に7671, 16576という数字が見えている。これは今比較している2つの配列の塩基番号（①の座標情報）に相当する。同じ位置をクリックしなければいけないわけではなく、ドットプロットのウィンドウが手前に表示されていればOK。
②かすかに環状の特徴が見えている。バグってログイン画
面になりがちだが（爆）、ウィンドウを広げて全体像を眺める②①

Aug 29-30 2017115

W14-3：バグった例

Aug 29-30 2017

バグってログイン画面になっても、気を取り直して再挑戦しよう。右下のようなドットプロットが得られるはず116

W14-4：ドットプロット全体像①

dotterのドットプロットは①左上が原点のようだ。②右上（と左下）にもくっきりと「①の対角線」と並行のラインが見えているので、45,853 bpのsequence3は環状と判断②

Aug 29-30 2017117

W14-4：ドットプロット全体像①

大まかには、全部で45,853 bpのうち、①最初と②最後の約5,000 bpが重複していると判断②

Aug 29-30 2017118

W14-5：dotter終了

Dotter実行結果を終了させるには、基本的に①該当するGUIの左上部分の×を押せばよい①①

Aug 29-30 2017119

W14-5：dotter終了②

Aug 29-30 2017

①赤い点線の枠内にカーソルを移動させるとメニューバーが見られるようになるので、②×。第3回W6-3①

120

W14-6：dotter終了後①コマンド打ち込み可能状態になっていることがわかります①

Aug 29-30 2017121

第7回原稿p106-107①dotter出力結果のドットプロットを眺めて
全体の傾向を把握する。次にBLASTを実行して正確なアラインメント情報を得る。クオリ
ティスコア分布情報（W11-9）とアラインメント
情報を用いて、重複除去の本番を行う①

Aug 29-30 2017

谷澤ら, 日本乳酸菌学会誌, 27: 101-110, 2016122

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017123

第7回原稿p107の左上①BLASTを（トリッキーな使い方だが）sequence3.faの重複領域の正確な情報を得る目的で利用する①

Aug 29-30 2017

Altschul et al., J Mol Biol., 215: 403-410, 1990

124

W15-1：makeblastdb①作業ディレクトリはどこでもよい。②makeblastdbのバージョンは2.2.28+①②

Aug 29-30 2017125

W15-1：makeblastdb①「-h」で大まかな利用法(usage)を確認。
②より詳細な説明は「-help」で出るようだ①②

Aug 29-30 2017126

①makeblastdb本番。入力はsequence3.fa。塩基配列であることを示すnuclを-dbtypeオプションで指定

W15-1：makeblastdb①

Aug 29-30 2017127

W15-1：makeblastdb①実行後は、sequence3.fa.n*というファイルが8個作成されている②

Aug 29-30 2017128

①blastnを実行。②DB側、③query側はともにsequence3.fa。④出力ファイル名はsequence3_blast.txt。計算は一瞬

W15-2：blastn

②

③

①
④

Aug 29-30 2017129

W15-2：blastn①lsで確認。確かに②出力ファイル(sequence3_blast.txt)があります①②

Aug 29-30 2017130

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017131

W15-3：結果を眺める

blastn実行結果ファイルを眺めるべく、①sequence3_blast.txtの最初の10行分（デフォルトが10行）を表示。②行数は3,852行①②

Aug 29-30 2017132

W15-3：結果を眺める

sequence3_blast.txtを適当なテキストエディタ
で概観。①1つ前のスライドの赤枠と同じ部分①

Aug 29-30 2017133

W15-3：結果を眺める①DB側の配列はsequence3.fa、②query側の配列もsequence3.fa。長さは45,853 bp①②

Aug 29-30 2017134

W15-3：結果を眺める①トップヒットのアラインメント結果の概要。
トップヒット(top hit)はsequence3の全長配列間で100%一致のアラインメントなので、②「Identities = 45853/45853 (100%)」②①

Aug 29-30 2017135

W15-3：結果を眺める

トップヒットのアラインメント結果は、①
dotterによるドットプロット（W14-4）の左上から右下にかけての対角線に相当します①

Aug 29-30 2017136

W15-4：ヒット数①BLAST結果ファイル中のヒット総数を把握したい場合は、”Score =”という文字列を含む行数をgrepで調べればよい①

Aug 29-30 2017137

W15-4：ヒット数①”Score =”という文字列を含む行を表示。
②その行数は10個。つまりヒット数は10①②

Aug 29-30 2017

138

参考

W15-5：grep -n①grep実行時に-nをつけることで検索文字列（この場合”Score =”）を含む行番号を表示。例えば
セカンドヒットは3,092行目、②サードヒットは3,425
行目などというのがすぐにわかる。③BLAST結果
ファイル(sequence3_blast.txt)は3,852行だった①②③

Aug 29-30 2017139

参考 ①grep –Aオプションで一致した行を含め後ろの3行分を表示。②1位は、一致領域が45,853 bpで100%一致、Gapsも0/45853。W14-4のドットプロットの③対角線の見栄えと一致

W15-6：grep -A①②③

Aug 29-30 2017

140

参考

W15-6：grep -A

同率2位の①と②の2つのヒットは、一致領域が4,901 bpとなっており、W14-4のドットプロット上の2つの赤四角の見栄えと一致する。どちらがどちらに相当するかは、
W15-5で得られたセカンドヒットの3,092行目以降、サード
ヒットの3,425行目以降のアラインメントを眺めればわかる①②

Aug 29-30 2017141

第7回原稿p107の左中①今はこのあたりの話。これから②のアラインメント結果を詳細に眺め、重複除去ポイントを探る①②

Aug 29-30 2017142

W15-7：lesslessコマンドでsequence3_blast.txtを開き、Score =で検索。 画面の横幅を広めにとっておいたほうがよい。第3回のW14-6-2に文字列検索のやり方あり①

Aug 29-30 2017143

lessコマンドでsequence3_blast.txtを開いた直後

W15-7：less

Aug 29-30 2017144

①「/Score =」と打って、Score =という文字列を検索

W15-7：less①

Aug 29-30 2017145

①トップヒットのものが最初に見える。②全長の45,853
bp全てで完全一致なので、③queryの1-60番目の塩基とDB側(Sbjct; Subjectの意味)の1-60番目の塩基だけで眺めても完全一致となっていることがわかる

W15-7：less①②③

Aug 29-30 2017146

W15-8：less「n」と打って、2番目に一致するScore =が先頭行にくる
ページを表示した結果。①query配列の40,967番目の塩基がDB側配列の1番目の塩基と一致していることを意味する。nとは逆方向に検索していきたい場合はN①

Aug 29-30 2017147

W15-8：less①query配列の40,967番目の塩基は、②のあたりのポジションに相当します②①

Aug 29-30 2017148

W15-8：less②や③のように、DB側(Sbjct)のところどころでGapが見られる。が、④全体で4,901 bpのアラインメントのうち31個だけGapがあった程度なので、実質的に無視でよい

④

②

③

Aug 29-30 2017149

W15-8：less「n」と打って、3番目に一致するScore =が先頭行にくるページを表示した結果。 ①query配列の1番目の塩基がDB側配列の40,967番目の塩基と一致していることを意味する。2番目と3番目はQueryとSbjctが入れ替わっているだけ
④

②

①

③

Aug 29-30 2017150

上矢印キーを10回押し、10行分だけページ上部に移動した結果画面。セカンドヒットのアライ
ンメント結果の最後のほうを確認するのが目的

W15-9：less

④

②

①

Aug 29-30 2017151

W16-1：トリム候補領域

上矢印キーをさらに押し続け、 （重複塩
基数が4900 bp程度なのでその半分の）2400 - 2500番目付近を眺める。具体的には①の赤枠分くらいを眺め、どこにもミスマッチやGapがないことを確認①

Aug 29-30 2017152

W16-1：トリム候補領域①のところでトリムすることにする。左端にする理由は、上が2450番目、下が43423番目の塩基だとすぐにわかるから①

Aug 29-30 2017153

W16-2：トリム後の配列①2450番目の塩基をトリム後の1塩基目にする場合は、 [2450, 43422 bp]を残せばよい。こうすることで、トリム後の塩基配列の最初のほうは①の赤枠のようになり、最後のほうは②のようになるはずである。③qで終了②①③

Aug 29-30 2017154

第7回原稿p107の右中

今は①のあたり。1つ前のスライドで決めた領域のトリム（重複除去）を行う。混乱してきたら、W13-1を復習。このあと②トリムを実行します(W16-3)①②

Aug 29-30 2017155

W16-3：トリム実行①まずはトリム後のFASTAファイル（ファイル名：sequence3_trimmed.fa）のdescription行を作成。W12-7とほぼ同じ①

Aug 29-30 2017156

W16-3：トリム実行①トリム実行本番。W13-2とほぼ同じ。②sequence3.faの最終行のみ取り出してパイプで流し、③2450-43422文字目を抽出した結果を、④sequence3_trimmed.faに追加書き込み①②

Aug 29-30 2017

③

④

157

W16-4：moreで確認①lsで確認。1 bp = 1 byte。ファイルサイズ
的に妥当な印象を受ける。②moreでも確認①②

Aug 29-30 2017158

W16-4：moreで確認①赤枠で示すトリム後の塩基配列の最初のほうは、W16-2と全く同じになっていることからうまくトリムできたと判断①

Aug 29-30 2017159

W16-4：moreで確認

スペースキーをガスガス押して最後まで表示し終わったところ。②最後の塩基配列の赤枠部分もW16-2と全く同じになっていることからうまくトリムできたと判断②

Aug 29-30 2017160

Contents（第7回後半分）
イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017161

W17-1：FASTQのトリム①FASTQファイル(sequence3.fq)の場合は、2行目(塩基配列情報の行)と4行目(クオリティ情報の行) についてのみheadとtailを組合せた操作（W16-3）を行えばよい。②得られるファイルはsequence3_trimmed.fq①②

Aug 29-30 2017162

W17-2：クオリティ分布

Aug 29-30 2017

FASTQファイル(sequence3_trimmed.fq)
を入力として、W11-9（スライド80）と同じようなクオリティスコア分布を作成

163

W17-2：クオリティ分布①コピペ実行後に得られるファイル。「eog sequence3_trimmed.png &」
で開けます

①

Aug 29-30 2017164

W17-2：クオリティ分布①pngファイルを眺めているところ。W11-9で見られていた両側の低クオリティ領域がうま
くトリムされていることがわかる。図3aと同じ①

Aug 29-30 2017165

W17-3：ドットプロット①sequence3.fa同士のドットプロットをdotter
で実行。トリム前（W14-4）と違って、配列末端部分の一致領域（②や③付近の対角線のプロット）がなくなっていることがわかる。うまく重複領域をトリムできていることを意味する①②③
Aug 29-30 2017

Sonnhammer and Durbin, Gene, 167: GC1-10, 1995

166

ここまでのまとめ

PacBioデータのde novoアセンブリ結果フ
ァイルをもとに、3番目に長い配列（約4.6万塩基のsequence3.fa）の重複配列除去を行った。時間はかかるが、2番目に長い配列（約8.7万塩基のsequence2.fa）も同じ要領で可能なので、トライしてみてください

イントロダクション（主に予習事項の確認）
 W10：multi-FASTAファイルの分割
 W11：FASTQファイルの分割とクオリティスコア分布
 環状化（ゲノム解読のfinishing作業の一部）
 アセンブリ結果として、最初と最後の末端部分が同じ配列の場合は、通常そのコンティグは環状と判断。それを確認するための基本的な考え方、手段、および環状化のノウハウを伝授 W12：seqinrパッケージを用いて、仮想環状コンティグのドットプロットで感覚をつかむ W13：重複配列の除去の感覚をつかむ（これが環状化作業の実体）
 W14：dotterプログラムで実際のPacBio出力結果に対して適用し、環状状態の概要を知る
 W15：Bio-Linux上のblastnで、環状の同一コンティグ同士をDB側とquery側にして実行 W16：正確なアラインメントを眺め、切断箇所をクオリティスコア分布と合わせて判断する
 W17：両端切断後のコンティグに対し、クオリティスコア分布やdotterを再度実行して確認

Aug 29-30 2017167

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017168

第8回原稿p188の左中

今は①のあたり。②でも書いているが、本番ではなく予備的なアノテーション①②
Aug 29-30 2017

169

参考

W4-1：DFAST

DFASTは乳酸菌に特化したゲノムアノテーションを行ってくれるウェブツール。入力はmulti-FASTAファイル。出力はGFF形式ファイルなどのアノテーション結果。
DFASTのトップ画面上で、①をクリック。となってはいるが、バクテリア（しかも乳酸菌のみ）に特化している
ため、エアーハンズオン（やったふり）でアノテーション
ツール利用法の概要を把握する程度で十分です①

Aug 29-30 2017

Tanizawa et al., Biosci Microbiota Food Health, 35: 173-184, 2016170

W4-1：DFAST①のところで、参照ボタンを押して、アノテーション
したいmulti-FASTAファイル(LH_hgap.fa)を指定①

Aug 29-30 2017171

W4-2：LH_hgap.fa①アノテーションしたいmulti-FASTAファイル(LH_hgap.fa)を共有フォルダにコピー。②共有フォ
ルダ内の状況はヒトそれぞれだが、最低限
LH_hgap.faが見えていればよい。（私のWindows環
境では）ホストOS上の絶対パスとして③のところに
LH_hgap.faが見えるので、これをアップロードします①②③

Aug 29-30 2017172

W4-3：アップロードと実行①参照ボタン、②に移動して、③目的のファイル(LH_hgap.fa)を指定して、④開く

②

③
①

④

Aug 29-30 2017173

①ジョブタイトル(hogegeee)をテキトーにつけて、②計算が終わったらメールでお知らせしてくれるので入力して、③ページ下部に移動

W4-3：アップロードと実行

①

③

②

Aug 29-30 2017174

①いろいろオプション指定できるがとりあ
えずここは無視して、②Run。③は設定値(デフォルトは200 bp)以下の短い配列
を除くため。それ以外の赤枠の属・種名等のオプションはアノテーション結果に影響を与えることはなく、後で変更することも可能であるため、規定値のままでOK

W4-3：アップロードと実行①③②
Aug 29-30 2017175

W4-3：アップロードと実行

Aug 29-30 2017

計算が始まったようだ。とりあえず計
算終了メールが来るまで思考停止

176

W4-4：計算終了

計算終了メールが届く。このときは、約5分で計算が終了。①にアクセス①

Aug 29-30 2017177

計算結果が見えている。①ページ下部に移動

W4-5：結果を眺める①

Aug 29-30 2017178

W4-5：結果を眺める

これがDFASTによるアノテーション結果の全体像。①このあたりの数値は、DDBJ Pipeline上でHGAPを実行したときの結果（第7回W9-2）と全く同じで妥当①

Aug 29-30 2017179

W4-5：結果を眺める①アノテーションファイルの一般的な形式である
GFFファイルをダウンロードしてエクセルなどで眺めてもよいが、ここでは②Featuresタブをクリックして、ウェブ上でアノテーション結果を眺める②①

Aug 29-30 2017180

W4-5：結果を眺める

こんな感じになります。①デフォルトはページあたり25エントリーになっているので、②の部分が25行分づつ表示されます①②

Aug 29-30 2017181

W4-5：結果を眺める

全エントリーを一気に表示させて説明したいので、①Allにする①

Aug 29-30 2017182

W4-5：結果を眺める

こんな感じになります。横幅が広くなるが細かいことは気にしない。①のあたりを見ることで、入力ファイル(LH_hgap.fa)の配列および座標順にアノテーションされた結果が表示されていることがわかる①

Aug 29-30 2017183

①と②を関連づけながら、アノテーション結果の全体像を把握する

W4-6：sequence1概観①

Aug 29-30 2017②

184

W4-6：sequence1概観①sequence1の左端（最初の3031 bpまで）
は「hypothetical proteinが多いなあ…」とか、②「prophage proteinがある」とか…

①②

Aug 29-30 2017185

W4-6：sequence1概観

赤枠はsequence1の①19649から②34418
bpの範囲（全部で2,289,497 bpの長さがあるので、このあたりもまだ左端といえる）。
赤下線で示すように、このあたりにもファージ(phage)関連のものがちらほら存在①②
Aug 29-30 2017186

W4-7：sequence2概観①ここからが、sequence2のアノテーション結果①

Aug 29-30 2017187

W4-7：sequence2概観

第7回でsequence2は環状のプラス
ミド配列だろうと予想していたが、①それを補強するアノテーション結果①

Aug 29-30 2017188

W4-8：sequence3概観①ここからが、sequence3のアノテーション結果①

Aug 29-30 2017189

W4-8：sequence3概観①赤下線で示す接合伝達(conjugaltransfer)関連遺伝子が多く見られることからも、sequence3がプラスミ
ドであることを裏付けている①

Aug 29-30 2017190

W4-9：sequence4概観①ここからが、sequence4のアノテーション結果①

Aug 29-30 2017191

①ここにtransposaseがあります

W4-9：sequence4概観①

Aug 29-30 2017

192

第8回原稿p188の右上①までで予備的なアノテー
ション結果の概観が終了①

Aug 29-30 2017193

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017194

①dotterでsequence4同士を比較

W5-2：dotter (sequence4)①

Aug 29-30 2017195

①この結果は、 (1 – 500)番目と(750 – 1,350)番目付
近の領域が似ていることを意味する。しかし両末端ではないので、sequence4は環状（プラスミド）ではないのだろうと判断。sequence2同士は第8回W5-1にあります（時間はかかるが実行可能）。sequence1同士は
第8回W5-3でトライしています（事実上実行不可能）

W5-2：dotter (sequence4)①

Aug 29-30 2017

196

第8回原稿p188の右下

今は①のあたり。sequence1をquery側の配列、sequence1-4をDB側の配列としてBlast実行する話へと移行していきます①

Aug 29-30 2017197

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017

198

第8回原稿p188の右下①sequence1をquery側の配列、
sequence1-4を含むLH_hgap.faをDB側の配列として、Blast実行する話へと移行。
②Blastだといっぺんに調べられて便利①②

Aug 29-30 2017199

W6-1：makeblastdb①makeblastdbの実行。第
7回W15-1と基本的に同じ①

Aug 29-30 2017200

W6-2：blastn①DB側をLH_hgap.fa、query側をsequence1.fa、出力ファイルをsequence1_blast.txtとしてblastnを実行。②出力ファイルは12MB。sequence1.fa同士のアラインメント結果を丸々含んでいるためであろう①②

Aug 29-30 2017

201

ヒット数

①”Score =”という文字列を含む行数は1,351個。つまりヒット数は1,351。これだけ多いと全体像の把握が困難なので、Blast結果閲覧専用のViewerを利用する①

Aug 29-30 2017202

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017203

第8回原稿p189の左上①Blast結果閲覧用Viewerの1つであるBlastViewerで説明していきます②

Aug 29-30 2017204

①Korilogというところが提供している
BlastViewerのウェブページ。②Windows版と③Macintosh版。ときどき②と③の上下が入れ替わるようなので注意。次のスライド以降はWindows版で説明

W7-2：BlastViewer

①

②

③

Aug 29-30 2017205

W7-2：BlastViewer

Windows 10上でのインストール画面。①ここに作成されるようだ。②Install。③デスクトップにBlastViewerのアイコンが作成されるはず。
④Finish。デフォルトでは⑤にチェックが入っているのでBlastViewer (ver. 2.2)が起動する

⑤

①

②
④
③

Aug 29-30 2017

206

Java関連でエラーが出たら①

BlastViewer (ver. 2.2)が起動できず、①Java関連でエラーが出たら、③どうにかして自力でインストール
してください。アグリバイオ貸与PC
はJavaインストール済み（なはず）②③

Aug 29-30 2017207

W7-2：BlastViewer①を見ることで、Blast結果ファイルを開く
際には、②を押せばよいことがわかる。
③をよく見るとBlast出力結果ファイルはXML形式しか受け付けていないようだ

②

①

③

Aug 29-30 2017208

W7-2：BlastViewer①ためしにW6-2で作成したテキスト形式のBlast結果ファイル(sequence1_blast.txt)
を読み込もうとしたらダメでした。②OK①②

Aug 29-30 2017

209

①XML形式のBlast結果ファイル(sequence1_blast.xml)を作成。②確かにあります

W7-3：Blast再実行①②

Aug 29-30 2017210

W7-4： BlastViewer①①共有フォルダ上の、②XML形式のBlast結
果ファイル(sequence1_blast.xml)を、③開く③①②

Aug 29-30 2017211

W7-4： BlastViewer①共有フォルダ上の、②XML形式のBlast結
果ファイル(sequence1_blast.xml)を、③開く①②③

Aug 29-30 2017212

W7-4： BlastViewer①

読み込み後の状態。①query側がsequence1.fa、②DB側がsequence1-4からなるLH_hgap.faだったことを思い出そう②

Aug 29-30 2017213

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017214

W7-5：解釈1①は、DB側のヒット数（配列類似領域の数；HSPの数)が、sequence1中に1347個、sequence4中に2個、sequence3中に1個、sequence2中に1個あったことを示す。sequence1-4の並びではないので、HSP数か
スコアの高い順でソートされているのだろうと妄想①

Aug 29-30 2017215

W7-5：解釈2①テキスト形式のBLAST結果ファイル(sequence1_blast.txt)
をlessで眺めて全体像を大まかに把握①

Aug 29-30 2017216

W7-5：解釈3

lessで開いた直後の状態。①用いたBLASTNのバージョン、②DB側、③query側の情報

①

②

③

Aug 29-30 2017217

W7-5：解釈4

1ページ分ほど下に移動(上下左右の矢印
キーで移動)。①確かにBlastViewerと同じ並び(sequence1, 4, 3, 2)になっている①

Aug 29-30 2017218

W7-5：解釈5①これがトップヒットの基本情報。query (sequence1)とDB(sequence1-4)の関係が分かっていれば、これがsequence1
vs. sequence1の100%一致の結果であることがわかる。スコ
アの計算方法がよくわかっていなくても、大まかに配列長(2,289,497 bp)の2倍程度の値がスコアっぽいなと学習する①

Aug 29-30 2017219

W7-5：解釈6となると、いくら①E-valueがそこそこ低くても、②sequence2にヒットしているsequence1の断片配列の領域は、スコアが124なので60塩基程度だろうと予想。また、③sequence3にヒットしているsequence1の断片配列の領域は、スコアが2113なので1000塩基程度だろうと予想①②

Aug 29-30 2017③220

W7-5：解釈7①sequence2でキーワード検索し、②のアラインメント結果を眺めます②①
Aug 29-30 2017221

最初はこんな感じになります

W7-5：解釈8

Aug 29-30 2017222

W7-5：解釈9

今は、テキスト形式のBLAST結果ファイル(sequence1_blast.txt)
をlessで眺めている。①sequence2でキーワード検索（nで順方向に検索、Nで逆方向に検索）し、アラインメントを表示させているところ。②ヒットしている領域は76塩基①②

Aug 29-30 2017223

①sequence3でキーワード検索

W7-5：解釈10①
Aug 29-30 2017224

W7-5：解釈11

Aug 29-30 2017

最初はこんな感じになるかもしれないが、「nで順方向に検索、Nで逆方向に検索」なので、Nと打ってページ上部を探すと…

225

W7-5：解釈12

こんな感じになって、①sequence3のアラインメント結果を見ることができます。②ヒットしている領域は1276塩基①②

Aug 29-30 2017226

W7-6：grep

トータルのヒット数が(1347 + 2 + 1 + 1) = 1351個だったということは、①「Score = 」を含む行数も1351個あるのだろうと予想（スライド215）。grepで確認すべく、qを押してlessから一旦抜ける①

Aug 29-30 2017227

W7-6：grep①「Score = 」を含む行数をgrepで調査。確かに1351個ある（スライド202）①

Aug 29-30 2017228

W7-6：grep①「Score = 」を含む最初の10行分を表示。②のように同じスコアのものが2個づつ表示されている状態を見て、「sequence1も両末端に重複領域がある環状コンティグなのだろう」と予想する

①

②
②
②
②

Aug 29-30 2017229

W7-6：grep①「Score = 」を含む最初の3行分を表示。②grep –Aオプションで一致した行を含め後ろの3行分を表示。③が今注目しているところ

①

②

③

Aug 29-30 2017230

W7-6：grep -A①Plus/Plusになっているので、おそらく両端に重複領域があるのだろう。sequence3同士のBlast実行結果（
第7回W15-6）も両端でPlus/Plusになっていたことを思い出せば納得できるだろう。本当に両端かどうかはBlastViewerなり、このテキスト形式のBlast結果ファイル(sequence1_blast.txt)を詳細に眺めて確認する①①

Aug 29-30 2017231

W7-6：grep -n①grep –nで検索文字列（この場合”> sequence”）を含む
行番号を表示。例えば②sequence4とヒットしたという結果情報は、入力ファイル(sequence1_blast.txt)の205785
行目からスタートしていることがわかる。第7回W15-5①②

Aug 29-30 2017232

①grep –nで検索文字列（この場合”Score =” ）を含む行番号も表示。②赤枠の上下関係（sequence1の最後の③205776とsequence3の④206509）から、sequence4とヒットした領域は2か所あり、いずれもスコア（約10320と8907）がそこそこ高い。スコアの合計（2万弱）とsequence4の長さ(11,372 bp)の関係から、sequence4の大部分の領域がsequence1と類似していると判断できる

W7-6：grep -n

①
③
②
④

Aug 29-30 2017233

行番号の関係性から、（queryであるsequence1と）①sequence3の一致が②のスコアであり、（sequence1と）③sequence2の一致が④のスコアであることがわかる。こんな感じでLinuxコマンドオ
プションを駆使して、BlastViewerのような可視化
ソフトで表示される内容の確認を（私は）行います

ちなみに

①

②

Aug 29-30 2017

③

④

234

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017

235

第8回原稿p189の左上①今はこのあたり。だんだんマニアックな話になっていくが、Blast結果とアノテーシ
ョン結果を眺めることで②のようなことが判明することもある、といったことをご自身の研究に活かせると…いいですね①②

Aug 29-30 2017236

①この話をします

第8回原稿p189の左中①

Aug 29-30 2017237

W8-1：seq1 vs. seq4

BlastViewerに戻り、①を押してsequence4とのヒット領域を眺める。ここでは、「sequence1（query側の配列）のどのあたりが、sequence4（DB側の配列）のどのあたりと一致しているか」を眺めようとしています。思考停止禁止!!

①

Aug 29-30 2017238

W8-1：seq1 vs. seq4

大まかな見方を説明。①3という数字は、②Accessionのところの数字と同じものであり、3’末端などという意味ではない。
③DB側配列(sequence4; 11,372 bp)、④
query側配列(sequence1; 2,289,497 bp)②④①

Aug 29-30 2017③

239

W8-1：seq1 vs. seq4

①DB側配列(sequence4; 11,372 bp)と、②query側配列(sequence1; 2,289,497 bp)は実際の長さは異なるが、長さを揃えて表示している②①

Aug 29-30 2017240

W8-1：seq1 vs. seq4①(sequence4; 11,372 bp)中には、配列類似領域
が2つ（HSPが2個）あったことを思い出そう。②をクリックして最初に見ているのは、③スコアが高いほうの「Score = 1.032e+04 bits (5588)」のHSP。④対応するquery側のアラインメント領域

①

②

④
③

Aug 29-30 2017241

W8-1：seq1 vs. seq4

①2つあるHSPのうち、②最初の1個目という意味。③次のHSPのアラインメントが見られる①①③
Aug 29-30 2017②

242

W8-1：seq1 vs. seq4

次のHSP（「Score = 8907 bits (4823)」）のアラインメント。①赤枠の位置が変わっているのがわかる。②対応するquery側のアラインメント領域もわずかに右側にシフトしていることがわかる②①

Aug 29-30 2017243

W8-1：seq1 vs. seq4

これは、①query側配列(sequence1)が②左から右の方向(Plus鎖)に並んでいるのに対して…

①②

Aug 29-30 2017244

W8-1：seq1 vs. seq4

これは、①query側配列(sequence1)が②左から右の方向(Plus鎖)に並んでいるのに対して、③DB
側配列(sequence4)が④右から左の方向（Minus
鎖）に並んでいることからも納得できる。テキスト
形式のBlast結果ファイル(sequence1_blast.txt)を詳細に眺めて確認してもいいだろう

①

③
②
④

Aug 29-30 2017245

①grep実行時に②-nと③-A 3を組み合わせて、“Score = ”と一致した行を含めた後ろの3行分を行
番号とともに表示させて、もう少し詳細情報を表示

W8-2：seq1 vs.
seq4
② ③
①

Aug 29-30 2017246

sequence4上の配列類似領域の①1つめと②2つめ。
③確かにストランド情報として、query側(sequence1)
がPlus鎖、DB側(sequence4)がMinus鎖となっている

W8-2：seq1 vs. seq4

①

③
②
③

Aug 29-30 2017247

W8-3：HSP1①1つめのHSP (以下、HSP1)に切り替えて、②Alignmentタ
ブをクリック。ここでアラインメント結果の詳細情報がわかる②①
Aug 29-30 2017248

W8-3：HSP1①query側(sequence1)の549,494番目と、②DB側(sequence4)の11,372番目から…③を一番右まで移動

①
②

①

②

Aug 29-30 2017③249

①query側(sequence1)の555,171番目と、②DB側(sequence4)の5,706番目の領域がHSP1。③を一番右まで移動させた結果で見ていますW8-3：HSP1

①
②

①

③
Aug 29-30 2017②
250

①2つめのHSP (以下、HSP2)に切り替えて、②一番左に移動

W8-4：HSP2②
Aug 29-30 2017①251W8-4：HSP2①query側(sequence1)の555,167番目と、②DB側(sequence4)の4,852番目から…③を一番右まで移動

①
②

①

②

Aug 29-30 2017③252

①query側(sequence1)の560,027番目と、②
DB側(sequence4)の1番目の領域がHSP2。③
を一番右まで移動させた結果で見ていますW8-4：HSP2

①
②

①

③
Aug 29-30 2017②253

sequence1とsequence4のアラインメント模式図。①
sequence1の一続きの領域[549494, 560027 bp]と、②領域[4853, 5705 bp]を除く、③sequence4の全長がほぼ一致。そして、④アラインメントされなかった領域[4853, 5705 bp]には、transposaseがコードされていた(W4-9)。第8回の図1と同じ

W8-5：模式図
5’3’549,494

560,027①

sequence1

HSP2

HSP1
sequence4

11,372

5,706

4,852
3’

③
1
5’

transposaseを含んだ領域 (IS)②④

Aug 29-30 2017254

第8回原稿p189の右上①

Aug 29-30 2017①このあたりまでの話が終了。
sequence4除外後、残るはsequence1の詳細な解析。②sequence1は（ちょっとややこしいが）環状染色体、が結論②

255

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017

256

第8回原稿p187

Aug 29-30 2017

①第8回原稿のアブストラクト。特にsequence1の解釈は赤下線部分を実感できます。得られた結果をパズルゲ
ームのように組み合わせて合理的な結論を導き出す考
え方は、バクテリア以外の分野のヒトにとっても有意義かも。取捨選択してご自身の研究に活かしてください257

W9-1：seq1 vs. seq1①

query側(sequence1)は固定で、①DB側がsequence1のBlast結果を表示。②配列類似領域(HSP)は1,347個あったことを思い出そう。③ス
コアトップのHSP (HSP1)が、④seq1 vs. seq1の全長が100%一致のHSPとなるのは当たり前②④③
Aug 29-30 2017258

W9-1：seq1 vs. seq1②

①や②や③を押して、スコアの上位33個(HSP1-33)までの全体像をまとめたのが…

③①
Aug 29-30 2017259

W9-2：HSP1-33
HSP番号

スコア

HSP1
4,227,900
HSP2 (HSP3)
11,624
HSP4 (HSP5)
10,754
HSP6 (HSP7)
10,516
HSP8 (HSP9)
10,429
HSP10 (HSP11)
9,679
HSP12 (HSP13)
9,524
HSP14 (HSP15)
9,362
HSP16 (HSP17)
8,385
HSP18 (HSP19)
7,301
HSP20 (HSP21)
6,186
HSP22 (HSP23)
6,047
HSP24 (HSP26)
6,043
HSP25 (HSP27)
6,043
HSP28 (HSP29)
6,038
HSP30 (HSP31)
5,050
HSP32 (HSP33)
4,021

これです。①括弧内のHSPは、カッコ外のHSPとqueryとDB側の一致領域が入れ替わっているだけで、実質的に同じものです。この表作成は1つ1つ確認しながらの手作業

一致領域(query側)
start
end
strand
1 2,289,497 plus
2,058,465 2,065,193 plus
37,32943,187 plus
2,059,080 2,065,193 plus
2,283,820 2,289,497 plus
2,273,804 2,279,100 plus
2,059,672 2,065,193 plus
1,088,170 1,093,274 plus
2,060,310 2,065,194 plus
2,060,901 2,065,194 plus
2,061,516 2,065,193 plus
2,275,808 2,279,168 plus
999,670 1,002,969 plus
1,088,173 1,091,466 plus
2,275,808 2,279,100 plus
2,062,154 2,065,193 plus
2,062,769 2,065,197 plus

一致領域(DB側)
query側 DB側のLengthstart
end
strandの領域 領域
1 2,289,497 plus 2,289,497
2,057,850 2,064,578 plus
6,732
1
5,860 plus
5,865 ①'①
2,057,850 2,063,963 plus
6,117
5,839
11,509 plus
5,679 ②'②
717,190 711,892 minus
5,301
2,057,830 2,063,348 plus
5,525
999,673 1,004,778 plus
5,106
2,057,850 2,062,735 plus
4,887
2,057,830 2,062,120 plus
4,296
2,057,830 2,061,504 plus
3,680
1,002,969 999,608 minus
3,365
711,886 715,184 plus
3,300
711,892 715,184 plus
3,294
1,091,466 1,088,173 minus
3,295
2,057,850 2,060,889 plus
3,042
2,057,850 2,060,279 plus
2,431

①

Aug 29-30 2017260

①の数値でソートした結果

W9-2：HSP1-33
HSP番号

スコア

HSP1
4,227,900
HSP4 (HSP5)
10,754
HSP24 (HSP26)
6,043
HSP14 (HSP15)
9,362
HSP25 (HSP27)
6,043
HSP2 (HSP3)
11,624
HSP6 (HSP7)
10,516
HSP12 (HSP13)
9,524
HSP16 (HSP17)
8,385
HSP18 (HSP19)
7,301
HSP20 (HSP21)
6,186
HSP30 (HSP31)
5,050
HSP32 (HSP33)
4,021
HSP10 (HSP11)
9,679
HSP22 (HSP23)
6,047
HSP28 (HSP29)
6,038
HSP8 (HSP9)
10,429一致領域(query側)
start
end
strand
1 2,289,497 plus
37,32943,187 plus
999,670 1,002,969 plus
1,088,170 1,093,274 plus
1,088,173 1,091,466 plus
2,058,465 2,065,193 plus
2,059,080 2,065,193 plus
2,059,672 2,065,193 plus
2,060,310 2,065,194 plus
2,060,901 2,065,194 plus
2,061,516 2,065,193 plus
2,062,154 2,065,193 plus
2,062,769 2,065,197 plus
2,273,804 2,279,100 plus
2,275,808 2,279,168 plus
2,275,808 2,279,100 plus
2,283,820 2,289,497 plus

一致領域(DB側)
query側 DB側のLengthstart
end
strandの領域 領域
1 2,289,497 plus 2,289,497
1
5,860 plus
5,865 ①'①
711,886 715,184 plus
3,300
999,673 1,004,778 plus
5,106
711,892 715,184 plus
3,294
2,057,850 2,064,578 plus
6,732
2,057,850 2,063,963 plus
6,117
2,057,830 2,063,348 plus
5,525
2,057,850 2,062,735 plus
4,887
2,057,830 2,062,120 plus
4,296
2,057,830 2,061,504 plus
3,680
2,057,850 2,060,889 plus
3,042
2,057,850 2,060,279 plus
2,431
717,190 711,892 minus
5,301
1,002,969 999,608 minus
3,365
1,091,466 1,088,173 minus
3,295
5,839
11,509 plus
5,679 ②'②①

Aug 29-30 2017261

①や②の一致領域が、後の議論対象です

W9-2：HSP1-33
一致領域(query側)
start
end
strandHSP14,227,900
1 2,289,497 plus
① 37,329
HSP4 (HSP5)
10,754
43,187 plus
HSP24 (HSP26)
6,043 999,670 1,002,969 plus
HSP14 (HSP15)
9,362 1,088,170 1,093,274 plus
HSP25 (HSP27)
6,043 1,088,173 1,091,466 plus
HSP2 (HSP3)
11,624 2,058,465 2,065,193 plus
HSP6 (HSP7)
10,516 2,059,080 2,065,193 plus
HSP12 (HSP13)
9,524 2,059,672 2,065,193 plus
HSP16 (HSP17)
8,385 2,060,310 2,065,194 plus
HSP18 (HSP19)
7,301 2,060,901 2,065,194 plus
HSP20 (HSP21)
6,186 2,061,516 2,065,193 plus
HSP30 (HSP31)
5,050 2,062,154 2,065,193 plus
HSP32 (HSP33)
4,021 2,062,769 2,065,197 plus
HSP10 (HSP11)
9,679 2,273,804 2,279,100 plus
HSP22 (HSP23)
6,047 2,275,808 2,279,168 plus
HSP28 (HSP29)
6,038 2,275,808 2,279,100 plus
② 2,283,820 2,289,497 plus
HSP8 (HSP9)
10,429HSP番号

Aug 29-30 2017スコア

一致領域(DB側)
query側 DB側のLengthstart
end
strandの領域 領域
1 2,289,497 plus 2,289,497
1
5,860 plus
5,865 ①'①
711,886 715,184 plus
3,300
999,673 1,004,778 plus
5,106
711,892 715,184 plus
3,294
2,057,850 2,064,578 plus
6,732
2,057,850 2,063,963 plus
6,117
2,057,830 2,063,348 plus
5,525
2,057,850 2,062,735 plus
4,887
2,057,830 2,062,120 plus
4,296
2,057,830 2,061,504 plus
3,680
2,057,850 2,060,889 plus
3,042
2,057,850 2,060,279 plus
2,431
717,190 711,892 minus
5,301
1,002,969 999,608 minus
3,365
1,091,466 1,088,173 minus
3,295
5,839
11,509 plus
5,679 ②'②

262

W9-2：HSP1-33
HSP番号

スコア

HSP1
4,227,900
HSP4 (HSP5)
10,754
HSP24 (HSP26)
6,043
HSP14 (HSP15)
9,362
HSP25 (HSP27)
6,043
HSP2 (HSP3)
11,624
HSP6 (HSP7)
10,516
HSP12 (HSP13)
9,524
HSP16 (HSP17)
8,385
HSP18 (HSP19)
7,301
HSP20 (HSP21)
6,186
HSP30 (HSP31)
5,050
HSP32 (HSP33)
4,021
HSP10 (HSP11)
9,679
HSP22 (HSP23)
6,047
HSP28 (HSP29)
6,038
HSP8 (HSP9)
10,429Aug 29-30 2017

右側の①’が領域[37329, 43187 bp]、①が[1, 5860 bp]、②’が [2283820, 2289497 bp]、②が[5839, 11509 bp]に相当。ここの「①’, ①, ②’, ②」が後の議論対象

一致領域(query側)
start
end
strand
1 2,289,497 plus
37,32943,187 plus
999,670 1,002,969 plus
1,088,170 1,093,274 plus
1,088,173 1,091,466 plus
2,058,465 2,065,193 plus
2,059,080 2,065,193 plus
2,059,672 2,065,193 plus
2,060,310 2,065,194 plus
2,060,901 2,065,194 plus
2,061,516 2,065,193 plus
2,062,154 2,065,193 plus
2,062,769 2,065,197 plus
2,273,804 2,279,100 plus
2,275,808 2,279,168 plus
2,275,808 2,279,100 plus
2,283,820 2,289,497 plus

一致領域(DB側)
query側 DB側のLengthstart
end
strandの領域 領域
1 2,289,497 plus 2,289,497
1
5,860 plus
5,865 ①'①
711,886 715,184 plus
3,300
999,673 1,004,778 plus
5,106
711,892 715,184 plus
3,294
2,057,850 2,064,578 plus
6,732
2,057,850 2,063,963 plus
6,117
2,057,830 2,063,348 plus
5,525
2,057,850 2,062,735 plus
4,887
2,057,830 2,062,120 plus
4,296
2,057,830 2,061,504 plus
3,680
2,057,850 2,060,889 plus
3,042
2,057,850 2,060,279 plus
2,431
717,190 711,892 minus
5,301
1,002,969 999,608 minus
3,365
1,091,466 1,088,173 minus
3,295
5,839
11,509 plus
5,679 ②'②

263

W9-3：アノテーション
一致領域(query側)
start
end
strandHSP14,227,900
1 2,289,497 plus
HSP4 (HSP5)
10,75437,32943,187 plus
HSP24 (HSP26)
6,043 999,670 1,002,969 plus
HSP14 (HSP15)
9,362③ 1,088,170 1,093,274 plus
HSP25 (HSP27)
6,043 1,088,173 1,091,466 plus
HSP2 (HSP3)
11,624 2,058,465 2,065,193 plus
HSP6 (HSP7)
10,516 2,059,080 2,065,193 plus
HSP12 (HSP13)
9,524 2,059,672 2,065,193 plus
HSP16 (HSP17)
8,385 2,060,310 2,065,194 plus
HSP18 (HSP19)
7,301 2,060,901 2,065,194 plus
HSP20 (HSP21)
6,186 2,061,516 2,065,193 plus
HSP30 (HSP31)
5,050 2,062,154 2,065,193 plus
HSP32 (HSP33)
4,021 2,062,769 2,065,197 plus
HSP10 (HSP11)
9,679 2,273,804 2,279,100 plus
HSP22 (HSP23)
6,047③ 2,275,808 2,279,168 plus
HSP28 (HSP29)
6,038 2,275,808 2,279,100 plus
HSP8 (HSP9)
10,429 2,283,820 2,289,497 plus
HSP番号

Aug 29-30 2017

スコア③のあたりはribosomal RNAがほとんど。DFASTで得られたアノテーション情報を追加しています

一致領域(DB側)
query側 DB側のLengthstart
end
strandの領域 領域
1 2,289,497 plus 2,289,497
1
5,860 plus
5,865 ①'①
711,886 715,184 plus
3,300 ribosomal RNA
999,673 1,004,778 plus
5,106 ribosomal RNA
711,892 715,184 plus
3,294 ribosomal RNA
2,057,850 2,064,578 plus
6,732
2,057,850 2,063,963 plus
6,117
2,057,830 2,063,348 plus
5,525
2,057,850 2,062,735 plus
4,887
2,057,830 2,062,120 plus
4,296
2,057,830 2,061,504 plus
3,680
2,057,850 2,060,889 plus
3,042
2,057,850 2,060,279 plus
2,431
717,190 711,892 minus
5,301 ribosomal RNA
1,002,969 999,608 minus
3,365 ribosomal RNA
1,091,466 1,088,173 minus
3,295 ribosomal RNA
5,839
11,509 plus
5,679 ②'②

264

W9-3：アノテーション④のあたりは、「adhesion exoprotein、
mucus-binding protein、hypothetical protein」となっており、はっきり言ってよくわからない

一致領域(query側)
start
end
strandHSP14,227,900
1 2,289,497 plus
HSP4 (HSP5)
10,75437,32943,187 plus
HSP24 (HSP26)
6,043 999,670 1,002,969 plus
HSP14 (HSP15)
9,362 1,088,170 1,093,274 plus
HSP25 (HSP27)
6,043 1,088,173 1,091,466 plus
HSP2 (HSP3)
11,624 2,058,465 2,065,193 plus
HSP6 (HSP7)
10,516 2,059,080 2,065,193 plus
HSP12 (HSP13)
9,524 2,059,672 2,065,193 plus
HSP16 (HSP17)
8,385 2,060,310 2,065,194 plus④
HSP18 (HSP19)
7,301 2,060,901 2,065,194 plus
HSP20 (HSP21)
6,186 2,061,516 2,065,193 plus
HSP30 (HSP31)
5,050 2,062,154 2,065,193 plus
HSP32 (HSP33)
4,021 2,062,769 2,065,197 plus
HSP10 (HSP11)
9,679 2,273,804 2,279,100 plus
HSP22 (HSP23)
6,047 2,275,808 2,279,168 plus
HSP28 (HSP29)
6,038 2,275,808 2,279,100 plus
HSP8 (HSP9)
10,429 2,283,820 2,289,497 plus

一致領域(DB側)
query側 DB側のLengthstart
end
strandの領域 領域
1 2,289,497 plus 2,289,497
1
5,860 plus
5,865 ①'①
711,886 715,184 plus
3,300 ribosomal RNA
999,673 1,004,778 plus
5,106 ribosomal RNA
711,892 715,184 plus
3,294 ribosomal RNA
2,057,850 2,064,578 plus
6,732
2,057,850 2,063,963 plus
6,117
2,057,830 2,063,348 plus
5,525
2,057,850 2,062,735 plus
4,887
2,057,830 2,062,120 plus
4,296
2,057,830 2,061,504 plus
3,680
2,057,850 2,060,889 plus
3,042
2,057,850 2,060,279 plus
2,431
717,190 711,892 minus
5,301 ribosomal RNA
1,002,969 999,608 minus
3,365 ribosomal RNA
1,091,466 1,088,173 minus
3,295 ribosomal RNA
5,839
11,509 plus
5,679 ②'②

HSP番号

Aug 29-30 2017

スコア

265

W10-1：sequence1
2,283,820
2,289,497(a)sequence1の末端付近の重複は、HSP8の一致領域（②と②’）に相当。これにHSP4の一致領域（①と①’）を含めた模式図。第8回の図2aと同じ。後の検証により、①の領域
は染色体上に存在しないことが確認された

sequence1末端部分②’
sequence1
先頭部分

Aug 29-30 2017

5,839

1

①

11,509②

37,32943,187①’

266

W10-2：アノテーション

HSP4の一致領域①[1, 5860 bp]のDFAST
アノテーション。HSP8の一致領域②[5839,
11509 bp]のアノテーションの一部。マニア
ックすぎるので、このあたりは飛ばします。第8回ウェブ資料に詳細情報あり

①

①

②

②

Aug 29-30 2017267

W10-3：染色体構造(b) sequence1の両末端である①と②’の領域をトリム。得られた領域[5839, 2283819 bp]の両末端を結合した環状コンティグが実際の染色体構造であると予想した

2,283,820
2,289,497(a)sequence1末端部分②’
sequence1
先頭部分

5,839

1

①

11,509②

37,32943,187①’(b)

プロファージ領域

Aug 29-30 2017268

①のあたりまでが終了

第8回原稿p190右上①

Aug 29-30 2017269

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017

270

①のあたりと…

第8回原稿p190右中①

Aug 29-30 2017271

①のあたりを行います

第8回原稿p191左上①

Aug 29-30 2017272

①アセンブリ結果ファイル(LH_hgap.fa)を入力として、重複除去を行う一連のコマンド。
②出力ファイルはLH_draft.fa。③ファイル
サイズの減少度合い的に妥当。ここでは示さないがlessで配列の最初と最後を狙い通りに抽出できているか確認しておこう

W11-2：ドラフト配列作成①②③

Aug 29-30 2017273

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017274

この後の展開（第8回）


W12：bwaでマッピング①自分は変異解析をやるつもりがなく、VCFフ
ァイルの理解はこれまで拒否していました。今回ドラフト配列の修正という局面で用いることになり、嫌々ながら必要最小限の勉強をしました （Bio-Linuxにプレインストールされている）bwaを利用し、PacBioデータから得られた

ドラフト配列をリファレンスとして、Illumina MiSeqデータをマップ。
 やり方のみを示し、この結果は使わない。後述のDDBJ Pipelineの結果のほうがVCFファイルまで出力してくれて便利だから。


W13：CyberduckのインストールとMiSeqデータのアップロード
 CyberduckはFTPソフト。DDBJ Pipelineへのデータアップロード用として利用。

W14：DDBJ Pipeline上でbwaを実行 W15：実行結果の概観とダウンロード
 W16～W17：マッピング結果ファイル（SAM/BAM形式）の説明
 W18：VCFファイルの説明 PacBioデータから得られたリファレンス配列と違っている部分を示したもの。リード数およびクオリティスコアの点でIllumina MiSeqデータのほうが優位。VCFファイルの結果をもとに手作業でドラフト配列を修正するための基礎知識や基本的な考え方を ①身につけることを目的として、VCFファイルの読み取り方を解説。VCF形式について
は、（平成27-28年度のNGSハンズオン講習会の）変異解析」でも紹介。

Aug 29-30 2017275

この後の展開（第8回）
①Tabletの解説も兼ねています。W12からW27の範囲は、PacBioデータしかない場合は行いません。MiSeqデータもあるので援用しています

W19：Tabletのインストールと利用 Tabletは、リファレンス配列へのマッピング結果をローカル環境で可視化するViewer。

VCFファイルでみた変異箇所は、ドラフト配列の修正候補箇所に相当する。MiSeqデータのマッピング結果のほうを採用するかどうかを、Tabletで修正候補箇所を実際に眺めて判断する。

W20：plasmid2上の変異箇所を確認
 W21：plasmid1上の変異箇所を確認①
 W22：chromosome上の変異箇所を確認
 W23：マップされたリードの詳細情報を知るTips
 W24：ドラフト配列修正方針のまとめとVCFファイルのおさらい W26：変異の反映…が意外と文字列置換でうまくいかないことを学ぶ W27：テキストエディタ(EmEditor)を用いてベタで修正 文字列検索と該当塩基位置情報を組み合わせて、ドラフト配列を修正

Aug 29-30 2017276

この後の展開（第9回）
W3：Ori-Finderで転写開始点を同定
 W4： 配列の“回転”

連載第9回以降は、Illumina MiSeqデータの有無に関係なく行う。①は、環状バクテリア
ゲノムの場合にしばしば慣例として行われる。GC skewで眺めたときに、時計の0時あたりで分布が切り替わるようにするのが主目的

 dnaA遺伝子が配列の先頭となるように“回転”

W5：INSDCフラットファイルの説明①

 entry, feature, qualifierからなる3つの階層構造の話など

W6：locus_tagの説明
 W7とW8：DFASTアノテーションの実行と結果の確認、Tips
 W9とW10：DFAST結果画面からのDDBJへの登録の概要
 W11とW12：描画ソフトDNAPlotterのインストールと利用 W13：grepを駆使してgenbank(.gbk)形式を理解
 W14：awkコマンドの利用（最初の配列のみ抽出）
 W15とW16：DNAPlotterを再度利用（プラスミド配も描画し） W17：描画結果を原著論文の図と比較
 W18：ovaファイルの再インストール関連情報
 W19：配列情報つきgffファイルを作成してDNAPlotterで描画Aug 29-30 2017277

この後の展開（第10回）

連載第10回は、de novoアセンブリ結果フ
ァイルをDFASTアノテーション結果とともに、DDBJに登録する一連の作業を解説

W2：DDBJの登録用アカウント作成
 W3：DDBJにログイン W4：NSSSとMSSの説明など（今回はMSSでの登録）
 W5とW６：BioProjectの登録。Locus tag prefixのおさらいなど W7とW8：BioSampleの登録。サンプル属性のテンプレートファイルの解説
 W9：アクセッション番号取得とDDBJによる編集結果の確認
 W10：MSS (Mass Submission System)での申込
 W11：DFAST上でDDBJ登録用ファイルを作成
 W12：DDBJとのメールのやりとり
 W13：公開
Aug 29-30 2017278

Contents（主に第8回）
W4：ゲノムアノテーション（DFAST）
 W5：dotterの実行（sequence4）
 W6：Blastの実行（DB配列はLH_hgap.fa、query配列はsequence1）
 W7：BlastViewerでBlast実行結果を眺める（W7-4まで）
 W7：lessとgrepでも確認（W7-6まで）
 W8：Blast実行結果のsequence1 vs. sequence4を眺める
 W9～W10：sequence1を詳細に調べる（Blastとアノテーション結果を併用）
 W11：乳酸菌ゲノム概要配列の作成
 この後の展開（第8回のW12以降と第9-10回の概要）
 W3：シェルスクリプトAug 29-30 2017279

W3-0：ファイル分割1①第8回W11-2(スライド273)では、3つの配列からなるmulti-FASTAファイル(LH_draft.fa)を作成した。よって、②配列数は3、③行数は6

①

②

③

Aug 29-30 2017280

W3-0：ファイル分割2

第7回W10-3からW10-4（スライド56-58）と同じようなファイル分割方法だと、①こんな感じになる。この例のように配列数が10個以下程
度なら、私は手作業でやります。しかし、配列
数が100とか1000個程度などの局面では、シェルスクリプトがおそらくよく使われます①

Aug 29-30 2017281

W3-1：シェルスクリプト①

Aug 29-30 2017

先程とは異なり、①contig[0-9].fa作成のためのスクリプトを統一的に記述している。②や③の記述法は一見無駄に見えるが、このような書き方でも目的を達成できるところがポイント

②

③

282

①で作成した、②contig[0-9].faと、先に作成したmongee[0-9].faのファイルサイズは確かに同じ

W3-1：シェルスクリプト

①②③

Aug 29-30 2017283

W3-1：シェルスクリプト①ファイルサイズだけでなく、中見も完全に同じであることがわかります①

Aug 29-30 2017284

目的を再確認

ここでは、配列数が100とか1000個程度からなるmultiFASTAファイルの分割を想定して、シェルスクリプトで記述するやり方を伝授しようとしています。①のようなコマンドを100とか1000個書かなくて済むので覚える価値があります①

Aug 29-30 2017285

W3-2：参考資料

参考資料は、例えば①平成27年度NGSハンズオン講
習会の、②7/24の講義資料（のスライド58以降）です①②
Aug 29-30 2017286

W3-3：戦略を練る

まず、①の部分は2で固定。理由は、出力ファイルがsingle-FASTA形式だから①

Aug 29-30 2017287

W3-3：戦略を練る②

Aug 29-30 2017

変わっているのは①と②の部分。①13までのループを回して、その2倍の値
として②を表現できるだろうと妄想する①

288

W3-4：基本形①シェルスクリプトの基本形ファイルJST-NBDC_1.sh
を、②wgetし、③moreで確認②③①

Aug 29-30 2017289

W3-5：解説①1から②3までの1刻みのループを回し、③iという変数で取り扱う。別にiでなくてもよいが、変数名として、i, j, kといった順番で使うヒトも一定数存在する①③

Aug 29-30 2017②290

W3-5：解説①変数iは、②$iとして取り扱う。つまり$を追加する。なぜ?という類のものではなく、お約束です①②

Aug 29-30 2017②291

W3-5：解説①

Aug 29-30 2017①「echo “…”」で囲んでいますが、これはLinuxコマンドのechoであり、“”内部の実際に実行したいコマンドがどんな感じになっているかを実行前に確認するためにつけています①

292

W3-6：echoで確認①

Aug 29-30 2017

①shコマンドでJST-NBDC_1.shを実行。
②echoで囲った中身に相当する、赤下線部分が表示されていることがわかる②

293

W3-7：head部分

実際に全体を実行してエラーに遭遇してもよいが一応説明。headコマンド部分の①は最初のy行という数値を指定するところだが、1*2とか3*2という掛け算の*が含まれている。このような指定法はダメ（だということをエラーに遭遇して学習する）①

Aug 29-30 2017294

W3-7：head部分①や②を実際に実行してみても、赤下線部分が原因
でエラーが出ます。エラーが出ないようにするには、③のように3*2の結果である6を与える必要があります①②③

Aug 29-30 2017295

リターンキーを押した後の状態

W3-7：head部分

Aug 29-30 2017296

W3-8：発展形①シェルスクリプトの発展形ファイルJST-NBDC_2.shを、②wgetし、③moreで確認②③①

Aug 29-30 2017297

W3-8：発展形

Aug 29-30 2017

JST-NBDC_1.shとの違いは、 赤下線部分のみ。jという別の変数を用いて$i*2を表現したいだけだが、シェルスクリプトの掟に従うと、このようになります298

W3-8：発展形①は$i*2の*に相当する部分ですが、「なんでもよい、という意味でのワイルドカードの*」と「掛け算の意味での*」を区別する必要があります。後者の意味として用いた
い場合に、②¥(入力できませんが、バックスラッシュです)を添えます。本当にこの説明でいいのかは怪しいです①②

Aug 29-30 2017299

W3-8：発展形

掟の続き。①や②のところにス
ペースを入れてはいけません①②

Aug 29-30 2017300

W3-8：発展形①

掟の続き。①のところにはスペースを入れないといけません①①

Aug 29-30 2017301

W3-8：発展形①①の記号は、②Shiftキーを押しな
がら、③@のキーを押すと出ます

①

③

②

Aug 29-30 2017302

W3-9：echoで確認①shコマンドでJST-NBDC_2.shを実行。②意図通りになっていることがわかる①②

Aug 29-30 2017303

W3-10：発展形2①シェルスクリプトの発展形2のファイル
JST-NBDC_3.shを、②wgetし、③moreで確認

②

③

①

Aug 29-30 2017304

JST-NBDC_2.shとの違いは、①と②の行の部分。①のecho行
頭に#を入れてコメントアウト（実行されないように）している。
削除するのと同じだが、エラーなど問題が起こったときの対処用などの目的で#をつけたまま残しておくことはよくやる。②は①の赤下線部分と同じもの。ここが実際に実行される部分

W3-10：発展形2②

Aug 29-30 2017①305

W3-11：実行①shコマンドでJST-NBDC_3.shを実行。echoを
コメントアウトしているので、意図通り何も表示されない。また、②イメージ通りの出力ファイル
test[0-9].faが作成されていることがわかる①②②

Aug 29-30 2017306

W3-12：確認①「ls –l」で詳細情報を表示。ファイルサイズの観点か
らは、contig[0-9].faとtest[0-9].faは全く同じ。②念のためdiffコマンドでも確認。何も表示されていないので中身も同一であることを確認したことになる。③違うもの同士だと、違いのある部分がババッと表示されます①②③

Aug 29-30 2017307

W3-12：確認

Aug 29-30 2017

リターンキーを押して、ババッと表示された後の状態です。diffコマンドの挙動については、自分で納得できるようなファイルを用意して確認するのが一番ですので、ぜひやってみてください308

W3-13：発展形3「expr シェル 遅い」などでググればわかりますが、
JST-NBDC_3.sh内で使われているexprは遅いらしいです。赤下線のように書くといいらしい。①高速版のJST-NBDC_4.shを②wgetし、③more。後は省略②③①

Aug 29-30 2017309