Help us understand the problem. What is going on with this article?

SPAdesを使ってゲノムアセンブリしてみる

初めに

本ページは岐阜大学3年ゲノム微生物学の講義の中で解説される de novo assemblyを実践するための参考サイトです。半期でゲノムアセンブリをこなすのはなかなか難しいので、今回は簡単な方法を紹介しています。基本的には、コピペと少しのタイピングで済むように作っていますが、チャレンジできる人はまずは全て写経すると仕組みが仕組みが良く分かると思います。他のページとの違いはど素人が書いているという点です。とにかく分かりやすいように書いているつもりです。本稿ではMacを利用します。免責事項として本稿の内容がすべて正しいとは限りません。その点理解して読み進めてください。

生物学の進歩

Wetを行う人間であっても、鳥瞰図的に物事が見られるBioinfomaticsは今後最低限の素養として必要になってくる学問であると考えています。どういうことかといえば、例えば研究のスピードがあげられるというのが分かりやすいかと思います。未知の生物のゲノム配列が分かるだけでどんな遺伝子を持っているかを大まかに知れますし、イントロン領域からも出ているRNAの配列を知ることもできます。Bioinfomaticsはあくまで計算科学なので、その点を補うのがWetの人の仕事です。結果を生かしてWetに行き、最終的に実験によって確かめるというのが重要なのかなと筆者は考えます。何れにせよ、やって損はない勉強です。一緒に頑張りましょう。

コマンド操作

De novo assemblyは基本的にCUI (コマンド操作, Character-based User Interface) です。Macではターミナルというアプリケーションを使うことで操作できます。ターミナルはLaunch Padから、「その他」に行くと見つけられます。

コマンド操作で重要なのは名前があっているか、コマンドが正しいかになります。例えば、下手に半角スペースがあるだけで動かないということが多々あります。ここに書いてあるコマンドは動作したコマンドになります。完全に同じコマンドになるようにしてください。ターミナルはUNIXと呼ばれるOSをベースにして動作しています。詳しく学びたい人は以下のサイトがオススメです。

https://codezine.jp/unixdic/i/UNIX基礎講座

コマンドの打ち方
$の後ろに打ちます。スペースにも意味があります。
ターミナルを開くと以下のような文字が表示されます。

MacBook-Air:~ RyoNiwa$

例えば、cdというコマンドを打つ際は以下のように打ちます。

MacBook-Air:~ RyoNiwa$ cd

もしもうまくできなかった時は

15分ググってください。検索ワードを変えて、英語に変えて、エラーをコピペして...いろんな手段でググってください。それでも解決できなかった時はv8024064@edu.gifu-u.ac.jpにメールください。助けられる範囲で返信します。

環境準備編

第1章ではパソコンの準備をします。パソコン内に最低限必要なものをインストールするという作業を行います。
今回はHomebrewと呼ばれるmacOS用のパッケージマネージャーをインストールして、単純に作業が進められるようにします。
本稿ではUNIXベースでデータを処理を行います。できる限りインストールを少なくすることで、簡単な方法を選んでいるつもりです。
本来であれば、使用するデータに応じて使用するソフトウェアを変えることはデータをよりよく準備する意味で肝要です。
もし、今回の処理を通して少しでも興味を持ってもらえた場合は継続的に学び、少し難しい課題にもトライして欲しいと考えます。

Homebrewのインストール

HomebrewはAppleのmacOSに含まれていないソフトウェアなどをダウンロードする際に利用するパッケージマネージャーです。
Homebrewは個別のディレクトリにパッケージをインストールし、それらへのシンボリックリンクを/usr/localに作ります。
難しく考えずに、MacOSで足りない部分を補強してくれるパッケージだと考えれば大丈夫です。以下のコマンドをターミナルにコピペするとインストールされます。

$ /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

例としてwgetコマンドをインストール

Homebrewを使って、早速wgetコマンドを使えるようにしてみます。wgetコマンドはLinuxで通常用いられるコマンドでターミナルでは使用できません。そこで、Homebrew内にwgetコマンドを入れて使えるようにしてみます。

