macでインフォマティクス

macでインフォマティクス

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

高速かつ高感度なRNA/DNAのアライナー HPG Aligner

 

 ハイスループットシーケンサーの最新世代は、前例のないスケールでデータを生成し、関連するシーケンシングコストが連続的に減少している。特に、トランスクリプトームの包括的なプロファイルを提供するRNAシーケンシング(RNA-seq)技術(論文より ref.1)は、従来のマイクロアレイをますます置き換えている(ref.2)。RNA-seqにおける一次データプロセシング(ゲノムリシークエンシングを含む他の大規模シーケンシング実験)は、リファレンスゲノムへのリードのマッピングを含む。この段階は、計算上高価なプロセスであり、加えて感度が重大な関心事である(ref.3)。さまざまなプログラムが存在し、その多くは、検索スピードを大幅に高速化させるindex技術であるBWT(Burrows-Wheeler Transform)アルゴリズムを実装している。

 近年提案されているマッピングアルゴリズムは、DNAリードの約98〜99%をマップする(これまでの利用可能なマッパーで得られた96%に対してわずか2%のゲイン)。 しかし、RNA-seqのシナリオはこれらの基準からは程遠い。転写物のシーケンシングリードのマッピングは、単純にゲノム配列にリードをマッピングするよりもはるかに複雑な問題である。真核生物のトランスクリプトームは複雑で、平均で1遺伝子あたり9つの転写物をもち(ref.12)、250bpのエキソンを遺伝子あたり12個以上持ち、その長さは何十万bpにもおよぶ。さらに新しい転写物アイソフォームが絶えず発見されている(ref.13)。シーケンシングエラーを伴う現実的なシナリオでは、 TopHatなど(ref.14)は、100 bpのリードの70%しか正確にマッピングせず、ようやくGSNAPによって最高の感度(82%)が達成されている(ref.15) 。この新しい報告によると、RUM(ref.17)、STAR(ref.18)またはMapSplice(ref.19)のような他のRNA-seqマッパーは、シーケンシングリード精度がさらに劣る。しかし、より最近の報告から(ref.20)、STARやMapSpliceのパフォーマンスは他の方法と比較してより良好で、これは以前のベンチマークと一致している(ref.17)。リード長が長くなると(シーケンシングマシンの連続的なアップグレードの自然な傾向)感度が低下し、感度とランタイムの点で重大な問題を引き起こす。ほとんどのBWTベースのマッパーでは問題がある。なぜなら、BWTベースのマッパーは少数のミスマッチしか許さないためである(変化が少ないシナリオで考えても現在のリード長では不十分) 。Bower2(ref.10)を使用するTopHat2(ref.14)のような新しいアライナーはマッピングステップを大幅に加速したが(感度のトレードオフが小さい)、ジャンクション検出の戦略は変わらない。この研究が検討されている間、BWTとFerragina-Manzini indexをスキームに組み合わせた新しいマッパーHISAT (ref.21)がpublishされた。他と同様に、スピードアップは感度とのトレードオフで達成されたようである。

 トランスクリプトーム解析のためにRNA-seqを使う関心の高さや定量のためのアルゴリズムの進歩、transcripts発見のためのアルゴリズムの進歩にも関わらず、基本的な問題であるリードのマッピング問題が解決されていない。

 ここでは、BWTによるマッピングとSmith-Waterman(SW)とのローカルアラインメントの組み合わせに基づいて、ショートリードとロングリード両方のハイクオリティなマッピングのための革新的なソリューションを提案する。これは、マッピング精度を大幅に高め(典型的なシナリオで現在のマッパーが60 %に対して95 %)、ランタイムを大幅に短縮する(TopHat2と比べて約15倍、MapSpliceで8倍、STARで2倍)。さらに、提案された戦略は、indelとミスマッチに対して非常に堅牢であることも示されている。この提案は、多数のミスマッチやindelを含むものであっても、ほぼすべてのリードをマッピングするシンプルで高速で洗練されたソリューションを提供する。この解決策はマッピングステップにおいてかなりの時間を節約し、その結果、シーケンスデータ処理の現在のパイプライン高速化に決定的に寄与する。さまざまなハイパフォーマンスコンピューティング (HPC) テクノロジを組み合わせたHPG Alignerは、ショートリードとロングリード両方で優れたパフォーマンスを示す。

 

チュートリアル

https://github.com/opencb/hpg-aligner/wiki/Tutorial

 

 

インストール

cent os6でテストした。

Github

https://github.com/opencb/hpg-aligner

yum install gcc 
yum install scons
yum install zlib-devel libcurl-devel libxml2-devel ncurses-devel gsl-devel check-devel


git clone https://github.com/opencb/hpg-aligner.git
cd hpg-aligner
git submodule update --init
cd lib/hpg-libs
git checkout master
cd ../..
git checkout master
scons
cd bin/

./hpg-aligner -h

$ ./hpg-aligner -h

Usage:

hpg-aligner {dna | rna | build-bwt-index | build-sa-index} [-f|--fq|--fastq=<file>] [-j|--fq2|--fastq2=<file>] [-z|--gzip] [-i|--index=<file>] [-o|--outdir=<file>] [--filter-read-mappings=<int>] [--filter-seed-mappings=<int>] [--min-cal-size=<int>] [-t|--cpu-threads=<int>] [--read-batch-size=<int>] [--sw-match=<double>] [--sw-mismatch=<double>] [--sw-gap-open=<double>] [--sw-gap-extend=<double>] [--min-score=<int>] [--paired-min-distance=<int>] [--paired-max-distance=<int>] [--report-best] [--report-all] [--report-n-best=<int>] [--report-n-hits=<int>] [--report-only-paired] [--prefix=<string>] [-l|--log-level=<int>] [-h|--help] [--bam-format] [--indel-realignment] [--recalibration] [-a|--adapter=<string>] [-v|--version] [--num-seeds=<int>]

-f, --fq, --fastq=<file>                          Reads file input. For more than one file: f1.fq,f2.fq,...

-j, --fq2, --fastq2=<file>                        Reads file input #2 (for paired mode)

-z, --gzip                                        FastQ input files are gzip

-i, --index=<file>                                Index directory name

-o, --outdir=<file>                               Output directory

--filter-read-mappings=<int>                      Reads that map in more than <n> locations are discarded

--filter-seed-mappings=<int>                      Seeds that map in more than <n> locations are discarded

--min-cal-size=<int>                              Minimum CAL size

-t, --cpu-threads=<int>                           Number of CPU Threads

--read-batch-size=<int>                           Batch Size

--sw-match=<double>                               Match value for Smith-Waterman algorithm

--sw-mismatch=<double>                            Mismatch value for Smith-Waterman algorithm

--sw-gap-open=<double>                            Gap open penalty for Smith-Waterman algorithm

--sw-gap-extend=<double>                          Gap extend penalty for Smith-Waterman algorithm

--min-score=<int>                                 Minimum score for valid mappings (0 to 100 for RNA)

--paired-min-distance=<int>                       Minimum distance between pairs

--paired-max-distance=<int>                       Maximum distance between pairs

--report-best                                     Report all alignments with best score

--report-all                                      Report all alignments

--report-n-best=<int>                             Report the <n> best alignments

--report-n-hits=<int>                             Report <n> hits

--report-only-paired                              Report only the paired reads

--prefix=<string>                                 File prefix name

-l, --log-level=<int>                             Log debug level

-h, --help                                        Help option

--bam-format                                      BAM output format (otherwise, SAM format. This option is only available for SA mode, BWT mode always report in BAM format), this option turn the process slow

--indel-realignment                               Indel-based realignment

--recalibration                                   Base quality score recalibration

-a, --adapter=<string>                            Adapter sequence in the read

-v, --version                                     Display the HPG Aligner version

--num-seeds=<int>                                 Number of seeds

hpg-aligner build-bwt-index -h

$ hpg-aligner build-bwt-index -h

Usage:

hpg-aligner {build-sa-index | build-bwt-index} [-i|--index=<file>] [-g|--ref-genome=<file>] [-h|--help] [-v|--version] [-r|--index-ratio=<int>]

-i, --index=<file>                                Index directory name

-g, --ref-genome=<file>                           Reference genome

-h, --help                                        Help option

-v, --version                                     Display HPG Aligner version

-r, --index-ratio=<int>                           BWT index compression ratio

 

ラン

 RNAマッピング

1、index作成。Burrows-Wheeler変換してメモリー使用量を抑える。

mkdir index_dir #indexディレクトリの作成
hpg-aligner build-bwt-index -g ref.fasta -i index_dir -r 8

 

2、マッピング。デフォルトでは全スレッドが使用される。

hpg-aligner rna -i index_dir -f pair1.fq -j pair2.fq -o output-dir

 output-dir/にbam出力される。

 

 

パフォーマンス比較。RNAはページの一番下にあります。

https://github.com/opencb/hpg-aligner/wiki/Results-and-benchmarks

 

引用

Highly sensitive and ultrafast read mapping for RNA-seq analysis

I. Medina, J. Tárraga, H. Martínez, S. Barrachina, M. I. Castillo, J. Paschall, J. Salavert-Torres, I. Blanquer-Espert, V. Hernández-García, E. S. Quintana-Ortí, and J. Dopazo

DNA Res. 2016 Apr; 23(2): 93–100.

 

