近年、配列決定されたゲノムの数と多様性が非常に増加している（論文より Reddy et al、2015）。 13,000以上の真核生物が配列決定されているか、配列決定の過程にあり、数百の植物や動物を含むより多くのものが計画されている。大部分のモデル生物は、高品質の参照配列（Arabidopsis Genome Initiative、2000; Lander et al、2001; Venter et al、2001; Waterston）を創出する目的で、幅広い国際コンソーシアムによって行われた大きなゲノムプロジェクトの傘下で配列決定されている。例えば 2002;ラットゲノム配列決定プロジェクトコンソーシアム、2004;国際ヒトゲノムシーケンシングコンソーシアム）。近年、第2世代シークエンシング技術はシーケンシングへのアクセスの増加に伴うシーケンシングコストの削減により、グループや個々の研究者が研究する種のゲノムを配列することを可能にしたため、新しいゲノムの生成のスピードを劇的に加速している。これらのプロジェクトのほとんどすべてが、ゲノムのドラフト版を作成する。
典型的なゲノムアセンブリプロセスの間、重複したリードを使用してコンティグを構築し、次にコンティグをメイトペアリードからの順序および方向情報を用いてより大きなスキャホールドに繋いでいく。ゲノム中の反復領域は、アセンブリアルゴリズムに重大な問題を提起する。これらのシーケンスを再構築できるようにするには、リピートのある隣接領域にアセンブリを固定できるように、メイトペアのインサートサイズがリピートの長さを超える必要がある。したがって、典型的なゲノムアセンブリプロジェクトでは、500 bpから8-10 kbに及ぶ複数のインサートサイズライブラリが必要になる。しかし、重要なロング・インサート・ライブラリーを構築することは、高価で労力を要する。ドラフトゲノムが生成されると、その解析における最初の最も重要なステップは、遺伝子機能および変異の下流の研究のための基礎を提供する遺伝子の検出である。 Deep RNA-seqは種の遺伝子を特徴付ける主要な手段となっており、アラインメントおよび転写産物のアセンブルツールを含むRNA-seq解析のための多数のツールが既に存在する（Wang et al、2010; Trapnell Floren and Salzberg 2013; Kim et al、2013; Mezlini et al、2013; Song and Florea 2013; Kim et al、2015）。また、種特異的遺伝子および選択的にスプライシングされた変異体についての重要な情報を提供し、新規タンパク質コード遺伝子および非コードRNAを含む。しかし、アセンブリのエラーやギャップは、遺伝子を断片化したり、エキソンを削除したり、遺伝子配列を局所的に変更したりすることによって、正確な遺伝子や転写注釈を妨げる可能性がある（Florea et al、2011）。したがって、特に遺伝子領域におけるアセンブリの品質を向上させるためにあらゆる努力がなされなければならない。
ヒトゲノムから始まるいくつかのゲノム配列決定プロジェクトは、アセンブリエラーを検出するため、またはアセンブリへのさらなるコンティグを補充するために、独立に生成された発現配列タグおよび全長メッセンジャーRNA（mRNA）からの遺伝子情報を使用している（Venterら、2001; Dalloul et a、2010; Shulaev et al、2011; Wegrzyn et al、2014）。しかし、ゲノムプロジェクトに体系的に適用され、生成される次世代シーケンシングデータを活用するツールは欠けている。伝統的なscaffolding法（Pop et al、2004; Dayarianら、2010; Boetzerら、2011; Gaoら、2011; Donmez and Brudno 2013）は、カバレッジの一様な分布に頼っており、RNA-seqデータに直接適用することはできない。 L_RNA_scaffolder（Xue et al、2013）は、Trinity（Haas et al。、2013）のようなツールで生成されたde novo RNAアセンブリを用いて遺伝子ベースのアプローチを適用している）、AGOUTI（Zhang et al。、2015）は、既知の遺伝子アノテーションの文脈でRNA-seqリードのアライメントを直接使用して新しい接続を検出する。しかしながら、L_RNA_scaffolderによって使用されるようなRNA-seq配列のデノボアセンブリーは困難であり、エラーを起こしやすく、また時間がかかる。キメラ転写産物の再構築は誤った骨格を導き、低発現遺伝子は部分的にのみ再構成または欠損する可能性があり、したがって影響は限られている。
RNA-seqリードのアライメントを使用して断片化ゲノムの新しいコンティグ連結を同定し、遺伝子とゲノムの完全性と正確性を同時に改善する新しいツールであるRascaf（RnA-SCAFfolder）を開発した。Rascafはエクソンブロックグラフを用いて遺伝子とその下のコンティグ関係を同時に表現し、最も重いコンティグパスを決定する。提案されたコンティグ連結は、種間の相補DNA（cDNA）およびタンパク質の証拠をデータベース検索することによって任意選択で確認できる。シミュレートされたデータとリアルデータの両方で評価され、同様のツールと比較し、Rascafはより正確で効率的であり、植物や動物の新しいゲノムアセンブリの品質を高めるために効果的に使用できることが示された。
Rascafは精度が非常に高く、誤ったアセンブリがわずかしか導入されれず、メモリの占有スペースが小さく、典型的なRNA-seqデータセットであれば一般的なワークステーションマシンで数分で動作する。Rascafは、いくつかのFragaria種およびRosaceae（Pyrus communis L.）のドラフトゲノムにおいて、1000〜10,000個の新規コンティグ連結を同定したことが示されている。
インストール
https://github.com/mourisl/Rascaf
git clone https://github.com/mourisl/rascaf.git
cd rascaf
make
> ./rascaf
r$ ./rascaf
usage: rascaf [options]
options:
-b STRING (required): the path to the coordinate-sorted alignment BAM file
-f STRING (recommended): the paths to the raw assembly fasta file(default: not used)
-o STRING : prefix of the output file (default: rascaf)
-bc STRING: the path to the alignment BAM file allowing clipping from non-spliced aligner (default: not used)
-ms INT: minimum support for connecting two contigs(default: 2)
-ml INT: minimum exonic length(default: 200)
-breakN INT: the least number of Ns to break a scaffold in the raw assembly (default: 1)
-k INT: the size of a kmer(<=32; <=0 if you do not want to use kmer. default: 23)
-cs : output the genomic sequence involved in connections in file $prefix_cs.fa (default: not used)
-v : verbose mode (default: false)
> ./rascaf-join
$ ./rascaf-join
usage: rascaf_join [OPTIONS]
OPTIONS:
-r STRING: the path to the rascaf output. Can use multiple of -r. (required)
-o STRING: the prefix of the output file. (default: rascaf_scaffold)
-ms INT: minimum support alignments for the connection (default: 2)
-ignoreGap: ignore the gap size, which do not consider the number of Ns between contigs (default: not used)
パスの通ったディレクトリにコピーしておく。
ラン
ランは２段階で行われる。
Githubでは、２つのRNA-seqデータA、Bを使いscaffoldsを作る例が示されている。まずAのfastqとBのfastqそれぞれを、別のDNAアセンブラでアセンブルして得たcontigファイルにマッピングしてbamを作る（マッピングには真核生物ならtophat2やhisat2を使う）。
それから、シングルRNA-seqデータからつながりを検出するrascafコマンドでA.bamとB.bamを解析する。
rascaf -b A.bam -f assembly.fa -o A
rascaf -b B.bam -f assembly.fa -o B
rascafの出力をマージし、検出された繋がりからscaffoldsを出力する。
rascaf-join -r A.out -r B.out -o assembly_scaffold
--テストラン--
サンプルデータを使ってみる。
rascaf -b sample/sample.bam -f sample/sample.fa
#rascaf.outができる。
rascaf-join -r rascaf.out -o output
#outout.faができる。
元のsample.fa ↓
$ seqkit stats sample/sample.fa
file format type num_seqs sum_len min_len avg_len max_len
sample/sample.fa FASTA DNA 2 200,000 100,000 100,000 100,000
rascafを使った後のoutput.fa ↓
$ seqkit stats output.fa
file format type num_seqs sum_len min_len avg_len max_len
sample.fa FASTA DNA 1 200,017 200,017 200,017 200,017
perlのラッパースクリプトで２段階ステップをワンライナーでランできるようです。詳細はgithubで確認してください。
別に紹介したP_RNA_scaffolderも確認してみてください。
引用
Rascaf: Improving Genome Assembly with RNA Sequencing Data.
Song L, Shankar DS, Florea L.
Plant Genome. 2016 Nov;9(3).