$ brew install wget

そうすると、頭の良さそうな動きをターミナルがします。
終わったら、wgetがインストールできたか確かめます。

$ wget -h

以下の文字が出てきたら成功。

RyonoiMac:~ ryoniwa$ wget -h
GNU Wget 1.20.3, 非対話的ネットワーク転送ソフト
使い方: wget [オプション]... [URL]...

長いオプションで不可欠な引数は短いオプションでも不可欠です。

できればMac、Windows関係なくできる方法を提案したいと考えていますが、まずはHomebrewでのやり方を進めていきます。すみません。

Javaの準備

MacはデフォルトでJavaが入っていません。今回はJavaを使うソフトウェアを使用しますので、Javaをあらかじめインストールしておきます。
古いJavaが入っている人も同じ方法でおそらくアップデートできますので、やっていきます。まずはJavaがあるか見ていきましょう。

$ java -version

ここで、コマンドがないと出たら、やはり入っていません。諦めて、https://support.apple.com/ja-jp/HT204036 を参照にJavaをインストールしてください。グラフィックできるので、ここからの作業は省略します。説明に沿って、ぽちぽちクリックするだけです。

データ準備

本章では、解析データの準備を行います。SRAと呼ばれるアーカイブからデータを落として、計算できる状態にします。
今回は乳酸菌のゲノムを解析することにしました。作業が煩雑ですので、じっくりやっていきましょう。

SRAからデータをダウンロードする

SRAはSequence Read Archiveの略で、その名の通りシーケンスデータをレポジトリしているデータベースです。基本的にはDDBJ SRA、 NCBI SRA、 EMBL-EBI ENAの3つのデータベースがあり、日本のデータベースからダウンロードするのが一番早いです。例えば、「ゲノムアセンブリした」というような論文にはAccession Numberと呼ばれる番号が記載されているので、これをデータベースで検索します。

ダウンロードするデータ

今回は適当に検索して、一番最初にヒットしたLactobacillus delbrueckiiのゲノムデータをアセンブルすることにしました。ただし、基本的にはアルファベットを変えてデータを落とせばいいだけなので、好きな乳酸菌やその他菌株があれば論文で探して落としてもらって構いません。

Complete Genome Sequencing of Lactobacillus plantarum ZLP001, a Potential Probiotic That Enhances Intestinal Epithelial Barrier Function and Defense Against Pathogens in Pigs

早速データをダウンロードしていきます。せっかくなので、日本サーバのDDBJ SRAで検索していきます。
1. GoogleでDDBJ SRAを検索して、トップでヒットしたサイトを開く。
2. 真ん中の方にある検索をクリック。
3. 一番上のAccessionに論文中に出てくるAccession Numberを入れる。(先の論文であればSRP102895)
4. 右のほうにあるExperimentのSRAをクリック、その後SRR5407012をクリック。
5. SRR5407012.sraをクリック。ダウンロードが始まるので少し待つ。

ここで、ダウンロードのフォルダを開くとSRR5407012.sraがあるはずです。ただ、ここで問題なのがSRR5407012.sraが特殊な形式であり、開くことができないという点にあります。そこで、通常のNGSのアウトプット形式であるfastqファイルに書き換えたいわけです。そのために、sra-toolkitという先のNCBI SRAが公開しているソフトをダウンロードして、パスを通す必要があります。

SRA Toolkitの準備

  1. まずは以下サイトにアクセスする。

https://www.ncbi.nlm.nih.gov/sra/docs/toolkitsoft/

  1. 続いて、MacOS 64bitのSRA toolkitをクリックしてダウンロードする。

ただし、通常のアプリケーションと異なり、SRA toolkitはグラフィカルユーザーインターフェースとかではないので、いわゆるパスを通すという作業が必要になります。パスを通すというのはプログラムの名前だけでプログラムをサクッと実行できるようにすることです。ざっくり言えば、ポケモンでピカチュウにカミナリという攻撃を覚えさせるのと同じです。