並列化に対応した高感度なアダプタートリミングツール PEAT

 

 次世代シークエンシング(NGS)プラットフォームでよく知られているシングルエンドシーケンシング技術からmodifyされたペアエンドシーケンシング技術は、ゲノミクスにおいてますます重要な役割を果たしている。 DNA(またはcDNA)断片の2つの鎖の5 '末端を配列決定することにより、DNA contentsだけでなく位置情報も提供するので、反復領域おけるアセンブリまたは構造変異を解決する強力なリソースである。さらに、それは、シーケンシングされたDNAの両方の鎖に関する情報を伝達し、エキソンジャンクションを分析し、他の多くのアプリケーションに応用できる[論文より ref.3]。

 典型的なイルミナのペアエンドシーケンシング技術は以下のように実施される:DNA断片の二本鎖はアダプターとバーコードの両方に連結され(multiplexingの場合)、二本鎖の5 '末端がフローセル表面に結合され、より良いヌクレオチド合成および蛍光イメージングのためのクラスターを生成する。各DNA断片の両方の鎖は、クラスターを再生することによってシーケンシングテンプレートとして機能することができ、その結果、対の末端のリードが生成される(論文 図1A)。イルミナのペアエンドシーケンスプロトコルによれば、ユーザーは200〜500塩基対(bp)の選択した長さのDNA断片を2億リード得ることができる。得られたペアエンドリードは、36 x 2,75 x 2、または100 x2 bpなどのマシン固有の配列長である。

 少なくとも予め指定された配列長(例えば100bp)と同じ長さのDNA断片については、両方の鎖の5 '末端から開始されるシーケンスプロセスは、アダプター配列前に終了するので、結果として実際のDNA情報のみ得ることができる。しかし、DNA断片が100bpよりも短い場合、シーケンサーは実際のDNAをリードスルーしてアダプターまでシーケンスする(論文 図1B)。その結果、生成されたペアリードには、不要なアダプタシーケンス(アダプタの汚染)が追加され、リファレンスへのマッピングのステップで排除される可能性がある。

 リードから「実際の」DNA部分を回復するためには、既存のアルゴリズムはリード3 '末端に対してアダプター配列間の正しいアライメントに依存する。このようなアルゴリズムは、もともとシングルエンドシーケンシングにおけるアダプター汚染を扱うために設計されたものである。ペアレントエンドリードは2セットのシングルエンドリードとみなされ、アダプターのトリミングは各セットごとに個別に実行される。実際には、FASTX [ref.4]、Cutadapt [ref.5]、ea-utils [ref.6]、TagCleaner [ref.7]、Trim_Galore [ref.8](ペアエンドシーケンシング用cutadapterのラッパー)、ペアエンドシーケンシングでは満足のいく結果が得られない。これは、今日のシーケンサーがリード3 '末端において有意に高い配列エラー率を有するためである(図1C)。言い換えれば、不十分なシーケンシング品質を有するリードセグメントに対してアダプター配列を適合させることによって推定トリミング位置を同定することは、誤りが起こりやすい。 AdapterRemoval [ref.9](紹介)、SeqPrep [ref.10]、Trimmomatic [ref.11]のような一部のトリマーは、ペアリードのトリムされていない部分の間の逆相補性を調べることによってアダプターを探すが、その精度は依然としてエラーの起こりやすい3 '末端でのアダプター配列スキャニングに依存する。 GATK ReadAdaptorTrimmer [ref.12]は、高精度を達成するために、ペアのforwardリードとreverseリードの重なりを探す。しかし、GATK [ref.13]のバリアントコール解析パイプラインのために特別に設計されているため、入力と出力はSAM / BAM形式[ref.14]でなければならず、他の下流分析には煩雑である。いくつかのトリマーの性能比較は、いくつかのベンチマーク研究によって調査されている[ref.15,16,17]。

 ここでは、ペアワイズシーケンシングのために特別に設計されたPEAT(Paired-End Adapter Trimmer)という、効率的で正確なアダプタートリミングアルゴリズムとその実装を提案する。 PEATはアダプター配列の入力を必要としない。これは、大規模な異なるアダプターを使用してライブラリーを処理する場合に特に便利である。 PEATは、従来のアプローチによって採用されたフィルタリング中の感度の損失を回避するために、リードのハイクオリティ部分の逆相補を直接スキャンする。 PEATと多くのアダプタトリミングツールを比較した。 PEATは、シミュレートされたベンチマークで比較的上手く動き、大きなリアルデータセットに適用された場合は高いスケーラビリティを示した。著者らはPEATを2つの公開リアルデータセット(1億5千万のシークエンシングリードを含む101 x 2のシーケンシングライブラリ)に適用し、何百万のアダプターつきリードが検出され、リカバリされ、リファレンスにバックマッピングされた。さらに、ChIP-seq、MNase-seq、およびRNA-seqのような異なる10種の公に利用可能なデータセットに対するPEATの影響を調べた。 下流分析をペアTのの処理の有無で比較すると、明らかなパターン変化が明らかになり、これはデータを偏向させ生物学的な意味を変える可能性がある。著者らは、すべてのアプリケーションでペアエンドリードを分析する際にアダプタの汚染にもっと注意を払う必要があることを提案している。

 

インストール

mac os 10.12でテストした。

Github

https://github.com/jhhung/PEAT

Githubリリースよりバイナリをダウンロードできる。

https://github.com/jhhung/PEAT/releases

chmod u+x PEAT

> PEAT 

$./PEAT 

 

 

*********************************************************************************

+----+

|PEAT|

+----+

  A integrated software that can do either paired-end and single-end adapter trimming operation.  

 

Usage:

 

 1. paired-end adapter trimming

> PEAT paired --help

  

 2. single-end adapter trimming

> PEAT single --help

*********************************************************************************

 

> PEAT paired --help 

$ PEAT paired --help

Error: the option '--input1' is required but missing

 

*********************************************************************************

+----------+

|Paired End|

+----------+

 

A software do paired-end adapter trimming operation.

It takes paired-end FastQ format input files (dual files), and reports adapter 

removed FastQ format output files (dual files).

 

>> PEAT paired

    Do paired-end adapter trimming operation with instruction like:

    bin/PEAT_linux paired -1 test_file/test_paired1.fq -2 test_file/test_paired2.fq 

 

*********************************************************************************

:

  -h [ --help ]                        display this help message and exit

  -1 [ --input1 ] arg                  The paired_1 input FastQ file (.fq) or 

                                       Gzip compressed FASTQ file (.fq.gz).

  -2 [ --input2 ] arg                  The paired_2 input FastQ file (.fq) or 

                                       Gzip compressed FASTQ file (.fq.gz).

  -o [ --output ] arg (=stdout)        Prefix for Output file name, stdout by 

                                       default. If you choose this option, you 

                                       couldn't use --output_1 and --output_2

  --output_1 arg (=stdout)             Prefix for Output file part1 name, 

                                       stdout by default 

  --output_2 arg (=stdout)             Prefix for Output file part2 name, 

                                       stdout by default 

  -n [ --thread ] arg (=1)             Number of thread to use; if the number 

                                       is larger than the core available, it 

                                       will be adjusted automatically

  -l [ --len ] arg (=30)               Minimum gene fragment length, i.e. the 

                                       fragment length for reverse complement 

                                       check, 30 bp by default

  -r [ --reverse_mis_rate ] arg (=0.4) Mismatch rate applied in first stage 

                                       reverse complement scan, 0.4 by default

  -g [ --gene_mis_rate ] arg (=0.6)    Mismatch rate applied in second stage 

                                       gene portion check, 0.6 by default

  -a [ --adapter_mis_rate ] arg (=0.4) Mismatch rate applied in second stage 

                                       adapter portion check, 0.4 by default

  --qtrim                              Quality trimmer; trim the last base of 

                                       the reads until the mean quality value 

                                       of the reads is larger than threshold

  -q [ --quality ] arg                 The quality type. Type any one of the 

                                       following quality type indicator: 

                                       ILLUMINA, PHRED, SANGER, SOLEXA

                                       Only for the option: --qtrim

  -t [ --threshold ] arg               The threshold (quality value) of the 

                                       quality trimmer, 30.0 by default

                                       Only for the option: --qtrim

  --out_gzip                           Compress the FASTQ output to Gzip file. 

                                       This option is required the option: -o 

                                       or --output_1/--output_2

  --verbose                            Output running process by stderr

  --adapter_contexts                   Output adapter contexts within the top 

                                       ten numbers in report.txt; You can use 

                                       the option: adapter_min_bp to select 

                                       adapter what you want; if you use this 

                                       option, the program becomes slower.

  --adapter_min_bp arg (=10)           Determine the mininal length of output 

                                       adapter contexts within the top 50 

                                       numbers in report.txt, 10 bp by default;

                                       Required the option: adapter_contexts

 

 

ラン

ペアエンドのアダプタートリミング。アダプターは自動で検出される。

PEAT paired -1 pair1.fq -2 pair2.fq -o out -n 4 --adapter_contexts
  • -1    The paired_1 input FastQ file (.fq) or Gzip compressed FASTQ file (.fq.gz).
  • -2    The paired_2 input FastQ file (.fq) or Gzip compressed FASTQ file (.fq.gz).
  • -o    Prefix for Output file name, stdout by  default. 
  • -n     Number of thread to use; if the number is larger than the core available, it will be adjusted automatically
  • --adapter_contexts    Output adapter contexts within the top ten numbers in report.txt; You can use the option: adapter_min_bp to select adapter what you want; if you use this option, the program becomes slower.

out_paired1.fqとout_paired2.fq、top10sレポートが出力される。

 

ライセンスはGPLv2になっています。

引用

PEAT: an intelligent and efficient paired-end sequencing adapter trimming algorithm.

Li YL, Weng JC, Hsiao CC, Chou MT, Tseng CW, Hung JH 

BMC Bioinformatics. 2015;16 Suppl 1:S2.

 

