ハイスループットDNAシーケンス技術の中で、Solexa / Illuminaプラットフォーム[ref.1]は、1回の実行で最大量のシーケンスデータを作成する[ref.2]。この技術の1つのユニークな特質は、与えられたDNA分子の両端からシーケンスリードを生成するその能力である。これはbiologicalな解釈のための多くの機会を提供する。例えば、ペアエンドリードをリファレンス配列とアライメントさせることによって、その全体をシーケンシングすることなく、ＤＮＡ分子の全範囲を推測することができる。
イルミナのシーケンスラン出力は、一連のFASTQファイルである。このファイルには、リード配列と対応するクオリティスコアが含まれている[ref.3]。 サンガーシーケンシングで最初に開発されたように、クオリティスコアは、以下の式［ref.４］を介して、所与のシーケンシングされた塩基が誤っているという確率から決定される。
したがって、クオリティスコアが40だと10,000分の1の確率で誤りがあるはずである。配列の変異検出では、偽陽性を減らすために[ref.5]最小のクオリティスコアを達成する塩基のみを考慮するプログラムが多く、ガン診断などの臨床変異検出では、公表されたガイドラインはしばしばそのような基準を組み入れている[ref.6、7]。 。しかしながら、イルミナのマシンによって生み出されたrawクオリティスコアは膨らんでいることが多くの研究によって示されている。すなわち、所与のクオリティスコアを有するシーケンシングされた塩基は、式１（ref.1）から予想されるよりも高いエラー率を有する、特にスケールの上限で[ref.8, 9, 10, 11]。
イルミナのシークエンシング技術の進歩によりリード長が増し、特定のライブラリーのペアエンドリードがかなりのシークエンスオーバーラップを持つようになっている。これらのオーバーラップした領域からＤＮＡ断片の全長をスパンした単一のリードにマージすることが可能である。このプロセスはエラーの訂正およびリードカバレッジの正確な決定を可能にし、そしてそれはターゲットバリアントリシーケンシング[ref.12]からメタバーコーディング（例えば16S / 18S rRNA研究）[ref.13]までの範囲の応用においてますます認められるようになった。
初期のマージプログラムはfastq-join [ref.14]であり、これは微生物生態分析パッケージQIIME [ref.15]のデフォルトオプションだった。そのパッケージの最新版、QIIME2は、オープンソースのVSEARCH [16]をマージに使用している。 3番目のマージプログラムであるPEAR [ref.17]は他の2つのプログラム（FLASH [ref.18]やCASPER [ref.19]のような他のプログラム）より大きなアドバンテージを持っており、"ありつぎ”のアライメントを行う。すなわち1つのリードの3'末端を延長してペアの５ '末端まで伸ばす（論文 図１ｂ参照 link）。リード長よりも短いDNA断片のシーケンシングでは、リードの3 '末端にアダプターの部分を含むリードをもたらすだろう。このようなリードペアは、dovetailed alignments（延長アライメント）を考慮していないプログラムではマージされない[ref.17]。
マージプログラム間のもう1つの重要な違いは、マージされたベースにクオリティスコアを割り当てる方法である。 Fastq-joinとFLASHは、2つリードの塩基がマッチした場合、より高い塩基の方のクオリティスコアが使用されるという単純な方式を使用する。塩基が一致しない場合は、クオリティスコアの高い方の塩基が選択され、クオリティスコアの差がマージされたスコアになる。PEARは、元のR1とR2の塩基が一致する場合、マージされたベースのクオリティスコアは両方の元の塩基が間違っている確率を反映するという推論で設計されている。
したがって、PEARはマッチする塩基のクオリティスコアを合計する。 VSEARCHは EdgarとFlyvbjerg [ref.20]によって開発されたより洗練されたモデルを利用するが、塩基をマッチングするために結果として得られるスキームはPEARのものとほとんど同じになる（論文追加ファイル1：図S1参照）。たとえば、クオリティスコアが40のマッチする塩基からマージされた塩基には、fastq-joinでは40がアサインされ、PEARでは80、VSEARCHでは85のクオリティスコアがアサインされる（プログラムによって設定された人工の上乗せクオリティスコアは無視される）。
欠点の可能性があるにもかかわらず、これらのクオリティスコアプロファイルのいずれも経験的にテストされていない。たとえば、PEARとVSEARCHのプロファイルは、式2に基づいている。 上記の1式は、上記に記載のように、イルミナシーケンシングでは不正確であることが実証されている。さらに、PEARとVSEARCHはどちらも、2回のリードでのシーケンスエラーが独立しているという前提のもとに設計されている。つまり、上記の式3は、2つのイベント（「R1の塩基が間違っている」、と「R2の塩基が間違っている」）が独立している場合のみ式2に従う。この仮定は、我々の知る限りでは検証されていない。
この論文では、まず、腸内細菌ファージΦX174「PhiX」から得られたリードを使用して、Illuminaのマシンによって作成されたベースラインクオリティスコアとエラー率の関係を評価する。このウイルスのゲノムは最初にシーケンシングされたゲノムDNAであり[ref.21]、コントロールとしてPhiX DNAの断片のライブラリが日常的にイルミナのシーケンシングランに加えられている。 PhiXリードの幅広いアクセス可能性に加えて、このライブラリのフラグメントのサイズは、長いIlluminaの実行によって生成されたほとんどのPhiX派生リードペアが十分なオーバーラップを持ち、よってマージする塩基のクオリティスコアプロファイルを作成するために使用できる。この作業を実行し、これらのプロファイルを新しいマージプログラムNGmergeに組み込んだ（図1）。 NGmergeが他のマージプログラムよりもエラーとあいまいな塩基（Ns）を正しく訂正し、塩基のエラー率を正確に反映したクオリティ値を持つマージリードを生成することを示す。
インストール
mac os10.14のminiconda3-4.3.30環境でテストした。
依存
本体 Github
#Anacondaを使っているならcondaで導入可能
conda install -c bioconda -y NGmerge
> NGmerge
$ NGmerge
Error! Need input/output files
Usage: NGmerge {-1 <file> -2 <file> -o <file>} [optional arguments]
Required arguments:
-1 <file> Input FASTQ file with reads from forward direction
-2 <file> Input FASTQ file with reads from reverse direction
-o <file> Output FASTQ file(s):
- in 'stitch' mode (def.), the file of merged reads
- in 'adapter-removal' mode (-a), the output files
will be <file>_1.fastq and <file>_2.fastq
Alignment parameters:
-m <int> Minimum overlap of the paired-end reads (def. 20)
-p <float> Mismatches to allow in the overlapped region
(a fraction of the overlap length; def. 0.10)
-a Use 'adapter-removal' mode (also sets -d option)
-d Option to check for dovetailing (with 3' overhangs)
-e <int> Minimum overlap of dovetailed alignments (def. 50)
-s Option to produce shortest stitched read
I/O options:
-l <file> Log file for stitching results of each read pair
-f <file> FASTQ files for reads that failed stitching
(output as <file>_1.fastq and <file>_2.fastq)
-c <file> Log file for dovetailed reads (adapter sequences)
-j <file> Log file for formatted alignments of merged reads
-z/-y Option to gzip (-z) or not (-y) FASTQ output(s)
-i Option to produce interleaved FASTQ output(s)
-w <file> Use given error profile for merged qual scores
-g Use 'fastq-join' method for merged qual scores
-q <int> FASTQ quality offset (def. 33)
-u <int> Maximum input quality score (0-based; def. 40)
-n <int> Number of threads to use (def. 1)
-v Option to print status updates/counts to stderr
実行方法
ペアエンドfastqをマージする
NGmerge -1 pair_1.fq.gz -2 pair_2.fq.gz -o merged.fq.gz
詳細はGithubで確認して下さい。
引用
NGmerge: merging paired-end reads via novel empirically-derived models of sequencing errors
John M. Gaspar
BMC Bioinformatics. 2018; 19: 536
関連