さて、SRA toolkitを早速使えるようにしていきます。その前に一旦SRA toolkitの入ったフォルダを解凍しデスクトップに移します。ここからは作業を丁寧に行っていきましょう。まずはTerminal上でSRA toolkitを開いてみます。ターミナルで表示される文字はユーザ名により異なる。作業はMacBook Airで行っているが、本文はiMacで書いているので、ユーザはiMacとなっているが気にしないでください。使うコマンドはたった3つです。残りのコマンドはコピペすれば大丈夫です。

cd: Change Directory
ディレクトリを変更。

$ cd [ディレクトリ名]

pwd: Print Working Directory
今どこのディレクトリにいるのかを教えてくれる。

$ pwd

ls: List segments
今いるディレクトリの中にあるファイルやディレクトリを教えてくれる。

$ ls

ちなみにTabボタンをいかにうまく使うかがコマンドを早く書くコツだと考えてください。
そのディレクトリ内で固有のワードになった時点でTabを押せば残りはパソコンが勝手に書いてくれます。

SRA Toolkitの場所の確認

$ cd Desktop/
$ ls

そうすると、sratoolkit.2.9.6-1-mac64というディレクトリが存在すると考えられます。
バージョン名が違うかもしれないが、おおよそこの名前であるだろうと思います。
続いて、中身を確認してみます。

iMac:Desktop ryoniwa$ cd sratoolkit.2.9.6-1-mac64/
iMac:sratoolkit.2.9.6-1-mac64 ryoniwa$ ls

CHANGES         README.md       schema
README-blastn       bin
README-vdb-config   example

iMac:sratoolkit.2.9.6-1-mac64 ryoniwa$ cd bin
iMac:bin ryoniwa$ pwd
/Users/ユーザー名/sratoolkit.2.9.6-1-mac64/bin

となる。このうち、binが重要なファイルです。このbinにパスを通すとSRA Toolkitが使えるようになります。そこで、早速パスを通していきます。以下を見て挑戦してください。

http://motomichi-works.hatenablog.com/entry/2016/03/02/103207

今回.bash_profileに記入するのは以下の通りです。

$ export PATH=$PATH:/Users/ryoniwa/Desktop/sratoolkit.2.9.6-1-mac64/bin/

うまくできているかを確認してみます。試しにSRA toolkitの一つの機能であるfastq-dumpのヘルプを見てみます。

$ fastq-dump -h

色々出てきたら成功です。

ファイルの処理 (形式変更)

続いて、ファイルを処理していきます。
先ほどダウンロードしたSRAファイルを全てデスクトップに移しておきます。
今のディレクトリがDesktopとなるようにします。どのディレクトリにいても、cd単発でコマンドを動かすと一番最初のディレクトリに戻れます。

$ cd
$ cd Desktop/

もしくはcd ..でも一つ前のディレクトリに戻れます。今のディレクトリはRyonoiMac:[今のディレクトリ]で確認可能です。
今回はペアエンドリードのデータなので、別々に出力させます。(--split-filesオプションを使用)

iMac:~ ryoniwa$ cd
iMac:~ ryoniwa$ cd Desktop/
iMac:Desktop ryoniwa$ fastq-dump --split-files SRR5407012.sra

アウトプットされるファイルは全てデスクトップで確認できます。(1番がフォワード、2番がリバース) ただし、このまま開くのはお勧めできません。普通のパソコンだとしばらく動かせなくなるからです。例えば、利用するファイルだと1 GB以上です。どうしても中身が見たければheadコマンドを使いましょう。上から10行を確認できます。