並列化に対応したアダプタートリミングツール AdapterRemoval 2

 

 Fossil material 由来などのごく短いDNA断片のハイスループットシーケンスでは、ライブラリーの調製中にインサートにライゲーションされたアダプター配列をシークエンシングする可能性がある[論文より ref.1]。このような汚染はよく知られた問題であり、下流分析にマイナスの影響を与える可能性がある[ref.2-5]。したがって、ワークフローの最初の部分には、通常、アダプタの汚染をフィルタリングまたは除去(トリミング)するステップが含まれている[ref.3]。(古代のクオリティの低いサンプルでも)ペアエンドリードののオーバーラップする領域を検出し、これらを品質認識的に崩壊(合併)させて鋳型分子全体を再構築することにより、精度をあげることができる[ref.6]。これは、短いインサートサイズが期待されるancient DNAシークエンシング、およびリード末端への死後のDNA修飾に起因する、すでに高い配列決定エラー率の減少が興味深い[ref.7 pubmed]場合に特に興味深い。

 AdapterRemoval v1 [ref.2]のオリジナルリリースでは、Needleman-Wunschアルゴリズムの修正版を使用して、アダプター配列と低品質塩基のトリミングのための使いやすいツールが提供されており、リード間のペアワイズアライメント(またはペア)および既知のアダプター配列を含まれる。シングルエンドリードの場合、このアライメントは、アダプター配列の5 '末端とrawシーケンシングリードの3'末端との間の最良のアラインメントを見出し、ペアエンドリードの場合には、mate1のリードと、mate2のリードの逆相補との間で最良のアラインメントを見出す。

 これらの凝集配列のアライメントは、各リードにおけるの3'末端の同定を可能にし、それにより外来性(アダプター)配列をトリミングすることを可能にする[ ref.2]。しかしながら、使用された改良Needleman-Wunschアルゴリズムはシーケンシング技術によって導入されたindelを考えないので、AdapterRemovalは主に、false indelsが少ないIllumina HTSプラットフォームから生成されたリード処理に適している。  AdapterRemovalはさらに、(オプションとして)重複するリードをマージするため、アダプターの位置合わせ手順の一部として重複フラグメントを使用する。これは、最高品質の塩基を選択し、重複位置で位置特異的スコアリングマトリックスとして扱うことで基本品質を再計算することによって達成される[ref.2]。

 しかし、AdapterRemoval v1は他の最新のツール[ref.8]と比較して、実行時間が比較的遅いという特徴がある。したがって、採用されたトリミング方法を変更することなく、スループットを向上させる目的で、AdapterRemoval v1の大幅な改訂を実施した。したがって、更新されたバージョンは、AdapterRemoval v1と同様の精度を示し、既存の分析パイプライン[ref.7]での使用に適切な互換品を提供する。

 AdapterRemoval v2は、AdapterRemoval v1で使用されるアライメントアルゴリズムを高速化するために、単一命令、複数データ(SIMD)命令(つまりコンシューマーグレードCPUで一般的にサポートされるSSE命令とSSE2命令)を使用することにより、複数のスレッドを使用してすべての操作のスループットをさらに向上させることができる。さらに、AdapterRemoval v2では、複数の異なる(ペアの)アダプタシーケンスを同時にトリミングし、リード(ペア)ごとに最適な一致を選択し、gzipおよびbzip2圧縮FASTQファイルを一貫して読み書きできる。 AdapterRemoval v2では、さらに、提供されたバーコードの単純な最大不一致数の比較を使用して、同時にdemultiplexingとアダプターのトリミングが可能になり、インサートの端を超えてシーケンスされるペアエンドリードも正しくトリミングされる。このようなリードは、ペアエンドの相方の逆相補バーコード配列に続いてアダプター配列によって終結され、両方とも除去されなければならない。最後に、AdapterRemoval v2は、インサートがオーバーラップするペアから推定アダプター配列を再構成することができる。これを使用して、実験エラーを検出し、文書化されていないデータセットを分析することができる。

 AdapterRemoval v2のパフォーマンスを評価するために、アダプターシーケンスを正しく調整する機能と、1つまたは複数のスレッドを使用する場合のスループットの点で、AdapterRemoval v1と現代のアダプタートリミングソフトウェアの選択肢を比較した。さらに、いくつかの他の重複シーケンスをマージするツールを使いAdapterRemovalの機能を比較した。後者については、比較をアダプター汚染の存在下でリードをマージするように設計されたプログラムに制限したが、アダプター汚染をほとんどまたはまったく含まないデータセットに適したいくつかの選択肢、例えばCOPE [9]、fastq-join [10]、FLASH [11]、およびXORRO(arXiv:1304.4620)に注意する。

 

 

 特徴(Githubより)

  • Trimming of adapters sequences from single-end and paired-end FASTQ reads.
  • Trimming of multiple, different adapters or adapter pairs.
  • Demultiplexing of single or double indexed reads, with or without trimming of adapter sequences.
  • Reconstruction of adapter sequences from paired-end reads, by the pairwise alignment of reads in the absence of a known adapter sequence.
  • Merging of overlapping read-pairs into higher-quality consensus sequences.
  • Multi-threading of all operations for increased throughput.
  • Reading and writing of gzip and bzip2 compressed files.
  • Reading and writing of interleaved FASTQ files.

 

インストール

mac os10.12でテストした。

本体 Github

https://github.com/MikkelSchubert/adapterremoval

git clone https://github.com/MikkelSchubert/adapterremoval.git 
cd adapterremoval
make
sudo make install

$ sudo make install

Password:

Building AdapterRemoval with gzip support: yes

Building AdapterRemoval with bzip2 support: yes

Building AdapterRemoval with pthreads support: yes

Linking executable 'build/AdapterRemoval'

clang: warning: argument unused during compilation: '-ansi' [-Wunused-command-line-argument]

Constructing man-page 'build/AdapterRemoval.1' from 'AdapterRemoval.pod'

Installing AdapterRemoval ..

  .. binary into /usr/local/bin/

  .. man-page into /usr/local/share/man/man1/

  .. README into /usr/local/share/adapterremoval/

  .. examples into /usr/local/share/adapterremoval/examples/

 

AdapterRemoval

$ AdapterRemoval

AdapterRemoval ver. 2.1.7

 

This program searches for and removes remnant adapter sequences from

your read data.  The program can analyze both single end and paired end

data.  For detailed explanation of the parameters, please refer to the

man page.  For comments, suggestions  and feedback please contact Stinus

Lindgreen (stinus@binf.ku.dk) and Mikkel Schubert (MikkelSch@gmail.com).

 

If you use the program, please cite the paper:

    Schubert, Lindgreen, and Orlando (2016). AdapterRemoval v2: rapid

    adapter trimming, identification, and read merging.

    BMC Research Notes, 12;9(1):88.

 

    http://bmcresnotes.biomedcentral.com/articles/10.1186/s13104-016-1900-2

 

 

Arguments:                           Description:

  --help                             Display this message.

  --version                          Print the version string.

 

  --file1 FILE                       Input file containing mate 1 reads or single-ended reads [REQUIRED].

  --file2 FILE                       Input file containing mate 2 reads [OPTIONAL].

 

FASTQ OPTIONS:

  --qualitybase BASE                 Quality base used to encode Phred scores in input; either 33, 64, or solexa

                                       [current: 33].

  --qualitybase-output BASE          Quality base used to encode Phred scores in output; either 33, 64, or solexa. By

                                       default, reads will be written in the same format as the that specified using

                                       --qualitybase.

  --qualitymax BASE                  Specifies the maximum Phred score expected in input files, and used when writing

                                       output. ASCII encoded values are limited to the characters '!' (ASCII = 33) to

                                       '~' (ASCII = 126), meaning that possible scores are 0 - 93 with offset 33, and

                                       0 - 62 for offset 64 and Solexa scores [default: 41].

  --mate-separator CHAR              Character separating the mate number (1 or 2) from the read name in FASTQ

                                       records [default: '/'].

  --interleaved                      This option enables both the --interleaved-input option and the

                                       --interleaved-output option [current: off].

  --interleaved-input                The (single) input file provided contains both the mate 1 and mate 2 reads, one

                                       pair after the other, with one mate 1 reads followed by one mate 2 read. This

                                       option is implied by the --interleaved option [current: off].

  --interleaved-output               If set, trimmed paired-end reads are written to a single file containing mate 1

                                       and mate 2 reads, one pair after the other. This option is implied by the

                                       --interleaved option [current: off].

 

OUTPUT FILES:

  --basename BASENAME                Default prefix for all output files for which no filename was explicitly set

                                       [current: your_output].

  --settings FILE                    Output file containing information on the parameters used in the run as well as

                                       overall statistics on the reads after trimming [default: BASENAME.settings]

  --output1 FILE                     Output file containing trimmed mate1 reads [default: BASENAME.pair1.truncated

                                       (PE), BASENAME.truncated (SE), or BASENAME.paired.truncated (interleaved PE)]

  --output2 FILE                     Output file containing trimmed mate 2 reads [default: BASENAME.pair2.truncated

                                       (only used in PE mode, but not if --interleaved-output is enabled)]

  --singleton FILE                   Output file to which containing paired reads for which the mate has been

                                       discarded [default: BASENAME.singleton.truncated]

  --outputcollapsed FILE             If --collapsed is set, contains overlapping mate-pairs which have been merged

                                       into a single read (PE mode) or reads for which the adapter was identified by

                                       a minimum overlap, indicating that the entire template molecule is present.

                                       This does not include which have subsequently been trimmed due to low-quality

                                       or ambiguous nucleotides [default: BASENAME.collapsed]

  --outputcollapsedtruncated FILE    Collapsed reads (see --outputcollapsed) which were trimmed due the presence of

                                       low-quality or ambiguous nucleotides [default: BASENAME.collapsed.truncated]

  --discarded FILE                   Contains reads discarded due to the --minlength, --maxlength or --maxns options

                                       [default: BASENAME.discarded]

 

OUTPUT COMPRESSION:

  --gzip                             Enable gzip compression [current: off]

  --gzip-level LEVEL                 Compression level, 0 - 9 [current: 6]

  --bzip2                            Enable bzip2 compression [current: off]

  --bzip2-level LEVEL                Compression level, 0 - 9 [current: 9]

 

