（2016-03-04 暫定版です）
前回の続きです。Illumina NGS (HiSeq, MiSeqなど）のデータとPacBioのデータの両方を使ってゲノムのde novo assembleを行う (hybrid assembly)方法について書きます。
ここではDGB2OLCのManualに示されているサンプルデータを用いて解析を行います。以下は解析の手順（ほぼマニュアルの翻訳）になります。
メモ：Step 1（Illumina readsのassemble）とStep 2（PacBioとIllumina contigのhybrid assemble）では、blasrとSparcは不要です。blasrがインストールできなくても、hybrid assembleのcontig（backbone_raw.fasta）を得ることはできます。Step 3（consensus配列の作成）にはblasrが必要になります。
サンプルデータのダウンロード（FTPサイト）
Illumina reads (fastq): ftp://qb.cshl.edu/schatz/ectools/w303/Illumina_500bp_2x300_R1.fastq
PacBio reads (fasta): ftp://qb.cshl.edu/schatz/ectools/w303/Pacbio.fasta.gz
Step 0: 下準備
まず、解析を行う適当なディレクトリを用意して、そのディレクトリに移動します。サンプルデータ（Illumina, PacBio）はそのディレクトリに置きます。
次に、実際に解析を行うディレクトリを1つ作成します。解析の結果、たくさんのファイルが出力されるので、適当なディレクトリを作っておかないと混乱します。
# 解析を行うディレクトリ // 日付 + sample mkdir 160304_sample # カレントディレクトリを移動 cd 160304_sample
Step 1: Illumina readsをassembleする
最初に、Illumina readsをassembleしてcontigを作成します。Assembleはどのようなソフトで行っても良いようですが（ようするに、contigのfastaファイルがあれば良い）、DBG2OLCが用意しているサンプルデータでは、SparseAssemblerを使用しています。
以下のコマンドを実行します。詳しくはManualを参照。
# SparseAssemblerによるIllumina readsのde novo assemble SparseAssembler LD 0 k 51 g 15 NodeCovTh 1 EdgeCovTh 0 GS 12000000 f ../Illumina_500bp_2x300_R1.fastq
解析がうまくいけば、いくつかのファイルができます。Contigのファイルは"Contigs.txt"になります。
Step 2: DBG2OLCを実行する
PacBio readsとIllumina readsから作られたcontigsを組み合わせたhybrid assemblyを行います。DBG2OLCを使用します。いくつかの変数（パラメーター）を設定する必要があります。
以下のようなコマンドで実行します。
# step2 // DBG2OLC k 17 AdaptiveTh 0.0001 KmerCovTh 2 MinOverlap 20 RemoveChimera 1 Contigs Contigs.txt f ../Pacbio.fasta
サンプルの例ではPacBioのデータが1つのファイルしかないですが、もし複数のsubreadsのファイルがある、また複数cell分のFASTA fileが存在する場合、catなどを使って1つのファイルに合体して使えば問題ないようです。
# たとえば、3つのFASTA fileを1つに合体する cat subread1.fasta subread2.fasta subread3.fasta > PacBio.fasta
以下、変数の説明です。
# 3つのメジャーな変数（あんまりわかってません。。すいません）
| AdaptiveTh | adaptive k-mer matchingのしきい値。もし"M < AdaptiveTh × Contig_Length"ならば、 そのcontigはlong readのanchorとしては使えない。（Mはcontig (Illumina)とlong read (PacBio)の間でmatchするk-mers） |
| KmerCovTh | fixed k-mer matchingのしきい値。あるcontigについて、もし"M < KmerCovTh"ならば、そのcontigはlong readのanchorとしては使われない。 |
| MinOverlap | Long readsのpair間で（許容する）最小のoverlap scoreの値。 |
Manualの[Miscellaneous]のところには、これらのパラメーター設定について、以下のような記述があります。
[Miscellaneous]
At this point, the parameters may be fine-tuned to get better performance. As with SparseAssembler, LD 1 can be used to load the compressed reads/anchored reads.
Suggested tuning range is provided here:
For 10x/20x PacBio data: KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01.
For 50x-100x PacBio data: KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02.（和訳）この点について、これらの（3つの）パラメーターを最適化することで、より良いパフォーマンスが得られる可能性がある。Compressed read/anchored readsを読み込むため、SparseAssemblerではLD 1にする設定が利用できる（<-- よくわかりません。すいません）。
10x/10x PacBio data: KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01にする
50x-100x PacBio data: KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02にする
また、実際に変数をどのように設定するかは「いくつかのパターンで実行してみて、より良い結果を出すものを探すべし（意訳）」とのことです。
その他の、あまりflexibleでない（重要性は相対的に低い）パラメーターについて、
| k-mer size | 17が良い |
| Contigs | Illumina readsのassmbleにより得られたcontigのファイル (FASTA) |
| MinLen | 最小のread length（しきい値？） |
| RemoveChimera | （たぶんPCR artefactで生じた）キメラのリードを除く。データのcoverageが10x以上の場合、1に設定する |
Step 3: ゲノムの「コンセンサス配列」を作成する
ここでblasrとSparcが必要になります。blasrおよびSparcのインストールについては前のエントリを参照して下さい。
(1) Sparcのプログラムおよびいくつかのスクリプトを、解析を行うディレクトリにコピーする
コピーでもシンボリックリンクでもいいですが、Sparcフォルダの中にある以下のものを解析を行うディレクトリに置きます。
compiled/Sparc utility/RunSparcConsensus.txt utility/SeqIO.py utility/split_and_run_sparc.sh utility/split_reads_by_backbone.py
(2) ctg_pb.fastaを作成する
"ctg_pb.fasta"はIllumina contigのファイルとPacBioのファイルを合体したファイルで、次の作業のinput fileになります。
cat Contig.txt ../PacBio.fasta > ctg_pb.fasta
(3) シェルスクリプトsplit_and_run_sparc.shを実行する
以下のコマンドを実行して、解析開始です。Consensus配列の結果は、"results"というフォルダ内にできます。
sh split_and_run_sparc.sh backbone_raw.fasta DBG2OLC_Consensus_info.txt ctg_pb.fasta results 2>cns_log.txt
ここでのinput fileは以下の4つです
(1) backbone_raw.fasta (by DBG2OLC)
(2) DBG2OLC_Consensus_info.txt (by DBG2OLC)
(3) DBG contigs (in fasta format)
(4) PacBio reads (in fasta format)
3と4を合体したのが"ctg_pb.fasta"になります。