iMac:Desktop ryoniwa$ head SRR5407012_1.fastq 
@SRR5407012.1 1 length=150
CTTTTTGATCGTATTATGGTAACGAAAAGTTCCATGCGCCATGGCCGCAATGCAACTTTTCACCAACTAACAATTAGTTGTTTACTTTTACGACTTGCTTGCCGTTGTAGTATCCACAACTTGGGCAAACACGGTGTGATACACGTAATT
+SRR5407012.1 1 length=150
AAFFFJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJAFJJJJJJJFJJJJJAJJJJJFJJFJJJJFAJFJJJJJJJJFFJJJAJJF-FFJJJJFFJF<FFFJJJJF<FJF
@SRR5407012.2 2 length=150
AACGGCTCGGATTGCGGTACAGGCAGCGCTAAGTGGTCACTTAGTTTTAAGTACTATTCATGCGCGCGGTGTCTACGGGATTATTCCCCGGTTGACGCAATTGGGAGTGGAAGCCGCTGTTTTGGAGCAGGCTTTAACTGCAATGTCCTA
+SRR5407012.2 2 length=150
AAFFFJFJJJFJJFJFJFJFJFFFJJJJFJJJ<JFFJJJJAFJJFAJJJJJJJAAFFJJJAJJJJFFJJJJJJJJJJJJJJJJJJJJJFFFJJJJFFJFJFFFJFFF<JFFJFJAFAAF)7AFJJJJJJJJFFJJJJJJJJJJFJFJFJJ
@SRR5407012.3 3 length=150
AACATAAGATCCCCTCCTCAATTGTCACTATAACGCCCTGACAACATGATAGCCACCAAGATGCACAGAGGCCCGAAATCTTTTCTATGCCTATCGATTCCAGCCTAATACTGAACAAGTCTGGAATGGTTGCCCTTGTCTATCATTCCA

その後、新しいフォルダを作成してください。(当方はTrainingという名前にしました) 全てのダウンロードしたファイルをフォルダ内に入れてください。これで全ての準備は整いました。いよいよ、データ解析に移りましょう。

クオリティートリミング

クオリティートリミングについては当方Githubをご覧ください。
https://github.com/Ryo-NIWA-GIFU-Oubi/RyoNIWA-Lecture/wiki/第3章-クオリティートリミング編--座学

クオリティートリミングをやってみる

早速クオリティートリミングを始めてみましょう。使うソフトウェアは先ほど利用したTrimmomaticです。TrimmomaticはJavaベースで動くソフトウェアになります。Trimmomaticに関してはコマンドを打つというよりもスクリプトを動かす形になります。Javaベースで動きます。先に述べたとおり、Javaを利用したことない人はまずJavaが入っているかを確認してください。

$ java -version

きちんとJavaがインストールされていれば、バージョンに関しての説明が出てきます。バージョンに関して出てこず、むしろコマンド分かりませんという説明が出てこれば、それはインストールからということになります。

-bash: java: command not found

javaのインストールは簡単です。Javaのサイトにアクセスして、適当なdmgファイルをダウンロードしてください。これを開いて、指示に従ってインストールするだけです。インストール後は念の為、先のコマンドでバージョンを確認してください。

ここからはインストールできているものとして、話を進めます。続いて、Trimmomaticをダウンロードします。以下サイトにアクセスしてください。

http://www.usadellab.org/cms/?page=trimmomatic

サイト内のBinaryというところをクリックするとTrimmomaticのことが書かれたjavaファイルをダウンロードできます。これをデスクトップ上に作ったTrainingフォルダに入れましょう。ここで、Trainingフォルダを見ると、当然なのですが、先に入れたFastqファイルとjavaのファイルが入っています。同じディレクトリに入れることで、1. 分かりやすくする、2. コードを書くのが楽になると思います。慣れてきたら、自由にディレクトリを作ってやるといいですが、最初はルールを作ってルール通りに進めると失敗が少ないです。

さて、いよいよTrimmomaticを動かしてクオリティートリミングしてみましょう。順番としては

  1. ディレクトリを移動して、使用しているディレクトリをTrainingにする (cd Training/)
  2. lsで中身を確認する (ls)
  3. Fastqとjavaファイルが入っていたら、javaファイルを動かしてFastqを処理する
  4. 処理後のファイルがあるかどうか確認する
