macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

long readのRNA seq向けマッピング評価ツール AlignQC

 

 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)]は、1DNA分子から非常に長いヌクレオチド配列を直接捕捉するための新しい生化学/生物物理学的方法を開発した。他の新たな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データの分析に関連する制限に対処するためのアプローチとして浮上している。たとえば、Pac​​Sioまたは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

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:

  xhtml_input           INPUT XHTML FILE

 

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で出力するとディレクトリに関連ファイルが全て出力される。

f:id:kazumaxneo:20190422130950j:plain

xhtmlファイルを開く。

f:id:kazumaxneo:20190422115541j:plain

f:id:kazumaxneo:20190422115557j:plain

f:id:kazumaxneo:20190422115600j:plain

 

 

f:id:kazumaxneo:20190422115645j:plain

 

 リファレンスゲノムがない場合

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で開いた。

f:id:kazumaxneo:20190422133128j:plain

 

引用

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年4月現在、査読中になっています。現在version2です。

 

*1

ショートリードではそこまででもないが、ラージゲノムのlong_reads.bamだとかなりメモリを使うためkillしてしまう。私の環境ではchr1つだけでもメモリが足りず落ちてしまったため、Mtゲノムだけ分析しています。注意して下さい。

 

関連