TRIMMING SETTINGS:

  --adapter1 SEQUENCE                Adapter sequence expected to be found in mate 1 reads [current:

                                       AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG].

  --adapter2 SEQUENCE                Adapter sequence expected to be found in mate 2 reads [current:

                                       AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT].

  --adapter-list FILENAME            Read table of white-space separated adapters pairs, used as if the first column

                                       was supplied to --adapter1, and the second column was supplied to --adapter2;

                                       only the first adapter in each pair is required SE trimming mode [current:

                                       <not set>].

 

  --mm MISMATCH_RATE                 Max error-rate when aligning reads and/or adapters. If > 1, the max error-rate

                                       is set to 1 / MISMATCH_RATE; if < 0, the defaults are used, otherwise the

                                       user-supplied value is used directly. [defaults: 1/3 for trimming; 1/10 when

                                       identifing adapters].

  --maxns MAX                        Reads containing more ambiguous bases (N) than this number after trimming are

                                       discarded [current: 1000].

  --shift N                          Consider alignments where up to N nucleotides are missing from the 5' termini

                                       [current: 2].

 

  --trimns                           If set, trim ambiguous bases (N) at 5'/3' termini [current: off]

  --trimqualities                    If set, trim bases at 5'/3' termini with quality scores <= to --minquality value

                                       [current: off]

  --minquality PHRED                 Inclusive minimum; see --trimqualities for details [current: 2]

  --minlength LENGTH                 Reads shorter than this length are discarded following trimming [current: 15].

  --maxlength LENGTH                 Reads longer than this length are discarded following trimming [current:

                                       4294967295].

  --collapse                         When set, paired ended read alignments of --minalignmentlength or more bases are

                                       combined into a single consensus sequence, representing the complete insert,

                                       and written to either basename.collapsed or basename.collapsed.truncated (if

                                       trimmed due to low-quality bases following collapse); for single-ended reads,

                                       putative complete inserts are identified as having at least

                                       --minalignmentlength bases overlap with the adapter sequence, and are written

                                       to the the same files [current: off].

  --minalignmentlength LENGTH        If --collapse is set, paired reads must overlap at least this number of bases to

                                       be collapsed, and single-ended reads must overlap at least this number of

                                       bases with the adapter to be considered complete template molecules [current:

                                       11].

  --minadapteroverlap LENGTH         In single-end mode, reads are only trimmed if the overlap between read and the

                                       adapter is at least X bases long, not counting ambiguous nucleotides (N); this

                                       is independant of the --minalignmentlength when using --collapse, allowing a

                                       conservative selection of putative complete inserts while ensuring that all

                                       possible adapter contamination is trimmed [current: 0].

 

DEMULTIPLEXING:

  --barcode-list FILENAME            List of barcodes or barcode pairs for single or double-indexed demultiplexing.

                                       Note that both indexes should be specified for both single-end and paired-end

                                       trimming, if double-indexed multiplexing was used, in order to ensure that the

                                       demultiplexed reads can be trimmed correctly [current: <not set>].

  --barcode-mm N                     Maximum number of mismatches allowed when counting mismatches in both the mate 1

                                       and the mate 2 barcode for paired reads.

  --barcode-mm-r1 N                  Maximum number of mismatches allowed for the mate 1 barcode; if not set, this

                                       value is equal to the '--barcode-mm' value; cannot be higher than the

                                       '--barcode-mm value'.

  --barcode-mm-r2 N                  Maximum number of mismatches allowed for the mate 2 barcode; if not set, this

                                       value is equal to the '--barcode-mm' value; cannot be higher than the

                                       '--barcode-mm value'.

 

MISC:

  --identify-adapters                Attempt to identify the adapter pair of PE reads, by searching for overlapping

                                       reads [current: off].

  --seed SEED                        Sets the RNG seed used when choosing between bases with equal Phred scores when

                                       collapsing. Note that runs are not deterministic if more than one thread is

                                       used. If not specified, a seed is generated using the current time.

  --threads THREADS                  Maximum number of threads [current: 1]

 

——

 

condaでも導入できます(Github)。

 

 

ラン

/usr/local/share/adapterremoval/examples/にあるテストシーケンスデータをランする。リンクを貼っておく。