iMac:~ RyoNiwa$ cd Desktop/
iMac:Desktop RyoNiwa$ cd TRAINING/
iMac:Desktop RyoNiwa$ ls

その後、以下のように書きます。trimmomatic-0.39.jarは2019年12月5日現在で最新バージョンですが、万が一更新されていた場合は数字を変えてください。それぞれの設定、何が起こるかは

https://bi.biopapyrus.jp/rnaseq/qc/trimmomatic.html

を見てください。非常に詳しく書いていただいています。

$ java -jar trimmomatic-0.39.jar \
PE \
-threads 2 \
-phred33 \
-trimlog log.txt \
ファイル名_1.fastq \
ファイル名_2.fastq \
paired_output_1.fq  \
unpaired_output_1.fq \
paired_output_2.fq \
unpaired_output_2.fq \
LEADING:20 \
TRAILING:20 \
SLIDINGWINDOW:4:25 \
MINLEN:100

基本的に paired_output_1.fq と paired_output_2.fq を利用して、次のクオリティコントロールあるいはマッピングなどの解析に用いれば大丈夫です。今回のトリミング方法については何も考えずに適当に作りました。実際には論文を参考にして考えるのがいいと思います。また、rawデータの場合、アダプターが付与されています。その際はILLUMINACLIPのオプションを追加して、アダプタートリムとクオリティートリムを同時に進行すると良いです。実行すると、どれくらいのリードがサバイブしたかなどが表示されます。「paired_output_1.fq unpaired_output_1.fq paired_output_2.fq unpaired_output_2.fq」の4つのfq (fastqと同じです) がTrainingフォルダに出力されたか確認してください。

続いて、きちんと処理されたかをFastQCにより確認していきます。

クオリティートリミングできているかを目視する

本来はクオリティーチェック、クオリティートリミング、再クオリティーチェックをするのが慣例となっています。今回はあくまでトライアルなので、最初のクオリティーチェックは省いていますが、自分でデータを取得した際は当然やるべきだと考えています。

クオリティーチェックを行う際、最強のソフトウェアが存在します。それがFastQCです。以下このサイトからの引用です。

FastQC はデータのQCに最もよく用いられている(と言っても過言でない)オープンソースのソフトウェアです。
シーケンス実験の品質をチェックするために便利な幾つかのプロットを生成してくれます。ちなみにjavaで書かれています。
出力されるプロットとしては、

サイクル(リードの塩基)あたりのシーケンスクオリティ
クオリティスコアの分布
同一シーケンスリードの重複率
などが挙げられます。 こういったデータから、データのクオリティを判断し、
以降の解析において、

全体的にクオリティーが高いためそのままデータを用いる
クオリティーの低い末端をトリミングしたデータを用いる
全体的にクオリティーが低いため、シーケンスをかけ直す
といった行動の指標にすることができます。

つまり、FastQCによりクオリティーレベルを確認して次のステップを決めるわけです。今回はデモなのである程度のクオリティーが確認できれば、de novoアセンブリに使います。Graphical User Interfaceでも利用できるソフトウェアなので、気軽にいきましょう。

FastQCはhomebrewを利用してインストールできます。以下コマンドをフォローしてください。

iMac:~ ryoniwa$ brew search fastqc
iMac:~ ryoniwa$ brew install fastqc

インストールが完了したところで、FastQCの使い方を説明します。使い方は2種類あります。一つはGraphical User Interfaceでポチポチクリックして進める方法、もう一つはコマンドで処理して進める方法です。今回は両方紹介します。

GUIでやる方法

FastQCを呼び出すのにコマンドを使います。

iMac:~ ryoniwa$ fastqc

すると、FastQCが立ち上がります。下記写真の通り、メニューバーからFile => Openを開き、ファイルを選択します。すると自動で解析が始まります。
fastqc.png

CUIでやる方法

コマンドでファイルを指定して、htmlファイルを吐き出させます。--nogroupというオプションを使うと解析量が増えるので解析としては重くなりますが、各塩基ポジションごとの解析が見れます。要は--nogroupをつけると各塩基ごとにデータが出てくる一方で、つけないと201~210塩基のところはこれぐらいのクオリティーのようにまとめられます。用途に応じて変えてください。今回はどちらでもいいです。-tは使用するスレッド数です。

