Pacific Bioscience（PacBio）が20111年に1分子リアルタイム（SMRT）シーケンス技術を商品化し、第3世代シーケンシング（TGS）が登場した。TGSプラットフォームには大きな技術的違いがある。これは、第2世代シーケンス（SGS）とは異なる。ペアエンド情報を考慮すると、メインのSGSプラットフォームであるイルミナは、各DNAフラグメントから50〜600 bpの情報を提供する。最長SGSリードを生成する454シーケンス（〜700bp）を含め、> 1000bpを提供するSGSプラットフォームはない。したがって、シーケンシング長が短いため、SGSの適用は遺伝子アイソフォームの再構築などの大規模または複雑なゲノムイベントに限定される。 TGSはリード長の長さによってこれらの困難な問題を解決する。
最も広く使用されているTGSプラットフォーム[PacBioおよびOxford Nanopore Technologies（ONT）]は、１DNA分子から非常に長いヌクレオチド配列を直接捕捉するための新しい生化学/生物物理学的方法を開発した。他の新たなTGSプラットフォーム（Moleculoおよび10X Genomics）は、同じDNA分子からのショートリードのアセンブリによる合成ロングリード（SLR）を生成することに基づいている。ここでは、PacBioとONTのデータ機能とそれらのトランスクリプトーム解析への応用に焦点を当てる。
PacBioは1DNA分子を捕獲し、Illuminaは増幅されたDNAフラグメントのクローン集団から増強されたシグナルを検出することを除いて、イルミナのシークエンシングと同様の合成によるシークエンシング戦略を採用する。raw PacBioデータのエラー率は13〜15％である。これは、1DNA分子からのS / N比が高くないためである（ref.3)。正確性を高めるために、PacBioプラットフォームはヘアピンアダプターを標的二本鎖DNAの両端に連結することにより環状DNAテンプレートを使う。ポリメラーゼが環状分子を繰り返し横断しそして複製するにつれて、DNA鋳型は連続的な長いリード（CLR）を生成するために複数回シーケンシングされる。 CLRは、アダプター配列を除去することによって複数のリード（「サブリード」）に分割でき、複数のサブリードはより高い精度で循環コンセンサス配列（「CCS」）リードを生成する。 CLRの平均長は> 10kbから最大60kbである。これはポリメラーゼの寿命による( ref.3)。したがって、CCSリードの長さと精度はフラグメントサイズによって異なる。PacBioシークエンシングは、ゲノム（例えば、新規のアセンブリ、構造変異の検出およびハプロタイプ決定）およびトランスクリプトーム（例えば、遺伝子アイソフォームの再構築および新規遺伝子/アイソフォームの発見）の研究に利用されてきた。
ヌクレオチド取り込みまたはハイブリダイゼーションを利用する他のシーケンシング技術と比較して、ONTは、特徴的な電流を測定することによって天然の一本鎖DNA（ssDNA）分子を直接シーケンシングする。塩基が分子モータータンパク質によってナノポアを貫通するにつれて変化する。ONT MinIONは、PacBio環状DNAテンプレートと同様のヘアピンライブラリー構造を使用する。DNAテンプレートとその相補鎖はヘアピンアダプターによって結合されている。したがって、DNAテンプレートはナノポアを通過し、続いてヘアピンが通過し、最後に補体が通過する。 rawリード、アダプターを取り外すことによって2つの「1D」リード（「テンプレート」と「補数」）に分割できる。 2回の「1D」リードのコンセンサスシーケンスは、より高い精度での「2D」リードである（ref.2）。 PacBioと同様のデータ機能を備えているため、多くの研究者がPacBioが適用されているアプリケーションでONTを利用したりテストしたりしている。
PacBioとONTプラットフォームはリード長が長いという利点を共有しているが、SGSと比較してシーケンスエラー率が高く、スループットが低いという同じ欠点もある(ref.3、14–16)。高いシーケンシングエラー率は、転写物の正確なシーケンシング、スプライス部位の同定およびSNP callingのような1ヌクレオチド分析に課題を投げかける。低処理量は、遺伝子／アイソフォーム存在量の推定などの定量分析にとって障害となる。 PacBio CCSとONT 2Dのコンセンサス戦略はエラー率を減らすことができるが、対応するリード長は短くなり、スループットは低下する。したがって、TGSデータとSGSデータとを統合するハイブリッドシーケンシング（「Hybrid-Seq」）は、SGSデータの支援によるTGSデータの分析に関連する制限に対処するためのアプローチとして浮上している。たとえば、PacSioまたはONTのロングリードをSGSのショートリードでエラー訂正すると、ロングリードの精度とマッピング性が向上する（ref.17–19)。 Hybrid-Seqは、ゲノム構築とトランスクリプトームの特徴付けに適用でき、全体的な性能と分解能を向上させることができる（ref.11–13,17)。
PacBioとONTのリード長の長さは、トランスクリプトーム研究、特に遺伝子アイソフォームの同定にとって非常に有益である。ヒトのトランスクリプトーム（ref.20–22)
に加えて、PacBioトランスクリプトシークエンシングプロトコルであるIso-Seqが、非モデル生物および特定の遺伝子/遺伝子ファミリーにおけるトランスクリプトームの複雑さを特徴付けるために広く使用されている（ref.23–31)。対照的に、ONTには標準的な転写産物シークエンシングプロトコルはなく、公に利用可能なパイロット試験はわずかしかない。 Minionを使って、Bolisetty et alはショウジョウバエの4つの遺伝子の非常に高いアイソフォーム多様性を発見した。そして、それは複雑な転写イベントを調査することにおけるONTの有用性を例示する（ref.32)。 Oikonomopoulos et alは Spike-In RNAを用いた92の転写産物の人工混合物を分析することにより、トランスクリプトームの定量におけるONTシーケンシングの安定性も実証した（ref,33）。 PacBioまたはONTを単独で使用するこれらの研究と比較して、Hybrid-Seqは、特にトランスクリプトーム全体の研究において、データサイズの要件を軽減し、出力を向上させることができる。例えば、一連のHybrid-Seq法（IDP、IDP-fusion、IDP-ASE）は、トランスクリプトーム研究をisoformレベル（例えば、遺伝子isoform再構築、融合遺伝子および対立遺伝子フェージング）までより高い感度および精度で改善するために開発された。（一部略）
この論文では、hESCのcDNAからPacBioとONTのデータを生成した。我々（著者ら）のツールAlignQC（http://www.healthcare.uiowa.edu/labs/au/AlignQC/）を使用して、rawデータ（subreadsと1D “template” reads）、またはコンセンサス（CCSと2D reads）のPacBioとONTデータの包括的な分析と比較を行った 。製造業者はサイズ選択を推奨しているので、PacBioシーケンシングをサイズ選択ライブラリで行った。サイズ選択はシーケンシング時には標準的ではなく、そしてONTについては行われなかったので、ONTライブラリーはサイズ選択しなかった。これらの技術は異なるライブラリー調製プロトコルに従うため、シーケンシングプラットフォーム自体が変動性を導入し得るのと同様に、これらの工程を変動性の潜在的な原因として考慮することが重要である。比較には、エラー率とエラーパターン、リード長、マッピング可能性と異常なアライメント、そして最新のシーケンスモデル（PacBio P6-C4とONT R9）と以前のバージョン（C2とR7）の間の技術改善が含まれる。また、PacBioとONT単独の機能を検証して比較し、スパイクイントランスクリプトのゴールドスタンダードセットを研究した。その後、アイソフォームの同定、定量、複雑なトランスクリプトームイベントの発見など、ヒトのトランスクリプトーム解析にロングリードオンリーおよび対応するHybrid-Seqアプローチを適用した。 2つの主要なTGSデータプラットフォームの特性の包括的な評価に加え、この論文は、PacBioとONTの応用とそれに対応するトランスクリプトーム解析のためのHybrid-Seqのガイドとして役立つだろう。
ここではAlignQCの使い方についてのみ簡単に紹介します。
このツールはlong read RNA seq のbam分析用に設計されていますが、short readのRNA seq、DNAのシーケンスのbamでも動作します。
インストール
ubuntu16.04のpython2.7環境でテストした（docker使用。ホストOS macos10.14）。
依存
Report Generation Requirements
- R
- python 2.7+
Report Viewing Requirements
- A Modern Web Browser
本体 Github
#pipで導入可能
pip install AlignQC
#anaconda環境ならcondaで導入可能（テスト時は失敗した）
conda create -n alignqc -c vacation alignqc
source activate alignqc
> alignqc -h
# alignqc -h
usage: alignqc [-h] {analyze,dump,compare}
Version 2.0.5 Review reports about alignments.
positional arguments:
{analyze,dump,compare}
MODE of program to run
optional arguments:
-h, --help show this help message and exit
> alignqc analyze --help
# alignqc analyze --help
usage: alignqc [-h] (-g GENOME | --no_genome)
(-t GTF | --gpd GPD | --no_transcriptome) [-o OUTPUT]
[--portable_output PORTABLE_OUTPUT]
[--output_folder OUTPUT_FOLDER] [--threads THREADS]
[--tempdir TEMPDIR | --specific_tempdir SPECIFIC_TEMPDIR]
[--min_intron_size MIN_INTRON_SIZE]
[--min_aligned_bases MIN_ALIGNED_BASES]
[--max_query_overlap MAX_QUERY_OVERLAP]
[--max_target_overlap MAX_TARGET_OVERLAP]
[--max_query_gap MAX_QUERY_GAP]
[--max_target_gap MAX_TARGET_GAP]
[--required_fractional_improvement REQUIRED_FRACTIONAL_IMPROVEMENT]
[--do_loci] [--min_depth MIN_DEPTH]
[--min_coverage_at_depth MIN_COVERAGE_AT_DEPTH]
[--min_exon_count MIN_EXON_COUNT]
[--locus_downsample LOCUS_DOWNSAMPLE]
[--alignment_error_scale ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE]
[--alignment_error_max_length ALIGNMENT_ERROR_MAX_LENGTH]
[--context_error_scale CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE]
[--context_error_stopping_point CONTEXT_ERROR_STOPPING_POINT]
[--samples_per_xval SAMPLES_PER_XVAL]
[--max_bias_data MAX_BIAS_DATA] [--rscript_path RSCRIPT_PATH]
input
Create an output report
optional arguments:
-h, --help show this help message and exit
Input parameters:
Required BAM file. If reference or annotation is not set, use
--no_reference or --no_annotation respectively to continue.
input INPUT BAM file
-g GENOME, --genome GENOME
Reference Fasta (default: None)
--no_genome No Reference Fasta (default: False)
-t GTF, --gtf GTF Reference transcriptome in GTF format, assumes gene_id
and transcript_id are defining gene and transcript
names (default: None)
--gpd GPD Reference transcriptome in genePred format (default:
None)
--no_transcriptome No annotation is available (default: False)
Output parameters:
At least one output parameter must be set
-o OUTPUT, --output OUTPUT
OUTPUT xhtml with data (default: None)
--portable_output PORTABLE_OUTPUT
OUTPUT file in a small xhtml format (default: None)
--output_folder OUTPUT_FOLDER
OUTPUT folder of all data (default: None)
Performance parameters:
--threads THREADS INT number of threads to run. Default is system cpu
count (default: 1)
Temporary folder parameters:
--tempdir TEMPDIR The temporary directory is made and destroyed here.
(default: /tmp)
--specific_tempdir SPECIFIC_TEMPDIR
This temporary directory will be used, but will remain
after executing. (default: None)
Alignment plot parameters:
--min_intron_size MIN_INTRON_SIZE
minimum intron size when smoothing indels (default:
68)
--min_aligned_bases MIN_ALIGNED_BASES
for analysizing alignment, minimum bases to consider
(default: 50)
--max_query_overlap MAX_QUERY_OVERLAP
for testing gapped alignment advantage (default: 10)
--max_target_overlap MAX_TARGET_OVERLAP
for testing gapped alignment advantage (default: 10)
--max_query_gap MAX_QUERY_GAP
for testing gapped alignment advantge (default: None)
--max_target_gap MAX_TARGET_GAP
for testing gapped alignment advantage (default:
500000)
--required_fractional_improvement REQUIRED_FRACTIONAL_IMPROVEMENT
require gapped alignment to be this much better (in
alignment length) than single alignment to consider
it. (default: 0.2)
Locus parameters:
Optionally produce plots and data regarding clusters of sequences
--do_loci this analysis is time consuming at the moment
(default: False)
--min_depth MIN_DEPTH
require this or more read depth to consider locus
(default: 1.5)
--min_coverage_at_depth MIN_COVERAGE_AT_DEPTH
require at leas this much of the read be covered at
min_depth (default: 0.8)
--min_exon_count MIN_EXON_COUNT
Require at least this many exons in a read to consider
assignment to a locus (default: 2)
--locus_downsample LOCUS_DOWNSAMPLE
Limit how deep to search loci (default: 100)
Alignment error parameters:
--alignment_error_scale ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE ALIGNMENT_ERROR_SCALE
<ins_min> <ins_max> <mismatch_min> <mismatch_max>
<del_min> <del_max> (default: None)
--alignment_error_max_length ALIGNMENT_ERROR_MAX_LENGTH
The maximum number of alignment bases to calculate
error from (default: 1000000)
Context error parameters:
--context_error_scale CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE CONTEXT_ERROR_SCALE
<ins_min> <ins_max> <mismatch_min> <mismatch_max>
<del_min> <del_max> (default: None)
--context_error_stopping_point CONTEXT_ERROR_STOPPING_POINT
Sample at least this number of each context (default:
5000)
Rarefraction plot parameters:
--samples_per_xval SAMPLES_PER_XVAL
Bias parameters:
--max_bias_data MAX_BIAS_DATA
Bias does not need too large of a dataset. By default
data will be downsampled for large datasets. (default:
500000)
Path parameters:
--rscript_path RSCRIPT_PATH
The location of the Rscript executable. Default is
installed in path (default: Rscript)
> alignqc dump --help
# alignqc dump --help
usage: alignqc [-h] [-o OUTPUT] [-v] (-l | -e EXTRACT)
[--tempdir TEMPDIR | --specific_tempdir SPECIFIC_TEMPDIR]
xhtml_input
Extract data from xhtml output of alignqc through the command line
positional arguments:
optional arguments:
-h, --help show this help message and exit
-o OUTPUT, --output OUTPUT
OUTPUTFILE or STDOUT if not set (default: -)
-v, --verbose Show all stderr messages (default: False)
-l, --list show available data (default: False)
-e EXTRACT, --extract EXTRACT
dump this data (default: None)
--tempdir TEMPDIR The temporary directory is made and destroyed here.
(default: /tmp)
--specific_tempdir SPECIFIC_TEMPDIR
This temporary directory will be used, but will remain
after executing. (default: None)
> alignqc compare --help
# alignqc compare --help
usage: alignqc [-h] -o OUTPUT [--threads THREADS]
[--tempdir TEMPDIR | --specific_tempdir SPECIFIC_TEMPDIR]
xhtml_inputs [xhtml_inputs ...]
positional arguments:
xhtml_inputs xhtml analysis files
optional arguments:
-h, --help show this help message and exit
-o OUTPUT, --output OUTPUT
OUTPUT directory (default: None)
--threads THREADS INT number of threads to run. Default is system cpu
count (default: 18)
--tempdir TEMPDIR The temporary directory is made and destroyed here.
(default: /tmp)
--specific_tempdir SPECIFIC_TEMPDIR
This temporary directory will be used, but will remain
after executing. (default: None)
実行方法
analyze
リファレンスゲノムがある場合
alignqc analyze long_reads.bam -g reference.fa.gz \
-t ref_transcriptome.gtf.gz -o output.xhtml
- --threads INT number of threads to run. Default is system cpu count (default: 1)
- -t Reference transcriptome in GTF format, assumes gene_id and transcript_id are defining gene and transcript names (default: None)
- -g Reference Fasta (default: None)
- -o OUTPUT xhtml with data (default: None)
- --portable_output OUTPUT file in a small xhtml format (default: None)
- --output_folder OUTPUT folder of all data (default: None)
出力フォーマットはxhtmlになっているが、拡張子をhtmlにして出力しても開ける。
--output_folderで出力するとディレクトリに関連ファイルが全て出力される。
xhtmlファイルを開く。
リファレンスゲノムがない場合
alignqc analysis long_reads.bam --no_genome --no_transcriptome \
-o output.html
- --no_transcriptome No annotation is available (default: False)
- --no_genome No Reference Fasta (default: False)
compare
複数結果を比較するstatisticsテーブルを出力する。
alignqc compare -o output_dir --threads 16 sample1/output_report.xhtml sample2/output_report.xhtml <sample3 sample4 ...>
出力txtをexcelで開いた。
引用
Comprehensive comparison of Pacific Biosciences and Oxford Nanopore Technologies and their applications to transcriptome analysis
Jason L Weirather, Mariateresa de Cesare, Yunhao Wang x, Paolo Piazza, Vittorio Sebastiano, Xiu-Jie Wang, David Buck, Kin Fai Au
F1000 Research
2019年４月現在、査読中になっています。現在version2です。
*1
ショートリードではそこまででもないが、ラージゲノムのlong_reads.bamだとかなりメモリを使うためkillしてしまう。私の環境ではchr１つだけでもメモリが足りず落ちてしまったため、Mtゲノムだけ分析しています。注意して下さい。
関連