ln -s /usr/local/share/adapterremoval/examples/* .

 

 

1、アダプター配列を探す。

AdapterRemoval --identify-adapters --file1 reads_1.fq --file2 reads_2.fq 
  • --identify-adapters    Attempt to identify the adapter pair of PE reads, by searching for overlapping  reads [current: off].
  • --file1    Input file containing mate 1 reads or single-ended reads [REQUIRED].  
  • --file2    Input file containing mate 2 reads [OPTIONAL].

 

2、シングルエンドのアダプタートリミング。(アダプターは自動で検出される。low qualityとNのトリミングも実行される)

AdapterRemoval --file1 reads_1.fq --basename output_single --trimns --trimqualities --gzip --threads 8
  • --basename    Default prefix for all output files for which no filename was explicitly set
  • --trimns    If set, trim ambiguous bases (N) at 5'/3' termini [current: off]
  • --gzip     Enable gzip compression [current: off]
  • --threads   Maximum number of threads [current: 1]

3つのファイルが出力される。アダプタートリミングされたリードは~truncated~というファイルに保存される。settings~にはシーケンスデータのstatisticsや実行条件がまとめられている。

 

3、ペアエンドのアダプタートリミング。

AdapterRemoval --file1 reads_1.fq --file2 reads_2.fq --basename output_paired --gzip --trimqualities --collapse --threads 8
  • --trimqualities    If set, trim bases at 5'/3' termini with quality scores <= to --minquality value   [current: off]

 7つのファイルが出力される。ペアエンドの両方のリードのアダプタートリミングがうまくいったリードはtruncated~に保存される。ペアエンドの片方がdiscardされ順番が壊れたリードはcollapse~に保存される。ペアエンドがマージされた場合はoutput_paired.collapsed~に保存され、--trimns や --trimqualitiesオプションの影響で トリムされてマージされたリードはoutput_paired.collapsed.truncatedに保存される。

f:id:kazumaxneo:20180516150936j:plain

 

4、ペアエンドのinterleave fastqのアダプタートリミング。

AdapterRemoval --interleaved --file1 interleaved.fq --basename output_interleaved

 

注:デフォルトではPhred+33 としてトリミングを実行する。昔のPhred+64のfastqなら"--qualitybase 64"をつけて実行する。"--qualitybase-output 33"もつけて実行するとPhred+33に修正して出力してくれる。

 

この他、Githubではv2.1で実装されたDemultiplexing機能を使い、Demultiplexingとアダプタートリミングを同時実行する方法も紹介されています。さらに、double indexのDemultiplexingや、アダプターが複数ある特殊なシーケンスデータのアダプター除去も可能です。詳細はGithubで確認してください。

引用

AdapterRemoval v2: rapid adapter trimming, identification, and read merging

Mikkel Schubert, Stinus Lindgreen, Ludovic Orlando

BMC Res Notes. 2016; 9: 88.

 

AdapterRemoval: easy cleaning of next-generation sequencing reads.

Lindgreen S

BMC Res Notes. 2012 Jul 2;5:337.

 

コード領域のリアライメントによってバリアントコールを改善する ABRA

 

 indel検出を制限するアラインメントエラーおよびリファレンスバイアスを克服するために、多数のリアライメントおよびアセンブリ方法が提案されている。ショートリードのマイクロアライナーは、局所的に組み立てられたバリアントグラフへリードを局所的に再調整する(Homer and Nelson、2010)。 PindelはPattern growth approarchを使用してindelsを検出する(Ye et al、2009)。 Dindelは候補をハプロタイプ候補に再編成し、ベイズ法を使用して最大50 bpの長さのindelsを呼び出す(Albers et al、2011)。GATKのIndelRealignerは、ローカルリアライメントを介して不一致塩基の数を最小限に抑えようとしている(DePristo et al、2011)。全ゲノムde novoアセンブリアプローチには、Fermi(Li、2012)およびCortex Var(Iqbal et al、2012)が含まれる。 SOAPIndelは、ペアリードの片方だけがマッピングされている領域で、ローカルアセンブリとコールを実行する(Li et al、2012)。 Clipping Reveals STructure(CREST)は、ソフトクリップされたリードとローカルアセンブリを実行して、体系的な構造変化を識別する(Wang et al、2010)。 Targeted Iterative Graph Routing Assembler(TIGRA)は、、推定ブレークポイントからターゲットのアセンブリによってコンティグを生成する(Chen et al、2014)。 さらにComplete Genomics(Carnevali et al、2012)およびFoundation Medicine(Framptonら、2013)によって独自のローカルアセンブリ方法も開発されている。

 新たに開発されたABRAというツールは、シーケンスアラインメント/マップ(SAM / BAM)ファイルを入力として受け入れ、リアライメントされたBAMファイルを出力として生成し、バリアントコールアルゴリズムやその他の下流解析の選択に柔軟性をもたらす。グローバルリアラインメントでは、アラインされていないか謝ってマップされたリードを正しいポジションに移動する。 ABRAは、元のリードのアラインメントには存在しない変異を検出し、存在する変異についてアレル頻度推定を改善する。 ABRAは、生殖細胞系および体細胞系の両方の検出を強化するために使用でき、ペアエンドおよびシングルエンドのデータでも機能する。

 

f:id:kazumaxneo:20180516102200j:plain

Githubより転載。

 

インストール

cent osでテストした。

依存

  • Building ABRA requires JDK 7, Maven and g++

本体 Github

https://github.com/mozack/abra

git clone https://github.com/mozack/abra.git 
cd abra/
make

ビルドに失敗したため、リリースにあるビルド済みのプログラムを使い(v. 0.97,linux 64bit)テストした。

https://github.com/mozack/abra/releases

java -jar abra.jar 

$ java -jar -Xmn4G abra.jar 

Java HotSpot(TM) 64-Bit Server VM warning: MaxNewSize (4194304k) is equal to or greater than the entire heap (4194304k).  A new max generation size of 4193792k will be used.

Starting 0.97 ...

Missing required input SAM/BAM file

Missing required input SAM/BAM file

Missing required reference

Missing required target regions

Missing required working directory

Option                                  Description                            

------                                  -----------                            

--adc [Integer]                         Skip regions with average depth        

                                          greater than this value (default:    

                                          100000)                              

--aur                                   Assemble unaligned reads (currently    

                                          disabled).                           

--bwa-ref                               BWA index prefix.  Use this only if    

                                          the bwa index prefix does not match  

                                          the ref option.                      

--ib                                    If specified, write intermediate data  

                                          to BAM file using the intel deflator 

                                          when available.  Use this to speed   

                                          up processing.                       

--in                                    Required list of input sam or bam file 

                                          (s) separated by comma               

--kmer                                  Optional assembly kmer size(delimit    

                                          with commas if multiple sizes        

                                          specified)                           

--lr                                    Search for potential larger local      

                                          repeats and output to specified file 

                                          (only for multiple samples)          

--mad <Integer>                         Regions with average depth exceeding   

                                          this value will be downsampled       

                                          (default: 250)                       

--mapq [Integer]                        Minimum mapping quality for a read to  

                                          be used in assembly and be eligible  

                                          for realignment (default: 20)        

--maxn [Integer]                        Maximum pre-pruned nodes in regional   

                                          assembly (default: 9000)             

--mbq [Integer]                         Minimum base quality for inclusion in  

                                          assembly.  This value is compared    

                                          against the sum of base qualities    

                                          per kmer position (default: 60)      

--mc-mapq [Integer]                     Minimum contig mapping quality         

                                          (default: 25)                        

--mcl [Integer]                         Assembly minimum contig length         

                                          (default: -1)                        

--mer <Double>                          Min edge pruning ratio.  Default value 

                                          is appropriate for relatively        

                                          sensitive somatic cases.  May be     

                                          increased for improved speed in      

                                          germline only cases. (default: 0.02) 

--mnf <Integer>                         Assembly minimum node frequency        

                                          (default: 2)                         

--mpc [Integer]                         Maximum number of potential contigs    

                                          for a region (default: 5000)         

--mur [Integer]                         Maximum number of unaligned reads to   

                                          assemble (default: 50000000)         

--no-debug                              Throttle down debug logging            

--out                                   Required list of output sam or bam file

                                          (s) separated by comma               

--rcf <Double>                          Minimum read candidate fraction for    

                                          triggering assembly (default: 0.01)  

--ref                                   Genome reference location              

--rna                                   Input RNA sam or bam file (currently   

                                          disabled)                            

--rna-out                               Output RNA sam or bam file (required   

                                          if RNA input file specified)         

--single                                Input is single end                    

--sv                                    Enable Structural Variation searching  

                                          (experimental, only supported for    

                                          paired end)                          

--target-kmers                          BED-like file containing target        

                                          regions with per region kmer sizes   

                                          in 4th column                        

--targets                               BED file containing target regions     

--threads <Integer>                     Number of threads (default: 4)         

--umnf [Integer]                        Assembly minimum unaligned node        

                                          frequency (default: 2)               

--working                               Working directory for intermediate     

                                          output.  Must not already exist

 

 

ラン

ランには、bam以外に、target領域のbedとFASTA、bwaのindexを用意する必要がある。target領域は遺伝子のコード領域を想定している。

bwa index ref.fa

bedはBEDtoolsのsortbedを使いポジションソートしておく。

sortBed -i input.bed > sort.bed

準備できたら実行する。

java -Xmn4G -jar abra.jar --in input.bam --out output.bam --ref ref.fa --targets sort.bed --threads 8 --working temp_dir
  • --in One or more input BAMs delimited by comma
  • --out One or more output BAM's corresponding to the set of input BAMs
  • --ref BWA indexed reference genome.
  • --targets BED file describing target assembly regions (Usually corresponds to capture targets) with optional kmer lengths for each region. Targets must be sorted by position in ascending order within each chromosome.

--workingで指定したワーキングディレクトリは上書きできないので、ランに失敗した時はワーキングディレクトリを消してから実行する。

 

 

引用

ABRA: improved coding indel detection via assembly-based realignment.

Mose LE, Wilkerson MD, Hayes DN, Perou CM, Parker JS

Bioinformatics. 2014 Oct;30(19):2813-5.

 

k-merを使いアライメントフリーでバリアントをコールする kestrel

 

 アライメントツールはエラーやバラツキを処理するように設計されているが、リファレンスとは大幅に異なるシーケンスリードを確実に正しい場所に割り当てることはできない。アラインメントの信頼性が低いと、バリアントコールの信頼性が低くなり、真のバリアントと偽のコールを区別することが困難になる。一部の領域ではリードがクリップされたり、まったくマッピングされないことすらある。その結果、アライメントのみに依存するバリアントコールアルゴリズムは、これらのイベントを特徴付けることができない。

 バクテリアゲノムは小さく、比較的単純だが、そのバラツキのために最も難しいターゲットの1つである。水平遺伝子伝達(HGT)は、リファレンスと比較してサンプルを有意に変化させる1つのメカニズムであり、しばしば薬剤耐性または病原性の増加をもたらす(Ochman et al、2000)。バクテリアタイピングおよびサーベイランス法は、既知の対立遺伝子(Li et al、2016)または複数のリファレンス配列(Li et al、2012)のデータベース、全ゲノムアセンブリ(Thomsen et al、2016)を必要とすることが多い。これらのアプローチは、分析を実行したり、対立遺伝子データベースの維持に必要な時間を増加させる。理想的な方法は、単一のリファレンスで細菌の多様性を扱うことで、リファレンスを変更することなく統一的な文脈でゲノムの変更を行うことができるだろう。

 一般的に使用されるいくつかの代替方法があるが、アライメントおよびバリアントコールプロトコルよりも能力が制限されているか、または計算リソースを消費している。Monoploidの生物ではde novoアセンブリの手法でバリアントをコールできるが、コンティグの状態ではカバレッジが"1"に減少するため、ミスコールを簡単に修正することはできなくなる(Olson et al。、2015)。典型的にアセンブリで使用されるde Bruijnグラフを構築するには、数ギガバイトのメモリまたはコンピューティングクラスタが必要になる可能性がある(Li et al、2010)。 Cortex(Iqbal et al、2012)は、de novoアセンブリアプローチの限界を克服するために設計されたが、依然として全ゲノムにわたって組み立てられたde Brujin graphに依存しています。 Platypus(Rimmer et al、2014)はローカルアセンブリを実行するが、バリアントの軌跡が先験的に分からない場合は、リードのマッピングとローカルアセンブリのオーバーヘッドが発生する。

 1つの可能なアプローチは、k-merスペクトルのセット内に含まれる情報をシーケンスデータにわたって使用することである。 K-mersは数値で表すことができ、シーケンスアライメントに依存せず、k-merの計数方法は高速である。今日まで、このようなアプローチ(Gardner and Slezak、2010; Gardner et al。、2013)では、まばらな間隔のSNPが修正されているが、高密度SNPおよびlarge insertionは、リファレンス配列とは非常に異なる多くのk-merを作製するので、そのようなアプローチがどのように効率的に機能するかはすぐには分からない。

著者らは、サンプルシーケンスデータに対してk-mer頻度に厳密に依存するアルゴリズムを構築し、リファレンス配列中のk-merを順序付けしてバリアントをコールする。この方法は、まず、期待されるk-mersの頻度分布の中断を使用して基準を超えるバリエーションの領域を見つけ、その領域に隣接する変更されていないk-mersから開始することによって、k-mer空間による検索を大幅に簡素化する。このアプローチを使用して、Kestrelがどのようにしてゲノムの改変領域を再構築し、高密度変異の領域を解決できるかを示す。

 

インストール

mac os10.12でテストした。

依存

Github

https://github.com/paudano/kestrel

Github リリースにlinux向けバイナリとPDFマニュアルがある。バイナリはosxでも動作する。

https://github.com/paudano/kestrel/releases/tag/1.0.1

./kestrel -h

$ ./kestrel -h

kestrel -f <INPUT_FORMAT> -i <INTERVAL_FILE> -k <KSIZE> -m <OUT_FORMAT>

    -o <OUT_FILE> -p <OUT_FILE> -r <REF_SEQUENCE> -s[<SAMPLE_NAME>]

    -w <WEIGHT_VEC> --alpha=<ALPHA> --ambiregions --ambivar --anchorboth

    --autoflank --byreference --byregion --charset=<CHARSET> --countrev

    --decaymin=<MIN_PROPORTION> --diffq=<QUANTILE> --filespersample=<N_FILES>

    --flank=<LENGTH> --free --hapfmt=<OUT_FORMAT> --lib=<LIB_FILE>

    --liburl=<LIB_URL> --logfile=<LOG_FILE> --loglevel=<LOG_LEVEL> --logstderr

    --logstdout --maxalignstates=<STATES> --maxhapstates=<STATES> --memcount

    --mincount=<MIN_COUNT> --mindiff=<COUNT_DIFF> --minmask=<MIN_MASK>

    --minsize=<MIN_SIZE> --noambigregions --noambivar --noanchorboth --nocountrev

    --nofree --nomemcount --norevregion --normikc --normrefdesc --noseqfilter

    --peakscan=<LENGTH> --quality=<SEQ_FILTER> --revregion --rmikc --rmrefdesc

    --scanlimitfactor=<FACTOR> --seqfilter=<SEQ_FILTER> --stdout --temploc=<TEMP>

    --varfilter=<FILTER_SPEC> SEQ_FILE [SEQ_FILE...]

kestrel -h[<TOPIC>]

 

-f --format <INPUT_FORMAT> (Default = auto)

    Set the input sequence format type. This option determines how the format

    files are read. This option may be set multiple times when reading files

    with different formats. See "count -hreader" for a full list of readers.

 

-h --help [<TOPIC>]

    Get command-line help. If "TOPIC" is specified, then help on a specific

    topic is shown. Try "--help=topics" (or "-htopics") for a list of topics

 

-i --interval <INTERVAL_FILE>

    Reads a file of intervals defining the regions over the reference sequences

    where variants should be detected. If no intervals are specified, variants

    are detected over the full length of each reference sequence. The file type

    is determined by the file name, such as "intervals.bed". BED files are

    supported by Kestrel, and others may be added.

 

-k --ksize <KSIZE> (Default = 31)

    Size of k-mers sequence data is translated to during analysis. If unsure,

    use the default value. If the sequencing error rate is very high, or if the

    reference is very short, a small (e.g. a single short gene), then a smaller

    k-mer size, such as 21, may be useful if the defalt value does not produce

    meaningful results.

 

-m --outfmt <OUT_FORMAT> (Default = vcf)

    Set output format.

 

-o --out <OUT_FILE>

    Set output file name.

 

-p --hapout <OUT_FILE>

    Set haplotype output file name.

 

-r --ref <REF_SEQUENCE>

    Add reference sequences variants will be called against. This can be any

    file that Kestrel can read. The format and character-set options apply to

    reference sequences, but not filters.

 

-s --sample [<SAMPLE_NAME>]

    Set the name of the sample that the next sample files are assigned to. If

    the argument (SAMPLE_NAME) is given, the name of the sample is set to this

    name. If the argument is not given, then the sample name is assigned from

    the name of the first file after this option. Any files on the command-line

    appearing before this option are assigned to a sample and will not be part

    of this sample. If --filespersample was used on the command-line before this

    option, it is reset and files are no longer automatically grouped.

 

-w --weight <WEIGHT_VEC> (Default = (10.0, -10.0, -40.0, -4.0, 0.0))

    Set the alignment weights as a comma-separated list of values. The order of

    weights is match, mismatch, gap-open, gap-extend, and initial score. If

    values are blank or the list has fewer than 5 elements, the missing values

    are assigned their default weight. Each value is a floating-point number,

    and it may be represented in exponential form (e.g. 1.0e2) or as an integer

    in hexadecimal or octal format. Optionally, the list may be surrounded by

    parenthesis or braces (angle, square, or curly).

 

--alpha <ALPHA> (Default = 0.80)

    Set the exponential decay alpha, which controls how quickly the recovery

    threshold declines to its minimum value (see --decaymin) in an active region

    search. Alpha is defined as the rate of decay for every k bases. At k bases

    from the left anchor, the threshold will have declined to alpha * range. At

    every k bases, the threshold will continue to decline at this rate.

 

--ambiregions (Default)

    Allow active regions to cover ambiguous bases, such as N.

 

--ambivar (Default)

    Allow variants over ambiguous bases, such as N.

 

--anchorboth (Default)

    Both ends of an active region (region with variants) must be bordered by

    unaltered k-mers or variants will not be called in it. This option may miss

    variants near the ends of a reference sequence, but it forces stronger

    evidence for the variants that are called.

 

--autoflank (Default)

    When extracting intervals from reference sequences, some bases are extracted

    on both sides of the interval whenever possible. This gives Kestrel more

    bases for active region detection, but it does not otherwise affect variant

    calls. This option tells Kestrel to pick the flank by multiplying 3.50 with

    the k-mer size.

 

--byreference (Default)

    If variant call regions were defined, variant call locations are still

    relative to the reference sequence and not the region.

 

--byregion

    If variant call regions were defined, variant call locations are relative to

    the region and not the reference sequence.

 

--charset <CHARSET> (Default = UTF-8)

    Character set encoding of input files. This option specifies the character

    set of all files following it. The default, "UTF-8", properly handles ASCII

    files, which is a safe assumption for most files. Latin-1 files with values

    greater than 127 will not be properly parsed.

 

--countrev (Default)

    Count reverse complement k-mers in region statistics. This should be set for

    most sequencing protocols.

 

--decaymin <MIN_PROPORTION> (Default = 0.55)

    Set the minimum value (asymptotic lower bound) of the exponential decay

    function used in active region detection as a proportion of the anchor k-mer

    count. If this value is 0.0, k-mer count recovery threshold may decline to

    1. If this value is 1.0, the decay function is not used and the detector

    falls back to finding a k-mer with a count within the difference threshold

    of the anchor k-mer count.

 

--diffq <QUANTILE> (Default = 0.9000)

    If set to a value greater than 0.0, then the k-mer count difference between

    two k-mers that triggers a correction attempt is found dynamically. The

    difference in k-mer counts between each pair of neighboring k-mers over an

    uncorrected reference region is found, and this quantile of is computed over

    those differences. For example, a value of 0.85 means that at most 15% (100%

    - 85%) of the k-mer count differences will be high enough. If this computed

    value is less than the minimum k-mer count difference (--mindiff), then that

    minimum is the difference threshold. This value may not be 1.0 or greater,

    and it may not be negative. If 0.0, the minimum count difference is always

    the minimum threshold (--mindiff).

 

--filespersample <N_FILES> (Default = 0)

    Set the number of input files per sample. For example, reading paired-end

    FASTQ files (2 files per sample) can be simplified by setting this value to

    2. Alternatively, samples can be separated by multiple -s (--sample)

    arguments. The default value, 0, will not automatically group input files.

    If -s (--sample) is read on the command-line, this value is set back to 0.

    Any sequence files found on the command-line before this option are assigned

    to a group.

 

--flank <LENGTH> (Default = k * 3.50)

    When extracting intervals from reference sequences, this many bases are

    extracted on both sides of the interval whenever possible. This gives

    Kestrel more bases for active region detection, but it does not otherwise

    affect variant calls. Set to 0 to disable flanks. If this option is not set,

    Kestrel will determine the appropriate length of flank by multiplying 3.50

    with the k-mer size.

 

--free

    Free resources between processing samples. This may reduce the memory

    footprint of Kestrel, but it may force expensive resources to be recreated

    and impact performance.

 

--hapfmt <OUT_FORMAT> (Default = sam)

    Set haplotype output format. Ignored if a haplotype output file is not set.

 

--lib <LIB_FILE>

    Load a library file. Kestrel can accept external components, and they must

    be packaged on a JAR file.

 

--liburl <LIB_URL>

    Load a library by its URL. Kestrel can accept external components, and they

    must be packaged on a JAR file. This option can access JAR files on the

    local system or stored anywhere the program can access and that can be

    represented as a URL.

 

--logfile <LOG_FILE> (Default = <STDERR>)

    Set log file name.

 

--loglevel <LOG_LEVEL> (Default = WARN)

    Set the log level. Valid levels are ALL, TRACE, DEBUG, INFO, WARN, ERROR,

    and OFF.

 

--logstderr

    Write log messages to standard error instead of a file. Unless redirected,

    this output is written to the the screen.

 

--logstdout

    Write log messages to standard output instead of a file. Unless redirected,

    this output is written to the the screen.

 

--maxalignstates <STATES> (Default = 15)

    Set the maximum number of alignment states. When haplotype assembly branches

    into more than one possible sequence, the state of one is saved while

    another is built. When the maximum number of saved states reaches this

    value, the least likely one is discarded.

 

--maxhapstates <STATES> (Default = 15)

    Set the maximum number of haplotypes for an active region. Alignments can

    generate more than one haplotype, and with noisy sequence data or

    paralogues, many haplotypes may be found. This options limits the amount of

    memory that can be consumed in these cases.

 

--memcount

    K-mer counts from each sample will be stored in memory. This option assumes

    that samples are relatively small or the machine has enough memory to handle

    the counts. Note that the JVM might need to be run with additional memory

    (-Xmx option) to support this option.

 

--mincount <MIN_COUNT> (Default = 5)

    Set the minimum k-mer count for processing samples. K-mers with a count less

    than this value will be discarded. Sequence read errors produce many

    erroneous k-mers, and this slows the process of variant calling

    significantly.

 

--mindiff <COUNT_DIFF> (Default = 5)

    Set the minimum k-mer count difference for identifying active regions. When

    the count between neighboring k-mer counts is this or greater, Kestrel will

    treat it as a region where a variant may occur.

 

--minmask <MIN_MASK> (Default = 0)

    Size of k-mer minimizers or 0 to disable processing by minimizers. The

    minimizer of a k-mer is determined by taking all sub-k-mers of a given size

    (set by this option) from a k-mer and its reverse complement and choosing

    the lesser of the sub-k-mers. Sub-k-mers are XORed with this mask while

    comparing them, but the minimizer is not XORed (it is still a sub-k-mer of

    the original k-mer). This option can be used to break up large minimizer

    groups due to low-complexity k-mers when minimizers are used.

 

--minsize <MIN_SIZE> (Default = 15)

    Minimizers group k-mers in the indexed k-mer count (IKC) file generated by

    Kestrel when reading sequences, and this parameter controls the size of the

    minimizer.

 

--noambigregions

    An active region may not span any base that is not A, C, G, T, or U.

 

--noambivar

    A variant may not span any base that is not A, C, G, T, or U.

 

--noanchorboth

    An active region (region with variants) must be bordered on at least one

    side by an unaltered k-mers, but it may extend to the end of the sequence.

    This will allow Kestrel to find variants less than a k-mer from the ends,

    but the evidence supporting these variants is weaker.

 

--nocountrev

    Do not include the reverse complement of k-mers in read depth estimates. If

    all sequence reads are in the same orientation as the reference, then this

    option should be used.

 

--nofree (Default)

    Retain resources between samples. This may use more memory, but it will

    avoid re-creating expensive resources between samples.

 

--nomemcount (Default)

    K-mer counts for each sample are offloaded to an indexed k-mer count file.

    This option reduces the memory demand of Kestrel.

 

--norevregion

    When set, regions variants are called on are always in the same orientation

    as the reference sequence. The stranded-ness of defined intervals is

    ignored.

 

--normikc

    Do not remove the indexed k-mer count (IKC) file for each sample.

 

--normrefdesc

    Use the full sequence name as it appears in the reference sequence file.

    FASTA files often include a description after the sequence name, and with

    this option, it becomes part of the full sequence name. If using an interval

    file, the full sequence name and description must match the sequence file.

 

--noseqfilter (Default)

    Turn off sequence filtering for all files following this option. If

    --seqfilter or --quality was specified, this option disables sequence

    filtering. These options together make it possible to specify filtering for

    some files and disable filtering for others.

 

--peakscan <LENGTH> (Default = 7)

    Reference regions with sequence homology in other regions of the genome may

    contain k-mers with artificially high frequencies from adding counts for

    k-mers that appear in both regions. This causes a peak in the k-mer

    frequencies over the reference, and it can trigger an erroneous

    active-region scan for variants. When encountering a difference, Kestrel

    will scan forward this number of k-mers looking for a peak in the k-mer

    frequencies. If the frequencies drop back down to the original range, the

    active-region scan is not performed. This keeps Kestrel from erroneously

    searching large regions of the reference. Setting this value to 0 disables

    peak detection.

 

--quality <SEQ_FILTER>

    This option is an alias for "seqfilter"

 

--revregion

    When set, reverse complement reference regions that occur on the negative

    strand. Only itervals defined with on the negative strand are altered.

 

--rmikc (Default)

    Remove the indexed k-mer count (IKC) for each sample after kestrel runs.

 

--rmrefdesc (Default)

    When set, remove the description from reference sequence names. The

    descirption is everything that occurs after the first whitespace character.

    FASTA files often have a sequence name and a long description separated by

    whitespace. This option ensures that the sequence name matches in the FASTA

    and an interval file, if used.

 

--scanlimitfactor <FACTOR> (Default = 5.0)

    Set a limit on how long an active region may be. This is computed by

    multiplying the k-mer size by this factor and adding the maximum length of a

    gap. The computed limit will be adjusted so that active regions are at least

    large enough to capture a SNP in cases where the maximum gap length is 0.

    Setting this to a low value or "min" will set the limit so that it is just

    large enough to catch SNPs and deletions, but it will miss large deletions

    if another variant is within the k-mer size window. Setting this to a high

    value or "max" lifts the restrictions on active region lengths, and this may

    cause the program to take an excessive amount of time and memory trying to

    solve arbitrarily long active regions.

 

--seqfilter <SEQ_FILTER>

    Filter sequences as they are read and before k-mers are extracted from them.

    Some sequence readers can filter or alter reads at runtime. The most common

    filter is a quality filter where low-quality bases are removed. The filter

    specification is a filter name followed by a colon (:) and arguments to the

    filter. If a filter name is not specified, then the "sanger" quality filter

    is assumed. For example, "sanger:10" and "10" will filter k-mers with any

    base quality score less than 10. The sequence filter specification is set

    for all files appearing on the command-line after this option. To turn off

    filtering once it has been set, files following --noseqfilter will have no

    filter specification.

 

--stdout (Default)

    Write output to standard output instead of a file. Unless redirected, this

    output is written to the the screen.

 

--temploc <TEMP> (Default = Output location)

    The location where segments are offloaded. This argument must be a directory

    or the location for a new directory. Parent directories will be created as

    needed.

 

--varfilter <FILTER_SPEC>

    Add a variant filter specification. The argument should be the name of the

    filter, a colon, and the filter arguments. The correct filter is loaded by

    name and the filter arguments are passed to it.

 

SEQ_FILE

    Input sequence file.

 

kestrel version: 1.0.1

 

ラン

リファレンスとfastqを指定して実行する。

kestrel -r ref.fa -k 31 -o variant.vcf input.fq
  • -k    Size of k-mers sequence data is translated to during analysis  (Default = 31)
  • -r     Reference sequences variants will be called against.

ペアエンド情報は使わないので、ペアエンドシーケンスデータを使うならマージしてから使う。

 

引用

Mapping-free variant calling using haplotype reconstruction from k-mer frequencies

Audano PA, Ravishankar S, Vannberg FO

Bioinformatics. 2018 May 15;34(10):1659-1665.

 

ONTのロングリードのアセンブリとポリッシュ

 

 2018年2月のNature CommunicationsにシロイヌナズナのゲノムをONTのロングリードを使ってアセンブリした論文( PCR-free paired-end readsでpolish)が出ている(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5803254/)。ONTの MinION R9.4 flowcell (FLO-MIN106)で48時間シーケンスし、平均サイズ11.4kbのリードが合計3.4Gb得られたらしい(およそx30)。それからminiasmを使い一般的な4 Cores, 16 Gb RAMハードウエア(15インチのmacbook proのようなマシン)でアセンブルして、62 contigs(N50 12.3 Mb)得た、となっている。良好なパフォーマンスである。ただしエラーを許容するminiasmでアセンブリすると大量に残ってしまうため、著者らは、それからRaconを3回かけてコンセンサス配列を出力し、さらに最後にMiseqのPCR free ショートリード250x2を使ってPilonでポリッシュして、最終出力としている。十分なcontiguityのcontigが得られるため、オーサーはIntroductionの最後に、ロングリードでドラフトを作り構造変化(SV)を検出するのが一番シンプルでstraightforwardな手法であろうと主張してる。

 

Oxford nanoporeの公式レポジトリでも、このフローの精度が調べられている。

nanoporetech/ont-assembly-polish

https://github.com/nanoporetech/ont-assembly-polish

 

こちらも参考になります。 

Assembly using miniasm+racon

http://inf-biox121.readthedocs.io/en/2016/Assembly/practicals/07_Assembly_using_minasm+racon.html

 

ONTのリードだけでアセンブルする手法は上記の論文以外に様々報告されており、精度を上げるアイデアも提案されている(リンク)。ここでは、上記の論文のようなminiasm(紹介)とracon(紹介)を使ったアセンブリとポリッシュのワークフローをまとめておく。Raconを複数回ランして精度が上がるように工夫している。

 

1、SRA_toolkitを使い上記のシーケンスデータ(ONT & illumina)をダウンロードする。

#ONTリードのダウンロード
prefetch ERR2173373 #このデータ
fastq-dump ~/ncbi/public/sra/ERR2173373.sra -O ONT #fastqに解凍

#illumina PCR free short read(250x20
prefetch ERR2173372 #このデータ
fastq-dump --split-files ~/ncbi/public/sra/ERR2173372.sra -O illumina 
#fastqに解凍

上記の論文では、この時点でナズナのリファレンスゲノムにリードをアライメントし、アライメントされなかったリードはlow qualityで、GC含量が大きく異なるようなリードであると見出している。

 

2、ダウンロードされたリードの分析

ショートリードはseqkit(紹介)で簡易チェック。

seqkit stats simulated.fa

$ seqkit stats ONT/ERR2173373.fastq 

file                  format  type  num_seqs        sum_len  min_len   avg_len  max_len

ONT/ERR2173373.fastq  FASTQ   DNA    300,071  3,421,779,258        5  11,403.2  269,087

 

ロングリードはNanofiltを使いサイズとクオリティを分析(紹介)。

NanoPlot --fastq ONT/ERR2173373.fastq --loglength -t 12

f:id:kazumaxneo:20180324212546j:plain

f:id:kazumaxneo:20180324212549j:plain

数千bp付近に妙なピークがある。 コピー数が多いクロロプラストの配列かもしれない(該当するサイズを抽出し、データベース検索して集計すれば調べられる)。

 

ここではアダプター除去や、クオリティトリミングは考えない。de novoアセンブリを実行する。

アセンブリ

3、minimap(紹介)でアセンブリ

minimap -Sw5 -L100 -m0 -t 20 ONT/ERR2173373.fastq ONT/ERR2173373.fastq |gzip -1 > minimap_out.paf.gz #論文と同じパラメータ

miniasm -f ONT/ERR2173373.fastq minimap_out.paf.gz > miniasm_out.gfa

awk '/^S/{print ">"$2"\n"$3}' miniasm_out.gfa | fold > raw_assembly.fasta

 トリミングなしでアセンブルすると、11,425,724bpのcontigが1つだけできた(染色体は5本ある、半数体ゲノムサイズはおよそ1.2億bp)。

 

ポリッシュ 

 4、Racon 1回目(Racon紹介)。Raconのランにはall reads vs all reads のマッピング部位を示したPAFファイル(Pairwise mApping Format) が必要なので、3でアセンブリしたcontigを使いminimapを動かす

minimap raw_assembly.fasta ONT/ERR2173373.fastq > racon1.paf

#Raconでコンセンサス配列を出力。
racon -t 20 ONT/ERR2173373.fastq racon1.paf raw_assembly.fasta racon_polished1.fasta

  

5、Racon 2回目。

minimap racon_polished1.fasta ONT/ERR2173373.fastq > racon2.paf
racon -t 20 ONT/ERR2173373.fastq racon2.paf racon_polished1.fasta racon_polished2.fasta

 

6 --Racon 3回目。

minimap racon_polished2.fasta ONT/ERR2173373.fastq > racon3.paf
racon -t 20 ONT/ERR2173373.fastq racon3.paf racon_polished2.fasta racon_polished3.fasta

 

7、Pilon(紹介)を使いilluminaのハイクオリティリードでポリッシュする。まずbamを作る。

#マッピング
bwa mem -t 20 racon_polished3.fasta illumina/ERR2173372_1.fastq illumina/ERR2173372_2.fastq |samtools view -@ 20 -bS - > aln.bam

#ソート
samtools sort -@ 20 aln.bam > aln_sorted.bam

#index
samtools index aln_sorted.bam

 

8、ポリッシュを実行。

#pilonの実行ファイルをダウンロード
wget https://github.com/broadinstitute/pilon/releases/download/v1.22/pilon-1.22.jar

java -Xmx16G -jar pilon-1.22.jar --genome racon_polished3.fasta --frags aln_sorted.bam --vcf --tracks --changes --verbose --threads 20 --outdir output_dir

 

引用

High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell

Todd P. Michael,Florian Jupe, Felix Bemm, Timothy Motley, Justin P. Sandoval, Christa Lanz, Olivier Loudet, Detlef Weigel, and Joseph R. Ecker

Nat Commun. 2018; 9: 541.

 

 

バリアントのコールと可視化のパイプライン MutScan

 

 次世代シーケンシング(NGS)は何千もの突然変異を検出することができる。しかし、一部のアプリケーションでは、これらのうちのほんのわずかなものが対象のターゲットである。 NGS技術によるがんの個人化された医療検査のようなアプリケーションでは、臨床医と遺伝カウンセラーは、通常、薬物治療可能な突然変異の検出に焦点を当てている[論文より ref.1](一部略)。これらの突然変異は、患者の無細胞腫瘍DNA(ctDNA)のディープシーケンスによって検出することができる[ref.3]。しかし、ctDNAシーケンスでコールされるバリアントの突然変異対立遺伝子頻度(MAF)は非常に低い。典型的には、MAFは通常5%以下であり、0.1%という低い値であることすらあり得る[ref.4 pubmed]。そのような低いMAFを用いた突然変異の検出の必要性は、ctDNAシーケンスデータを分析する高感度な方法の開発を推進する[ref.4]。

 NGSデータの通常の突然変異検出パイプラインは、通常、各ステップごとに異なるツールを含む。例えば、著者たちの通常の腫瘍変異コールパイプラインでは、データ前処理のためのAfter [ref.5 紹介]、マッピングのためのBWA [ref.6]、パイプアップ生成のためのSamtools [ref.7]、バリアントコールのためのVarScan2 [ref.8]など多くの補助ツールが必要である。これらのステップで使用されるさまざまなツールは、適用されるフィルタが異なるために情報が失われる可能性があり、最終的に特にMAFが低いものが誤検出される可能性がある。このタイプのデータ分析によるfalse negativesは、患者のよりよい治療の機会を逃す可能性があるため、臨床応用では許容できない。

 対照的に、高価だが効果のない治療法を導入する可能性があるため、鍵突然変異のfalse positives検出も避けるべきである。false positivesによる間違った治療法は重大な副作用を引き起こすことさえある[ref.10 pubmed]。従来のNGSパイプラインは多くの置換とINDELを検出することができるが、必然的に誤検出を引き起こす。特に、アライナーのリファレンスゲノムへの不正確なマッピングのために、ゲノムの高度に反復する領域において偽陽性突然変異が検出され得る。この誤ったコール頻度を減らすには、すべての重要な突然変異を検証する必要がある[ref.11]。バリアントの視覚化は、突然変異の信頼性を手動で確認する重要な方法である。 IGV [ref.12]やGenomeBrowseなどのツールを使用してバリアントのビジュアライゼーションを行うことができるが、これらのツールには低速で非効率なBAMファイルの操作が必要である。特に、非常にディープなシーケンスデータにおいて低いMAF変異を視覚化すると、IGVまたはGenomeBrowseは、突然変異したリードを何千ものリード中に配置することが困難になるため不便である。したがって、高速で軽量でクラウドに優しいバリアントの視覚化ツールが必要である。

 ここで紹介されているツールMutScanは、これらの問題に対処するために特別に設計されている。エラー耐性のある文字列検索アルゴリズムを基に構築されており、ローリングハッシュ[ref.13]とブルームフィルタ[14]を使用して速度を最適化している。 MutScanは、CSVファイルまたはプログラムであらかじめ定義されているターゲット変異を検出するためリファレンスフリーモードでも実行できる。 VCFファイルとそれに対応するリファレンスゲノムのFastAファイルを提供することで、MutScanはこのVCF内のすべてのバリアントをスキャンし、各バリアント用のHTMLページをレンダリングすることによって視覚化することができる。

 

 

f:id:kazumaxneo:20180513183822j:plain

MutScanのワークフロー。論文より転載。

 

特徴

  • Ultra sensitive, guarantee that all reads supporting the mutations will be detected
  • Can be 50X+ faster than normal pipeline (i.e. BWA + Samtools + GATK/VarScan/Mutect).
  • Very easy to use and need nothing else. No alignment, no reference genome, no variant call, no...
  • Contains built-in most actionable mutation points for cancer-related mutations, like EGFR p.L858R, BRAF p.V600E...
  • Beautiful and informative HTML report with informative pileup visualization.
  • Multi-threading support.
  • Supports both single-end and pair-end data.
  • For pair-end data, MutScan will try to merge each pair, and do quality adjustment and error correction.
  • Able to scan the mutations in a VCF file, which can be used to visualize called variants.
  • Can be used to filter false-positive mutations. i.e. MutScan can handle highly repetive sequence to avoid false INDEL calling.

 

想定されるシナリオ

  • you are interested in some certain mutations (like cancer drugable mutations), and want to check whether the given FastQ files contain them.
  • you have no enough confidence with the mutations called by your pipeline, so you want to visualize and validate them to avoid false positive calling.
  • you worry that your pipeline uses too strict filtering and may cause some false negative, so you want to check that in a fast way.
  • you want to visualize the called mutation and take a screenshot with its clear pipeup information.
  • you called a lot of INDEL mutations, and you worry that mainly they are false positives (especially in highly repetive region)
  • you want to validate and visualize every record in the VCF called by your pipeline.

ガン原遺伝子の探索ではbuilt-inのデータベースを使えるので、raw シーケンスデータから新規に変異解析を実行できます。それ以外のケースでは、あらかじめ別のツールで解析してvcfを得る必要があります。本ツールにvcfを入力することで、簡単なフィルタリングと結果の可視化を行えます。

 

GithubにSample reportが用意されている。

http://opengene.org/MutScan/report.html

 

インストール

mac os 10.12でテストした。

Github

GitHub - OpenGene/MutScan: Detect and visualize target mutations by scanning FastQ files directly

git clone https://github.com/OpenGene/mutscan.git
cd mutscan
make
sudo make install

mutscan

$ mutscan 

usage: mutscan --read1=string [options] ... 

options:

  -1, --read1         read1 file name (string)

  -2, --read2         read2 file name (string [=])

  -m, --mutation      mutation file name, can be a CSV format or a VCF format (string [=])

  -r, --ref           reference fasta file name (only needed when mutation file is a VCF) (string [=])

  -h, --html          filename of html report, default is mutscan.html in work directory (string [=mutscan.html])

  -t, --thread        worker thread number, default is 4 (int [=4])

  -S, --support       min read support for reporting a mutation, default is 2 (int [=2])

  -k, --mark          when mutation file is a vcf file, --mark means only process the records with FILTER column is M

  -l, --legacy        use legacy mode, usually much slower but may be able to find a little more reads in certain case

  -s, --standalone    output standalone HTML report with single file. Don't use this option when scanning too many target mutations (i.e. >1000 mutations)

      --simplified    simplified mode uses less RAM but reports less information. This option can be auto/on/off, by default it's auto, which means automatically enabled when processing large FASTQ with large VCF. (string [=auto])

  -v, --verbose       enable verbose mode, more information will be output in STDERR

  -?, --help          print this message

——

 

ラン

テストデータのダウンロード。

http://opengene.org/dataset.html

がん原遺伝子を探索するなら、built-inのデータベースと称号するため、 ペアエンドデータを指定するだけでランできる(MutScan contains a built-in list with most actionable gene mutations for cancer diagnosis [18]. 論文より )。

mutscan -1 R1.fq.gz -2 R2.fq.gz -t 8
  • -1     read1 file name (string)  
  • -2     read2 file name (string [=])
  • -t      worker thread number, default is 4 (int [=4])

シングルエンドのシーケンスデータは"-1"で指定する。出力ファイル名を指定するなら-hフラグを立てて、"-h output.html"などと書く。

ラン結果はhtmlで出力される。

f:id:kazumaxneo:20180513175338j:plain

htmlを開く。コール部位がまとめられている。

f:id:kazumaxneo:20180513190302j:plain

1つ開いてみる。真ん中のCが変異部位。6リード変異をサポートしている。

f:id:kazumaxneo:20180513190630j:plain

書いてある通り文字の色でベースクオリティを表している。赤がlow qualityなベース。

リードをクリックするとraw sequence readが表示される。

f:id:kazumaxneo:20180513190847j:plain

 デフォルトでは最低2以上バリアントをサポートするリードがないと出力しない。フィルタリング感度を変えるには-Sフラグをつける(e.g., "-S 1")。

 

リダイレクトするとプレーンテキストで出力される。

mutscan -1 R1.fq.gz -2 R2.fq.gz -t 8 > result.txt

> head result.txt 

 $ head result.txt 

 

---------------

NRAS-neg-1-115258748-2-c.34G>T-p.G12C-COSM562 GGATTGTCAGTGCGCTTTTCCCAACACCAC A TGCTCCAACCACCACCAGTTTGTACTCAGT chr1

1, pos: 86, distance: 0, reverse

@NB551106:59:HTFV3BGX2:1:11102:9413:4074 1:N:0:AGTCAA

GGGCCTCACCTCTATGGTGGGATCATATTCATCTACAAAGTGGTTCTGGATTAGCT GGATTGTCAGTGCGCTTTTCCCAACACCAC A TGCTCCAACCACCACCAGTTTGTACTCAGT CATTTCACACCAGCAAGAACCTGTTGGAAACC

-

AAEEEAEEEEEEEEEAAEAAAEEEEEE/EEEAEEEEAEEEEEEEAEEEAEEEEEEE AEEAEEEEEEEEEEEEEEEAEEEEEEEEEE E AEEEEEEEEEEEEEEEEE/EEEEEEEEEEE EEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA

2, pos: 86, distance: 0, reverse

@NB551106:59:HTFV3BGX2:1:11102:9413:4074 1:N:0:AGTCAA

 

がん原遺伝子のデータベース以外を可視化するには、vcfを指定してランする。vcf指定時はリファレンスファイルも指定する必要がある(.csvならリファレンスは必要ない)。

mutscan -1 R1.fq.gz -2 R2.fq.gz -t 8 -r ref.fa -m input.vcf
  • -m      mutation file name, can be a CSV format or a VCF format
  • -r       reference fasta file name (only needed when mutation file is a VCF) (string [=])

vcfを与える場合、built-inのデータベースに依存しないため、どのような変異でも分析して可視化できる(記載はないがヒト以外でも動作する)。SVも可視化できそうだが、50~100-bp以上のSVでは横長になって視認性が悪いので、SV専用の可視化ツールの方が適していると思われる(samplotSVPV)。

 

引用

MutScan: fast detection and visualization of target mutations by scanning FASTQ data

Shifu Chen,Tanxiao Huang, Tiexiang Wen, Hong Li, Mingyan Xu, and Jia Gu

BMC Bioinformatics. 2018; 19: 16.