iMac:Training ryoniwa$ fastqc -t 4 -o fastqc_results/ データ.fastq
iMac:Training ryoniwa$ fastqc -t 4 -o fastqc_results/ --nogroup データ.fastq

この方法でやると自動的にfastqc_resultsというフォルダがTrainingフォルダに作られます。htmlファイルが作られますので、Google Chromeなどに読み込ませると見やすいです。

結果の確認

解析が終わると、データが確認できます。確認の方法はhttps://bi.biopapyrus.jp/rnaseq/qc/fastqc.htmlに載っています。様々な項目がありますが、注意して確認すべき項目はBasic Statistics、Per base sequence quality、Per sequence quality scoresです。

Per base sequence qualityは上記サイトによれば

リードの各位置におけるクオリティスコアの分布が示されている。横軸がリード上の位置を表し、縦軸はクオリティスコアを表す。スコアの分布はボックスプロットによって描かれる。

スコアの下側四分位点(1/4 分位点)が 10 未満または中央値が 25 未満の位置が存在する場合 Warning とされる。また、スコアの下側四分位点(1/4 分位点)が 5 未満または中央値が 20 未満の位置が存在する場合 Warning とされる。

を表記します。したがって、今回の場合クオリティスコアの分布がQ25を上回っていれば、きちんとクオリティーカットされていることになります。

Per sequence quality scoresは上記サイトによれば

各クオリティスコアを持つリードがどれぐらい存在するかを示している。横軸がクオリティスコアを示し、縦軸はリード数を示す。スコアの高いところにリード数が多くなることが相応しいと考えられる。

リード数の最も多いスコアが 27 未満ならば Warning とされる。また、20 未満ならば Failure とされる。

を表記します。したがって、今回の場合クオリティスコアの分布がQ25を下回るリードがなければ、きちんとクオリティーカットされていることになります。

自分の決めたThresholdを下回るリードが存在しなければ、信頼できるデータとするというのが基本的な考え方だと理解しています。他にも項目があります。FastQCは幸いにもググるとたくさん情報が出てきますので、気になる人は自学してください。

ここでデータを確認し、クオリティートリミングが上手くできたとして次に進めます。

SPAdes

http://kazumaxneo.hatenablog.com/entry/2018/06/13/232022に詳しく書かれています。また、論文もあります。当方は難しかったので半分しか読んでおりません。

Bankevich A., Nurk S., Antipov D., Gurevich A., Dvorkin M., Kulikov A. S., Lesin V., Nikolenko S., Pham S., Prjibelski A., Pyshkin A., Sirotkin A., Vyahhi N., Tesler G., Alekseyev M. A., Pevzner P. A. SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing. Journal of Computational Biology, 2012

SPAdesはde brujin graphと呼ばれるアルゴリズムを搭載したde novoアセンブリソフトウェアになります。De Brujin Graphを搭載したことがいかにすごいかはこのサイトを見れば分かります。我々が頭の中で、ゲノムアセンブルってこんな感じかなと想像するのは別の方法でOverlap Layout Censensus法 (OLC法) と呼ばれます。両者の違いは根底にある計算方法です。

OLC法ではハミルトンパス問題が使われるのに対して、De Brujin Graph法ではオイラーパス問題が使われています。当方も詳しいことは分かりませんが、コンピュータの計算上ではオイラーパス問題の方が極めて都合がいいそうです。違いとしてハミルトンパス問題では設定した領域分共通しているリードをつなげて長くするという考え方ですが、オイラーパス問題ではリードを細かくして、繋げれるだけ繋げて、最終的に一番長く繋がるルートを探すという考え方になります。

SPAdesは複数のK-mer (リードを細かくする際の長さ) を用意し、同時に解析させて、最終的に混ぜて長く繋がるように設計されているようです。多くのゲノムアセンブラは更新が止まっていますが、SPAdesは未だ更新が続いています。そういう意味では信頼できるソフトウェアだなと考えています。

SPAdesのインストール

SPAdesもHomebrewでインストールできます。

$ brew search spades
$ brew install spades
$ spades -h

SPAdesの使い方

SPAdes genome assembler v3.13.0

Usage: /usr/local/bin/spades.py [options] -o <output_dir>

Basic options:
-o  <output_dir>    directory to store all the resulting files (required)
--sc            this flag is required for MDA (single-cell) data
--meta          this flag is required for metagenomic sample data
--rna           this flag is required for RNA-Seq data 
--plasmid       runs plasmidSPAdes pipeline for plasmid detection 
--iontorrent        this flag is required for IonTorrent data
--test          runs SPAdes on toy dataset
-h/--help       prints this usage message
-v/--version        prints version

Input data:
--12    <filename>  file with interlaced forward and reverse paired-end reads
-1  <filename>  file with forward paired-end reads
-2  <filename>  file with reverse paired-end reads
-s  <filename>  file with unpaired reads
--merged    <filename>  file with merged forward and reverse paired-end reads
--pe<#>-12  <filename>  file with interlaced reads for paired-end library number <#> (<#> = 1,2,...,9)
--pe<#>-1   <filename>  file with forward reads for paired-end library number <#> (<#> = 1,2,...,9)
--pe<#>-2   <filename>  file with reverse reads for paired-end library number <#> (<#> = 1,2,...,9)
--pe<#>-s   <filename>  file with unpaired reads for paired-end library number <#> (<#> = 1,2,...,9)
--pe<#>-m   <filename>  file with merged reads for paired-end library number <#> (<#> = 1,2,...,9)
--pe<#>-<or>    orientation of reads for paired-end library number <#> (<#> = 1,2,...,9; <or> = fr, rf, ff)
--s<#>      <filename>  file with unpaired reads for single reads library number <#> (<#> = 1,2,...,9)
--mp<#>-12  <filename>  file with interlaced reads for mate-pair library number <#> (<#> = 1,2,..,9)
--mp<#>-1   <filename>  file with forward reads for mate-pair library number <#> (<#> = 1,2,..,9)
--mp<#>-2   <filename>  file with reverse reads for mate-pair library number <#> (<#> = 1,2,..,9)
--mp<#>-s   <filename>  file with unpaired reads for mate-pair library number <#> (<#> = 1,2,..,9)
--mp<#>-<or>    orientation of reads for mate-pair library number <#> (<#> = 1,2,..,9; <or> = fr, rf, ff)
--hqmp<#>-12    <filename>  file with interlaced reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)
--hqmp<#>-1 <filename>  file with forward reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)
--hqmp<#>-2 <filename>  file with reverse reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)
--hqmp<#>-s <filename>  file with unpaired reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)
--hqmp<#>-<or>  orientation of reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9; <or> = fr, rf, ff)
--nxmate<#>-1   <filename>  file with forward reads for Lucigen NxMate library number <#> (<#> = 1,2,..,9)
--nxmate<#>-2   <filename>  file with reverse reads for Lucigen NxMate library number <#> (<#> = 1,2,..,9)
--sanger    <filename>  file with Sanger reads
--pacbio    <filename>  file with PacBio reads
--nanopore  <filename>  file with Nanopore reads
--tslr  <filename>  file with TSLR-contigs
--trusted-contigs   <filename>  file with trusted contigs
--untrusted-contigs <filename>  file with untrusted contigs

Pipeline options:
--only-error-correction runs only read error correction (without assembling)
--only-assembler    runs only assembling (without read error correction)
--careful       tries to reduce number of mismatches and short indels
--continue      continue run from the last available check-point
--restart-from  <cp>    restart run with updated options and from the specified check-point ('ec', 'as', 'k<int>', 'mc', 'last')
--disable-gzip-output   forces error correction not to compress the corrected reads
--disable-rr        disables repeat resolution stage of assembling

Advanced options:
--dataset   <filename>  file with dataset description in YAML format
-t/--threads    <int>       number of threads
                [default: 16]
-m/--memory <int>       RAM limit for SPAdes in Gb (terminates if exceeded)
                [default: 250]
--tmp-dir   <dirname>   directory for temporary files
                [default: <output_dir>/tmp]
-k      <int,int,...>   comma-separated list of k-mer sizes (must be odd and
                less than 128) [default: 'auto']
--cov-cutoff    <float>     coverage cutoff value (a positive float number, or 'auto', or 'off') [default: 'off']
--phred-offset  <33 or 64>  PHRED quality offset in the input reads (33 or 64)
                [default: auto-detect]

SPAdesは設定が極めて複雑です。パラメータが多い理由はそれだけゲノムアセンブリのパラメータ設定が肝要であるかを表していると考えていいと思います。少し数字を触るだけで、かなりデータが変わります。今回はトライアルなので、難しく考えずに以下に当てはまるように書き換えてください。output_dirは現在使用しているディレクトリ上に出来ます。

$ spades.py -1 paired_output_1.fq -2 paired_output_2.fq -k auto --careful -o output_dir

万が一、PacBioやNanoporeを含むデータを選んだ場合は--pacbioや--nanoporeのオプションを利用してください。Short readシーケンスで足りない部分を補えるので、確かにロングリードを入れるアセンブリの方が早くコンプリートになりやすいです。解析が終了すると、たくさんのファイルがoutput_dirに出来ます。そのうち、scaffolds.fastaが最終生産物になります。

続いて、scaffolds.fastaがどういう状態なのかを確認していきます。

QUAST

QUASTはアセンブリの性能や精度を評価するツールになります。つまり、どれくらいアセンブルできたかどうかを簡単に確認できます。今回であれば、SPAdesを利用してアセンブルしていますが、このツールを利用することでどれくらいコンティグができたかなどが確認できます。

QUASTの使い方

QUASTもFastQC同様、GUIとCUIの両方で利用できます。QUASTに関してはGUIで利用するのが楽かと思います。

GUIでQUASTを行う。

まずはhttp://quast.bioinf.spbau.ruを開きます。
データをアップロードして、メールアドレスを打ち、しばらく待ってください。大抵10分くらいで終わります。
リファレンスゲノムがあると非常に良いですが、今回はトライアルなので省きます。

CUIでQUASTを行う。

QUASTはhomebrewでインストールできます。

$ brew install quast
$ quast scaffolds.fasta (-R reference.fasta) -o results -t 4

resultsフォルダがCurrent Directoryにできます。その中にPDFファイルとHTMLファイルができます。

結果を確認する

http://kazumaxneo.hatenablog.com/entry/2017/08/22/231007を見てください。結果については非常に分かりやすいかと思います。説明を省きます。

DFASTによるアノテーション

https://dfast.nig.ac.jp/dfc/を見てください。
動画の場合はhttps://togotv.dbcls.jp/20180318.htmlをご覧ください。十分に説明があるので、あまり書きませんがDFASTはWeb上で使えるアノテーションパイプラインです。他のアノテーションソフトのProkkaなどに比べて、データが選べる、データが見やすいので当方はDFASTを好んで使っています。Prokkaも使ってますが、DFASTの方がしっかりアノテーションしてくれる感じがしています。

使い方

説明するまでもなく、超簡単です。https://dfast.nig.ac.jp/dfc/を開いてアップロード、メールアドレスを打って、待つだけです。コンピュータに強ければ、dfast-coreを入れてもらうのがいいと思います。For文を使うと同時に多サンプルが処理できますので、当方はdfast-coreしか使っていません。dfast-coreにやり方載っています。

結果の見方

調べたい遺伝子を探すもよし、SnapGeneやIGVを使って環状ゲノムマップを書くもよし、もう少し勉強してComparative Genomicsするもよし。ここから先がサイエンスです。(小並)

終わりに

本講義の内容は以上になります。お疲れ様でした。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした