macでインフォマティクス

macでインフォマティクス

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

ターゲットアンプリコンシーケンシングのプライマーを除く cutPrimers

 

 リードからのプライマーの切断は、ターゲットアンプリコンのNGSデータを処理する上で重要なステップである。既存のツールは、リードから1つまたはいくつかのプライマー/アダプター配列を切断し、そして出現するそれらすべてを除去するように適合されている。また、既存のツールのほとんどは、kmerを使用しており、プライマーの一部または研究された遺伝子配列を有するプライマーのみを切断することができる。このため、そのようなプログラムを使用すると、誤ったトリミング、カバレッジの縮小、および誤検出増加につながる。本著者らはリードから任意の数のプライマーを正確に切断するcutPrimersという新しいツールを開発した。 BRCA1 / 2遺伝子の研究中に得られたシークエンスリードを用いて、cutPrimersをcutadapt、AlienTrimmer、およびBBDukと比較した。それらすべてが、cutPrimerでトリミングされたリードと比較して、少なくとも2つのアンプリコンのカバレッジが許容できないレベル(<30リード)に減少するようにリードをトリミングした。同時に、TrimmomaticとAlienTrimmerはすべてのプライマー配列をカットしたため、残りのリードの長さは予想以上に短くなった。

 

インストール

ubuntu16.04でテストした。

依存

Linux

cutPrimers works on the Python3+ and requires the following packages:

本体 Github

git clone https://github.com/aakechin/cutPrimers.git
cd cutPrimers/

> python3 cutPrimers.py -h

# python3 cutPrimers.py -h

usage: cutPrimers.py [-h] [--readsFile_r1 READSFILE1]

                     [--readsFile_r2 READSFILE2] [--bam-file BAMFILE]

                     [--coordinates-file COORDSFILE]

                     [--out-bam-file OUTBAMFILE]

                     [--out-untrimmed-bam-file OUTUNTRIMMEDBAMFILE]

                     [--minimal-read-length MINREADLEN] [--hard-clipping]

                     --primersFileR1_5 PRIMERSFILER1_5

                     [--primersFileR2_5 PRIMERSFILER2_5]

                     [--primersFileR1_3 PRIMERSFILER1_3]

                     [--primersFileR2_3 PRIMERSFILER2_3]

                     [--trimmedReadsR1 TRIMMEDREADSR1]

                     [--trimmedReadsR2 TRIMMEDREADSR2]

                     [--untrimmedReadsR1 UNTRIMMEDREADSR1]

                     [--untrimmedReadsR2 UNTRIMMEDREADSR2]

                     [--primersStatistics PRIMERSSTATISTICS]

                     [--error-number ERRNUMBER]

                     [--primer-location-buffer PRIMERLOCBUF]

                     [--min-primer3-length MINPRIMER3LEN] [--primer3-absent]

                     [--identify-dimers IDIMER] [--threads THREADS]

 

This script cuts primers from reads sequences

 

optional arguments:

  -h, --help            show this help message and exit

  --readsFile_r1 READSFILE1, -r1 READSFILE1

                        file with R1 reads of one sample

  --readsFile_r2 READSFILE2, -r2 READSFILE2

                        file with R2 reads of one sample

  --bam-file BAMFILE, -bam BAMFILE

                        BAM-file in which you want to cut primers from reads

  --coordinates-file COORDSFILE, -coord COORDSFILE

                        file with coordinates of amplicons in the BED-format

                        (without column names and locations of primers):

                        chromosome | start | end. It is necessary for cutting

                        primer sequences from BAM-file. Its order should be

                        the same as for files with primer sequences

  --out-bam-file OUTBAMFILE, -outbam OUTBAMFILE

                        name of file for output BAM-file with reads

  --out-untrimmed-bam-file OUTUNTRIMMEDBAMFILE, -outbam2 OUTUNTRIMMEDBAMFILE

                        name of file for output BAM-file with untrimmed reads.

                        It is not required. If you do not use this parameter,

                        all untrimmed reads will be lost

  --minimal-read-length MINREADLEN, -minlen MINREADLEN

                        minimal length of read after trimming. Default: 10

  --hard-clipping, -hard

                        use this parameter, if you want to trim reads wuith

                        hard clipping. By default, primer sequences are

                        trimmed with soft-clipping

  --primersFileR1_5 PRIMERSFILER1_5, -pr15 PRIMERSFILER1_5

                        fasta-file with sequences of primers on the 5'-end of

                        R1 reads

  --primersFileR2_5 PRIMERSFILER2_5, -pr25 PRIMERSFILER2_5

                        fasta-file with sequences of primers on the 5'-end of

                        R2 reads. Do not use this parameter if you have

                        single-end reads

  --primersFileR1_3 PRIMERSFILER1_3, -pr13 PRIMERSFILER1_3

                        fasta-file with sequences of primers on the 3'-end of

                        R1 reads. It is not required. But if it is determined,

                        -pr23 is necessary

  --primersFileR2_3 PRIMERSFILER2_3, -pr23 PRIMERSFILER2_3

                        fasta-file with sequences of primers on the 3'-end of

                        R2 reads

  --trimmedReadsR1 TRIMMEDREADSR1, -tr1 TRIMMEDREADSR1

                        name of file for trimmed R1 reads

  --trimmedReadsR2 TRIMMEDREADSR2, -tr2 TRIMMEDREADSR2

                        name of file for trimmed R2 reads

  --untrimmedReadsR1 UNTRIMMEDREADSR1, -utr1 UNTRIMMEDREADSR1

                        name of file for untrimmed R1 reads. If you want to

                        write reads that has not been trimmed to the same file

                        as trimmed reads, type the same name

  --untrimmedReadsR2 UNTRIMMEDREADSR2, -utr2 UNTRIMMEDREADSR2

                        name of file for untrimmed R2 reads. If you want to

                        write reads that has not been trimmed to the same file

                        as trimmed reads, type the same name

  --primersStatistics PRIMERSSTATISTICS, -stat PRIMERSSTATISTICS

                        name of file for statistics of errors in primers. This

                        works only for paired-end reads with primers at 3'-

                        and 5'-ends

  --error-number ERRNUMBER, -err ERRNUMBER

                        number of errors (substitutions, insertions,

                        deletions) that allowed during searching primer

                        sequence in a read sequence. Default: 5

  --primer-location-buffer PRIMERLOCBUF, -plb PRIMERLOCBUF

                        Buffer of primer location in the read from the start

                        or end of read. If this value is zero, than cutPrimers

                        will search for primer sequence in the region of the

                        longest primer length. Default: 10

  --min-primer3-length MINPRIMER3LEN, -primer3len MINPRIMER3LEN

                        Minimal length of primer on the 3'-end to trim. Use

                        this parameter, if you are ready to trim only part of

                        primer sequence of the 3'-end of read

  --primer3-absent, -primer3

                        if primer at the 3'-end may be absent, use this

                        parameter

  --identify-dimers IDIMER, -idimer IDIMER

                        use this parameter if you want to get statistics of

                        homo- and heterodimer formation. Choose file to which

                        statistics of primer-dimers will be written. This

                        parameter may slightly decrease the speed of analysis

  --threads THREADS, -t THREADS

                        number of threads

 

 

テストラン

入力のfastqを指定する。

cd cutPrimers/
python3 cutPrimers.py \
-r1 example/1_S1_L001_R1_001.fastq.gz \
-r2 example/1_S1_L001_R2_001.fastq.gz \
-pr15 example/primers_R1_5.fa \
-pr25 example/primers_R2_5.fa \
-pr13 example/primers_R1_3.fa \
-pr23 example/primers_R2_3.fa \
-tr1 example/1_r1_trimmed.fastq.gz \
-tr2 example/1_r2_trimmed.fastq.gz \
-utr1 example/1_r1_untrimmed.fastq.gz \
-utr2 example/1_r2_untrimmed.fastq.gz \
-t 2
  • -r1     file with R1 reads of one sample
  • -r2     file with R2 reads of one sample
  • -pr15    fasta-file with sequences of primers on the 5'-end of R1 reads
  • -pr25    fasta-file with sequences of primers on the 5'-end of R2 reads. Do not use this parameter if you have single-end reads
  • -pr13     fasta-file with sequences of primers on the 3'-end of R1 reads. It is not required. But if it is determined, -pr23 is necessary
  • -pr23    fasta-file with sequences of primers on the 3'-end of R2 reads
  • -tr1      name of file for trimmed R1 reads
  • -tr2      name of file for trimmed R2 reads
  • -utr1    name of file for untrimmed R1 reads. If you want to write reads that has not been trimmed to the same file as trimmed reads, type the same name
  • -utr2    name of file for untrimmed R2 reads. If you want to write reads that has not been trimmed to the same file as trimmed reads, type the same name
  • -t     number of threads

出力のtrimmed/untrimmed fastqはgzip圧縮される。

 

bamを指定してprimerを除くこともできます。helpを確認して下さい。

 

引用

cutPrimers: A New Tool for Accurate Cutting of Primers from Reads of Targeted Next Generation Sequencing
Kechin A, Boyarskikh U, Kel A, Filipenko M

J Comput Biol. 2017 Nov;24(11):1138-1143

 

メタゲノムのraw fastqからantibiotic resistance genesを再構成する fARGene

2019 5/20 関連ツール追記

2022/06/02 インストール追記

 

 抗生物質耐性菌による感染は世界的に増加しており、公衆衛生に大きな脅威をもたらしている[ref.1]。抗生物質耐性は細菌種の固有の特性である場合があるが、その臨床的意味において、それは既存の染色体DNAの突然変異によって、またはより一般的には遺伝子の水平移動によって獲得された形質である[ref.2]。環境コミュニティは非常に多様な抗生物質耐性遺伝子(ARG)を抱えていて[ref.3、4、5]、ヒトまたは動物において直接または共生細菌を介して病原体に広がる可能性がある。実際、臨床的に関連性のあるARGの多くは環境細菌に由来すると考えられており、一緒になって、大きく、そしてほとんど探求されていない抵抗性貯留層を構成している[ref.6]。さらに、高濃度の抗生物質にさらされる環境で特に顕著な現象[ref.8、9、10]で、強い選択圧がARGの豊富さと多様性を豊かにし得る[ref.7]。それゆえ、多くの特徴付けされていないARGを含む、環境および共生細菌群集におけるレジストソームを調査することが不可欠である。これにより、ARGのヒト病原体への進化および動員の背後にあるプロセスの理解が深まり、臨床現場に到達する前の早期監視および封じ込め行動が促進される。

 ショットガンメタゲノミクスは、収集されたゲノムのランダムフラグメントのシーケンシングを通して細菌群集の分析を可能にする[ref.11, 12]。環境細菌種の1%未満が標準的な方法で培養可能と見なされる。メタゲノムはそれゆえ、培養ベースのアプローチと比較して、そのARGを含むコミュニティの補完的な見方を提供する[ref.13]。大規模並列シーケンシングのキャパシティ増加に伴い、ショットガンメタゲノミクスの実行が急速に安価になり、それに応じてデータのアクセスが容易になった[ref.14]。しかし、メタゲノムデータは非常に細分化されており、幅広いアーティファクトとノイズを特徴としている[ref.15]。これは、特により豊富さが少ない細菌種および株にとって、より大きなゲノム領域のアセンブリを困難にする。これはARGに特に当てはまり、ARGは、例えばインテグロン、トランスポゾン、およびプラスミド中のリペティティブエレメントおよび複数の文脈を有する領域にしばしば存在する。したがって、ARGはメタゲノムデータから再構築するのが難しいことでよく知られる[ref.8、16、17]。したがって、メタゲノムデータセットにおけるこれまでに特徴付けられていなかったARGの多様性は、既存の分析パイプラインでは見過ごされがちである。

 細菌DNA配列中の抗生物質耐性遺伝子の同定は、ARGリファレンスデータベースに対する相同性検索を通してしばしば行われる[ref.18-20]。 ARDB [21]、SARG [ref.22]、CARD [ref.23]、ResFinder [ref.20]など、いくつかのリファレンスデータベースがこの目的のために設計されている。これらおよび他のリファレンスデータベースが何千もの特徴付けられたARGを含んでいても、それらは全抵抗性のごく一部を反映しているだけである[ref.6]。さらに、アノテーションパイプラインの大多数は既知のARGを見つけるために開発されたものであり、既知のARGとの類似性が低い機能耐性遺伝子を正確に同定するために最適化されていないか、あるいは不可能であることが多い。これは、フラグメント長が短いためにアノテーションを付けることが特に困難なメタゲノムデータに特に当てはまる。[ref.24-27]。例えばARGs-OAP [ref.22]、GROOT [ref.28]、AmrPlusPlus [ref.29]、MEGAN [ref.30]、ARIBA [ref.31]など、いくつかの方法がこの目的のために開発されている。ただし、これらの方法では、既知のARGを使用してデータベースに厳密に一致させることができるリードのみが考慮される。それゆえ、完全に新規な遺伝子からメタゲノムリードを見つける既存の方法の正確さは未知であり、そして潜在的に低い。 ARGを検出するための他のアプローチは、感度を高めるために隠れマルコフモデルを使用するResfams[ref.32]、およびタンパク質構造についての情報を組み込むために機械学習を適用するPCM [ref.33]である。しかしながら、これらの方法は、より長く完全にアセンブリされた遺伝子配列について働くよう設計されている。最近、ショットガンのメタゲノムデータから直接新しいARGを見つける新しい方法deepARGが発表された。 deepARGは人工ニューラルネットワークをベースにした分類器を使用しているため、短いシーケンスフラグメントでも直接処理できる[ref.34]。しかしながら、deepARGは同定されたフラグメントを全長ARG配列にアセンブリする機能性を欠いている。したがって、まとめると、新規のARGからメタゲノムリードを同定し、次いでそれらの全配列を正確に再構築できるレジストソームの探索方法はない。

 ここではfARGeneというメタゲノムデータから直接ARGの同定と再構成を行う新しい方法を提示する。fARGeneは、既知のARGとの配列類似性が低い場合でも、以前には特徴付けられていなかった耐性遺伝子を正確に識別するために最適化された確率的遺伝子モデルを使用する。この方法は、最初にARGに由来する可能性があるフラグメントを同定することによって機能し、次いでこれは次に全長遺伝子に再構築される。この方法は全メタゲノムの完全なアセンブリを必要としないため計算上効率的であり、非常に大きなシーケンスデータセットにも適用可能である。この方法を実証するために、4つすべてのβ-ラクタマーゼアンブラークラス(A、B、C、およびD)を表すARGモデルを作成し、50億を超えるシーケンシングリード含む5つのメタゲノムデータセットに適用した。結果はARGを221の再構成し、そのうち58はこれまで報告されていなかった(NCBI  Genbank中のいずれの遺伝子に対しても70%以下の配列類似性)。さらに、以前にfARGeneによって再構築され大腸菌で発現された38の新規ARGのうち、80%超が機能的で大腸菌に耐性を付与した。最後に、本著者らはfARGeneがdeepARGおよび新規β-ラクタマーゼを検出する他の5つの方法と比較して有意により高い感度を有することを示す。 fARGeneはあらゆるクラスのARGに適用可能で、MITライセンスの下でGitHubhttps://github.com/fannyhb/fargene)から、ドキュメントとチュートリアルと共に、自由に利用できる。

 

tutorial

https://github.com/fannyhb/fargene/blob/master/tutorial/tutorial.md

 

インストール

依存

mamba create -n fARGene
conda activate fARGene
mamba install -c bioconda -y seqtk emboss hmmer spades trim-galore prodigal

#ORF_finderはbinaryをダウンロードし、実行権を付け、PATHの通ったディレクトリに移動する
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/TOOLS/ORFfinder/linux-i64/
chmod +x ORFfinder
mv /usr/local/bin/

本体 Github

(Gitubより)

setup.py will look for and try to install numpy and matplotlib so make sure that you either:

have these packages installed
run setup.py as root or with sudo
install the program in a conda environment.

git clone https://github.com/fannyhb/fargene.git
cd fargene
python setup.py install

fargene -h

# fargene -h

usage: fargene [-h] --infiles INFILES [INFILES ...] --hmm-model HMM_MODEL

               [--score LONG_SCORE] [--meta] [--meta-score META_SCORE]

               [--output OUTDIR] [--force] [--tmp-dir TMP_DIR] [--protein]

               [--processes PROCESSES] [--min-orf-length MIN_ORF_LENGTH]

               [--retrieve-whole] [--no-orf-predict] [--no-quality-filtering]

               [--no-assembly] [--orf-finder] [--store-peptides] [--rerun]

               [--amino-dir AMINO_DIR] [--fasta-dir FASTA_DIR]

               [--translation-format TRANS_FORMAT] [--loglevel {DEBUG,INFO}]

               [--logfile LOGFILE]

 

Searches and retrieves new and previously known genes from fragmented

metagenomic data and genomes. Copyright (c) Fanny Berglund 2018.

 

optional arguments:

  -h, --help            show this help message and exit

  --infiles INFILES [INFILES ...], -i INFILES [INFILES ...]

                        Input file(s) to be searched. Could either be in FASTA

                        or FASTQ format.

  --hmm-model HMM_MODEL

                        The Hidden Markov Model that should be used to analyse

                        the data. Could either be one of the pre-defined

                        models or the path to a custom HMM.

  --score LONG_SCORE, -sl LONG_SCORE

                        The threshold score for a sequence to be classified as

                        a (almost) complete gene (default: None).

  --meta                If the input data is paired end metagenomic data

                        (default: False).

  --meta-score META_SCORE, -sm META_SCORE

                        The threshold score for a fragment to be classified as

                        a positive. Expressed as score per amino acid

                        (default: None).

  --output OUTDIR, -o OUTDIR

                        The output directory for the whole run (default:

                        ./fargene_output).

  --force, -f           Overwrite output directory if it exists (default:

                        False).

  --tmp-dir TMP_DIR     Directory for (sometimes large) intermediate files.

                        (default: OUT_DIR/tmpdir)

  --protein             If the input sequence(s) is amino acids (default:

                        False).

  --processes PROCESSES, -p PROCESSES

                        Number of processes to be used when processing

                        metagenomic data (default: 1).

  --min-orf-length MIN_ORF_LENGTH

                        The minimal length for a retrieved predicted ORF (nt).

                        (default: 90% of the length of the chosen hmm.)

  --retrieve-whole      Use this flag if the whole sequence where a hit is

                        detected should be retrieved (default: False).

  --no-orf-predict      Do not perform ORF prediction.

  --no-quality-filtering

                        Use if no quality control should be performed on the

                        metagenomic data (default: False).

  --no-assembly         Use if you want to skip the assembly and retrieval of

                        contigs for metagenomic data (default: False).

  --orf-finder          Use NCBI ORFfinder instead of prodigal for ORF

                        prediction of genomes/contigs (default: False).

  --store-peptides, -sp

                        Store the translated sequences. Useful if you plan to

                        redo the analysis using a different model and want to

                        skip the preprocessing steps (default: False).

  --rerun               Use of you want to redo the analysis or do the

                        analysis using a different model and have kept either

                        the nucletide or amino acid sequences. Please note

                        that this only works if the input data is the same for

                        both runs (default: False).

  --amino-dir AMINO_DIR

                        Where the amino acid sequences generated by the method

                        are located. Only to be used in combination with

                        --rerun

  --fasta-dir FASTA_DIR

                        Where the nucleotide sequences in FASTA generated by

                        previous runs of the method are located. Only to be

                        used in combination with --rerun

  --translation-format TRANS_FORMAT

                        The translation format that transeq should use.

                        (default: pearson)

  --loglevel {DEBUG,INFO}

                        Set logging level (default: INFO).

  --logfile LOGFILE     Logfile (default: fargene_analysis.log).

導入に手間取ったのでdocker imageも上げておきます。他のツールのテストもしてたためライブラリをかなり入れていてサイズが肥大しています。使用する場合は注意して下さい。

docker pull kazumax/fargene
source ${HOME}/.profile
fargene -h

 

実行方法

 raw fastqを指定する。ランには

fargene -i input*.fastq --meta --hmm-model class_b_1_2 -o output 

(Gitubより) hmm-model can be any of the pre-defined models:

  • Class A beta-lactamases --hmm-model class_a
  • Subclass B1 and B2 beta-lactamases --hmm-model class_b_1_2
  • Subclass B3 beta-lactamases --hmm-model class_b_3
  • Class C beta-lactamases --hmm-model class_c
  • Class D beta-lactamases --hmm-model class_d_1 --hmm-model class_d_2
  • qnr --hmm-model qnr

 predicted-orfs*が再構築されたARGになる。周辺配列を調べるなら、retrieved-contigs.fastaを使う。

 

またはアセンブリして得たcontig配列(binnig済みのもの)を指定する。

fargene -i contig.fasta --hmm-model class_b_1_2 -o _output

 

protein配列も利用できる。

fargene -i contig_proteome.fasta --hmm-model class_b_1_2 --protein -o _output

 出力の詳細はGithubで説明されています。

 

Identification and reconstruction of novel antibiotic resistance genes from metagenomes

Fanny Berglund, Tobias Österlund, Fredrik Boulund, Nachiket P. Marathe, D. G. Joakim Larsso, Erik Kristiansson
Microbiome 2019 7:52

関連


 

 

 

ロングリードを使ってscaffoldsのgap closingを行うLR_Gapcloser

 

 次世代シークエンシング( NGS)技術は、デノボアセンブリによるゲノム配列の低コストおよび高速構築を可能にする。 NGS技術の利点と共に、この10年間で、多くのゲノムプロジェクト(例えば、10Kゲノムプロジェクト[ref.1]や100K病原体ゲノムプロジェクト[ref.2])が開始され、多数の種のゲノムがアセンブリされた[ref.3、4]。 。しかし、シーケンスバイアス[redf.5]、リピート領域[ref.6]、ヘテロクロマチン[ref.7]などの要因により、一部の領域はアセンブリが困難または不可能になり、ギャップや断片化されたゲノムアセンブリが発生する。

 ギャップクロージャのプロセスは、ゲノムアセンブリの完全性と連続性を高めるための最後の、しかし最も重要なステップである。完全なゲノムを得るために、GapFiller [ref.8]、GapCloser [ref.9]、Sealer [ref.10]、GapBlaster [ref.11]、GapReduce [ref.12]、およびGap2Seq [ref.13]を含むいくつかのギャップクロージャアプローチがNGSリードまたは事前にアセンブリされたコンティグ[ref.14]の隙間を埋めるために利用される。しかしながら、これらのツールは、ギャップクロージャプロセス中に高い割合のミスアセンブリを起こす[ref.15]。さらに、これらのツールを使用してすべてのギャップ、特に大きなギャップを埋めることは困難である。例えばPacific Biosciences(PacBio)やNanoporeプラットフォームなどの第3世代シーケンシング(TGS)技術としても知られる長い1分子シーケンシング技術は、これらのギャップを埋める可能性を秘めた長い偏りのないリードを生み出し、完全なゲノムアセンブリを達成する。 PBJelly [ref.18]とGMcloser [ref.15]は、ギャップを埋めるためにPacBioリードを採用している。 PBJellyは、基本的なローカルアラインメント[ref.19]を使用してロングリードをリファレンスアセンブリにアラインメントさせ、サポートするリードを選択し、ローカルギャップアセンブリを実行し、ギャップを埋め正確なアセンブリを決定する。 GMcloserは、scaffoldsをサブコンティグに分割し、MUMmer [ref.20]またはblastn(basicローカルアライメント検索ツール)を使用してロングリードをサブコンティグに位置合わせする。 GMcloserは、尤度ベースの分類子を使用して、ロングリードをscaffoldsのギャップに正しく割り当てる。ただし、実行時間が長い、メモリ使用量が多い、クロージャパフォーマンスが低いなどの欠点があるため、特に大規模で複雑なゲノムではアプリケーションが制限される。したがって、ゲノムアセンブリのギャップを埋めるために、高速でメモリ効率の良いギャップクロージャアプローチが必要である。

 ここでは、ロングリードを使用してアセンブリ内のギャップを効率的かつ迅速に埋めるLR_Gapcloserを開発した。以前のギャップ解消ツールと比較して、多くの注目すべき利点があり、これにはより高いギャップクロージャー性能、より少ない実行時間、より少ないピークメモリー、そしてより少ないミスアセンブリが含まれる。大きくて複雑なゲノムおよび繰り返し由来のギャップの両方に対してさえ、このツールはより良い性能を示した。最後に、NGSとTGS技術の両方を使用してシーケンスされたゲノムについて、異なるハイブリッドアセンブリ戦略の連続性と正確性を評価し、TGSベースとNGSベースのアセンブリをLR_Gapcloserと組み合わせて高品質のアセンブリを作成する最適ハイブリッド戦略を提案した。

 

 

f:id:kazumaxneo:20190510222450p:plain

The primary steps in LR_Gapcloser.  論文より転載

  

インストール 

ubuntu16.04でテストした。

依存

  • Perl and Bioperl should be installed on the system.
  • GLIBC 2.14 should be installed.

git clone https://github.com/CAFS-bioinformatics/LR_Gapcloser.git
cd LR_Gapcloser/src/

> ./LR_Gapcloser.sh

$ ./LR_Gapcloser.sh 

ls: /proc/81616/fd/: No such file or directory

Usage:sh LR_Gapcloser.sh -i Scaffold_file -l Corrected-PacBio-read_file

  -i    the scaffold file that contains gaps, represented by a string of N          [         required ]

  -l    the raw and error-corrected long reads used to close gaps. The file should                      

        be fasta format.                                                               [         required ]

  -s    sequencing platform: pacbio [p] or nanopore [n]                             [ default:       p ]

  -t    number of threads (for machines with multiple processors), used in the bwa                      

        mem alignment processes and the following coverage filteration.             [ default:       5 ]

  -c    the coverage threshold to select high-quality alignments                    [ default:     0.8 ]

  -a    the deviation between gap length and filled sequence length                 [ default:     0.2 ]

  -m    to select the reliable tags for gap-closure, the maximal allowed                                

        distance from alignment region to gap boundary (bp)                         [ default:     600 ]

  -n    the number of files that all tags were divided into                         [ default:       5 ]

  -g    the length of tags that a long read would be divided into (bp)              [ default:     300 ]

  -v    the minimal tag alignment length around each boundary of a gap (bp)         [ default:     300 ]

  -r    number of iteration                                                         [ default:       3 ]

  -o    name of output directory                                                    [ default: ./Output]

 

 

実行方法

pacbio

bash LR_Gapcloser.sh -i input_scaffolds.fasta -l Corrected-PacBio-reads.fasta -s p
  •  -l    the raw and error-corrected long reads used to close gaps. The file should be fasta format
  • -i     the scaffold file that contains gaps, represented by a string of N [required]
  • -s    sequencing platform: pacbio [p] or nanopore [n]  [ default: p ]

 

nanopore

bash LR_Gapcloser.sh -i input_scaffolds.fasta -n nanopore.fasta -s n
  •  -l    the raw and error-corrected long reads used to close gaps. The file should be fasta format
  • -i     the scaffold file that contains gaps, represented by a string of N [required]
  • -s    sequencing platform: pacbio [p] or nanopore [n]  [ default: p ]

 

引用

LR_Gapcloser: a tiling path-based gap closer that uses long reads to complete genome assembly

Gui-Cai Xu Tian-Jun Xu Rui Zhu Yan Zhang Shang-Qi Li Hong-Wei Wang Jiong-Tang Li

GigaScience, Volume 8, Issue 1, January 2019, giy157

 

この論文でも使用することが推奨されています。


 

メタゲノムのアセンブリ配列からプラスミド配列を予測する PlasFlow

 

 プラスミドは、変化する環境条件下で急速な進化とそれらの宿主の適応を促進するmobile genetic elementsである(ref1,2)。プラスミドは、宿主細胞内で自律的に複製するの染色体外のDNA断片であり、細菌種において広く存在している。既知のプラスミドの大部分は環状DNAとして存在し、これはアルカリlysis法を用いた容易な単離を可能にする重要な特徴である。しかしながら、主にBorelliaStreptomycesNocardiaおよびRhodococcus属からのプラスミドもあり、それらは線状である(ref.3)。プラスミドの重要な特徴は、それらが通常一次代謝プロセスに一般的に割り当てられる遺伝子は欠いているが、むしろ宿主の環境適合性を改善するか異化または抵抗機能をコードする遺伝子を保有することである(ref.4–6)。さらに、それらは多様な分類群からの異なる種間の遺伝子水平伝播に寄与することができ、これはプラスミドを重大な生態学的影響を与える因子にする(ref.7)。従って、プラスミド指向の研究は、多様な環境で起こる過程をよりよく理解するために重要であり、そして新しいプラスミドを同定しそしてシーケンシングするための方法が必要とされている。

 特定の表現型に基づいた宿主細菌の分析中に偶然に多数のプラスミドが同定されている(ref.8)。しかしながら、この同定方法は非常に面倒であり、培養とは無関係に特定の環境試料の全プラスミドDNA含有量、いわゆる「plasmidome」への適切な洞察を提供することができない(ref.3)。培養依存性plasmidome研究(ref.9-12)は、特定の細菌または細菌群のmobile genetic elementsの理解に大きく貢献しているが、これらの方法では培養不可能な生物からのプラスミドは手が届かない。Exogenous plasmid isolation(ref.13、14)およびTransposon-aided Capture(TRACA、(ref.6、15〜16))を含むこの問題に取り組むいくつかのアプローチが開発されているが、これらの方法には重大な制限がある。低処理量の外性プラスミド単離のためには、捕捉されるプラスミドは接合性(または少なくとも可動性)であり、レシピエント細胞内で安定的に複製される必要がある(ref.3)。これにより、移動できないプラスミドのヒット数が大幅に減る。TRACAはよりハイスループットだが、サイズが2〜10 kbの小さなプラスミドしか捕捉できない(ref.3,17)。環境からプラスミドを直接単離することは可能だが、それは一般に小さいものに限定されている(ref.17)。以前に環境サンプルからの間接的プラスミド単離とそれに続くエキソヌクレアーゼ処理およびPhi29増幅を使用する技術が開発されている(ref.18–20)が、この方法は高いバクテリアバイオマス含有量および酵素操作を妨害する汚染物質が少ないサンプルに有用である。

 メタゲノムデータセット中のプラスミド配列を同定するためのバイオインフォマティック方法も開発されているが、分子方法と同様に、それらは一般に環状エレメントの同定を目的としている(ref.21)。シーケンシング技術の迅速な発展及びシーケンシングコストの削減は、メタゲノムプロジェクト数の増加をもたらし、その結果、シーケンスデータベースは絶えず成長してきた。与えられたメタゲノムはほとんど染色体とプラスミドのmixtureであり、両者の関係は染色体にとって非常に有利であり続けているので(ref.3)、多くのプラスミド配列がシーケンシングされたメタゲノムに含まれているが未確認のままである。メタゲノムデータからプラスミドをアセンブリする試みがいくつかなされている(ref.6,22)が、それらは面倒で計算集約的なアプローチに頼っている。最近、de Bruijnアセンブリグラフから環状コンティグを回収することを目的とした、Recyclerと同様に、プラスミドをアセンブリすることができるSPAdesアセンブラバージョンが開発された(ref.23)。もう1つの簡単なアプローチはPlasmidFinder(ref25)に実装されている。これは、細菌ゲノムの集まりの中のプラスミドレプリコンの同定を目的としたWebベースのプログラムである。このプログラムは、可能性のあるプラスミドを同定するため、明確に定義されたレプリコン配列に対して類似性検索を行う。これはプラスミド由来の配列を高精度で予測することができるが、その主な制限はデータベースのサイズにあり、それは大部分が腸内細菌科レプリコンからなる。この側面はメタゲノム研究のためのその使用法を著しく制限している。データベース類似性検索アプローチはPlasmid Constellation Network (PLACNET)において拡張され、これはリファレンスデータベースに対して配列を比較するためにBLASTを使用し、次いで再構築したプラスミドのネットワーク分析を行う。しかしながら、このプログラムは得られたシーケンスクラスターの手動キュレーションに依存しているので、いかなる自動アノテーションパイプラインでの使用も妨げられている。さらに、PLACNETから得られた結果は完全に再現可能ではなく、研究者の専門知識に依存している。

 ショットガンメタゲノムデータからプラスミドを同定する別の興味深いアプローチは、任意のメタゲノムにおけるプラスミド配列の自動検出を可能にする機械学習技術の使用である。ゲノムシグネチャに基づく方法は、プラスミド - 宿主類似性がゲノム%G + C含有量と相関することを明らかにした(ref.27, 28)。プラスミドゲノムのバーコードもまた、おそらく細胞培養物間でのプラスミドの頻繁な移動により引き起こされる類似の選択圧のために、類似の特徴を有する傾向がある(ref.29)。ゲノムシグネチャに基づいてプラスミド配列を同定する方法は、cBarソフトウェア(ref.30)で実装された。

 最近、シーケンシングデータからのプラスミド再構築を目的とした利用可能なツールの性能が比較され(ref.31)、特に大型(>50kb)プラスミドには依然として限界があることが示された。しかしながら、そのレビューは複雑なメタゲノムプロジェクトのシーケンシングデータよりむしろ単一のバクテリア単離株に焦点を合わせていた。その結果、メタゲノムデータセットにおけるプラスミド同定に向けられたアルゴリズムの系統的な比較が欠如している。本研究では、メタゲノムコンティグ内の細菌プラスミド配列を予測する新しいアプローチPlasFlowを紹介し、それを他の利用可能なツールと比較する。 9,565の細菌の染色体とプラスミドからのシーケンスのゲノムシグネチャを使用して、異なる系統から染色体とプラスミドのシーケンスを分離するためにディープニューラルネットワークモデルを訓練した。本著者らのアプローチは、試験データ上のプラスミド予測において96%の分類精度を達成し、これは他の以前に開発されたツールよりも著しく優れている。実際のメタゲノムデータセットに対して実施された試験は、メタゲノムデータを用いたplasmidome分析に対するその多用性を明らかにした。

 

f:id:kazumaxneo:20190510174626j:plain

Flowchart describing the training and classification procedures implemented in the PlasFlow. 論文より転載

 


インストール

オーサーの指示に従い、condaでpython3.5の仮想環境を構築してテストした(mac os 10.12使用)。

依存

  • Python 3.5
  • Python packages:
  • Scikit-learn 0.18.1
  • Numpy
  • Pandas
  • TensorFlow 0.10.0
  • rpy2 >= 2.8
  • scipy
  • biopython
  • dateutil >= 2.5

Conda is recommended option for installation as it properly resolve all dependencies (including R and Biostrings) and allows for installation without messing with other packages installed.

Githubより)

本体  github

conda config --add channels bioconda
conda config --add channels conda-forge

conda create --name plasflow python=3.5
source activate plasflow
#Mac users should install Tensorflow at this step
conda install -c jjhelmus tensorflow=0.10.0rc0
#then,
conda install -c smaegol plasflow

#Perl scripts included with PlasFlow requires few Perl modules.
conda install -c bioconda perl perl-bioperl perl-getopt-long

#To deactivate an active environment, use:
source deactivate

もしくはpipでも導入できる。Github参照。 

PlasFlow.py -h

$ PlasFlow.py -h

usage: PlasFlow.py [-h] --input INPUTFILE --output OUTPUTFILE

                   [--threshold THRESHOLD] [--labels LABELS] [--models MODELS]

                   [--batch_size BATCH_SIZE]

 

PlasFlow v.1.1 - predicting plasmid sequences in metagenomic data using genome

signatures. PlasFlow is based on the TensorFlow artificial neural network

classifier

 

optional arguments:

  -h, --help            show this help message and exit

  --input INPUTFILE     Input fasta file with sequences to classify (required)

  --output OUTPUTFILE   Output file with classification results (required)

  --threshold THRESHOLD

                        Threshold for probability filtering (default=0.7)

  --labels LABELS       Custom labels file

  --models MODELS       Custom models localization

  --batch_size BATCH_SIZE

                        Batch size for large datasets

 

> filter_sequences_by_length.pl  

$ filter_sequences_by_length.pl  

Filters sequences in multifasta file by length

 

Usage of program:

 

 /Users/user/.pyenv/versions/miniconda3-4.3.21/envs/plasflow/bin/filter_sequences_by_length.pl [options] 

-input multifasta file with sequences to be filtered

-output output file for filtered sequences 

-thresh sequence length threshold [bp]

 

テストラン

 必要に応じてサイズでフィルタリングしてから使用する。1000bp以上のみ抽出。

git clone https://github.com/smaegol/PlasFlow.git
cd PlasFlow/test/

filter_sequences_by_length.pl \
-input Citrobacter_freundii_strain_CAV1321_scaffolds.fasta \
-output filtered_output.fasta -thresh 1000

 

PlasFlowを実行する。 scaffolds.fastaを指定する。

PlasFlow.py --input Citrobacter_freundii_strain_CAV1321_scaffolds.fasta --output test.plasflow_predictions.tsv --threshold 0.7

$ PlasFlow.py --input Citrobacter_freundii_strain_CAV1321_scaffolds.fasta  --output test.plasflow_predictions.tsv --threshold 0.7

Importing sequences

Imported  78  sequences

Succesfully read previously calculated kmer frequencies for kmer 5

Predicting labels using kmer 5  frequencies

Succesfully read previously calculated kmer frequencies for kmer 6

Predicting labels using kmer 6  frequencies

Succesfully read previously calculated kmer frequencies for kmer 7

Predicting labels using kmer 7  frequencies

Voting...

Filtering by probability threshold 0.7

 

Resulting plasmid sequences prediction:

   chromosome  plasmid  unclassified

0          34       33            11

 

Resulting taxonomical assignment:

   unclassified  Proteobacteria

0            10              68

 

出力

f:id:kazumaxneo:20190505230742p:plain

test.plasflow_predictions.tsv

f:id:kazumaxneo:20190505230645p:plain

 

最近発表されたmetaplasmidSPAdesと組み合わせると良いかもしれませんね。

 

引用

PlasFlow: predicting plasmid sequences in metagenomic data using genome signatures

Krawczyk PS, Lipinski L, Dziembowski A

Nucleic Acids Res. 2018 Apr 6;46(6):e35

 

関連


 

 

 

 

 

 

Illumina、454、およびPacBioのSmith-Watermanアライメントによる高感度なアライナー InDelFixer

 2019 5/30 インストール追記

 

InDelFixerは454、Illumina、およびPacBioデータ用の高感度なアライナーである。完全なSmith-Watermanアライメントを採用している。事前の高速k-merマッチングによって次世代シーケンス(NGS)および第3世代のリードを一連のリファレンスシーケンスにアライメントする。

 

インストール

mac os10.14のjava1.8環境でjarファイルをダウンロードしてテストした。

Github

リリースから.jarファイルをダウンロードする。condaでも導入できる。

#bioconda (link)
conda install -c bioconda indelfixer

java -jar InDelFixer.jar

$ java -jar InDelFixer.jar 

 

InDelFixer version: 1.1

Get latest version from http://bit.ly/indelfixer

 

USAGE: java -jar InDelFixer.jar options...

 

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

 === GENERAL options ===

  -o PATH : Path to the output directory (default: current directory).

  -i PATH : Path to the NGS input file (FASTA, FASTQ or SFF format) [REQUIRED].

  -ir PATH : Path to the second paired end file (FASTQ) [ONLY REQUIRED if first file is also fastq].

  -g PATH : Path to the reference genomes file (FASTA format) [REQUIRED].

  -r interval : Region on the reference genome (i.e. 342-944).

  -k INT : Kmer size (default 10).

  -v INT : Kmer offset (default 2).

  -cut INT : Cut given number of bases (primer) at 5' and 3' (default 0).

  -refine INT : Computes a consensus sequence from alignment and re-aligns against that.

  Refinement is repeated as many times as specified.

  -mcc INT : Minimal coverage to replace a reference base in the consensus (default 1).

  -rmDel : Removes conserved gaps from consensus sequence during refinement.

  -sensitive : More sensitive but slower alignment.

  -fix : Fill frame-shift causing deletions with consensus sequence.

  -noHashing : No fast kmer-matching to find approximate mapping region. Please use with PacBio data!

  -realign DOUBLE : Reads are aligned to the whole reference sequence,

  if the relative mismatch rate is above the given threshold (default 0.1).

 

 === FILTER ===

  -l INT : Minimal read-length prior alignment (default 0).

  -la INT : Minimal read-length after alignment (default 0).

  -ins DOUBLE : The maximum percentage of insertions allowed [range 0.0 - 1.0] (default 1.0).

  -del DOUBLE : The maximum percentage of deletions allowed [range 0.0 - 1.0] (default 1.0).

  -sub DOUBLE : The maximum percentage of substitutions allowed [range 0.0 - 1.0] (default 0.5).

  -maxDel INT : The maximum number of consecutive deletions allowed (default no filtering).

  -q INT : Minimal average Phred score of the aligned read (default 20).

 

 === GAP costs ===

  -gop : Gap opening costs for Smith-Waterman (default 30).

  -gex : Gap extension costs for Smith-Waterman (default 3).

 

 === GAP costs predefined ===

  -454 : 10 open / 1 extend

  -illumina : 30 open / 3 extend

  -pacbio : 5 open / 3 extend

 

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

 === EXAMPLES ===

  454/Roche : java -jar InDelFixer.jar -i libCase102.fastq -g referenceGenomes.fasta -454

  PacBio : java -jar InDelFixer.jar -i libCase102.ccs.fastq -g referenceGenomes.fasta -noHashing -pacbio

  Illumina : java -jar InDelFixer.jar -i libCase102_R1.fastq -ir libCase102_R2.fastq -g referenceGenomes.fasta -illumina

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

 

 

インストール

ペアエンドfastq

java -jar InDelFixer.jar -i pair_R1.fastq -ir pair_R2.fastq -g ref.fasta 
  • -o     Path to the output directory (default: current directory).
  • -i      Path to the NGS input file (FASTA, FASTQ or SFF format) [REQUIRED].
  • -ir     Path to the second paired end file (FASTQ) [ONLY REQUIRED if first file is also fastq].
  • -g     Path to the reference genomes file (FASTA format) [REQUIRED].

 

pacbioのccs.fasta

java -jar InDelFixer.jar -i libCase102.fasta -g ref.fasta -noHashing

 

 

引用

https://github.com/cbg-ethz/InDelFixer

 

関連

kazumaxneo.hatenablog.com

 

 

メタゲノムのraw fastqから高速なtaxonomy assignmentを行う FOCUS

 

 微生物は他のどの細胞生物よりも豊富であり(Whitman、Coleman&Wiebe、1998年)、どの生物が存在し、それらが何をしているのかを理解することが重要である(Handelsman、2004)。多くの環境では、微生物群集の大多数は培養できず、メタゲノムは未培養のゲノムを直接探索し、それらのDNAのみを用いて微生物群集の多様性を理解するための強力なツールである(Sharon&Banfield、2013)。

 微生物群集を理解することは、生物学の多くの分野で重要である。例えば、メタゲノムは、海洋動物(Trindade-Silva et al、2012)または病状(Belda-Ferre et al、2012)に関連する微生物の分類学的および機能的シグネチャを区別することができる。大量のシーケンシング、リード長の短さ、およびシーケンシングエラーは、メタゲノムに存在する生物の多様性を同定することを困難にしている(Mande、Mohammed&Ghosh、2012)。このために多くのプログラムが存在し、それらは相同性または構成に基づいている。

 相同性に基づくプログラムは、通常、BLASTプログラムを使用して(Altschul et al、1997)、大規模データベースの出力における最も良いヒットを識別する。MG-RAST(Meyer et al、2008)では、メタゲノム試料を分類するために、配列を一組のデータベースにアライメントさせる。 MetaPhlAn(Segata et al、2012)およびGenomePeek(K McNair、R Edwards、未発表データ(published in 2015))は、BLAST検索を高速にすることを可能にする、マーカー遺伝子、例えばユニーククレードおよびハウスキーピング遺伝子のみを含む簡約データベースを使用する。 PhymmBL(Brady&Salzberg、2011)は、補間マルコフモデルを用いてBLASTの結果を改善する。 GASiC(Lindner&Renard、2013)は、Bowtie(Langmead et al、2009)とリファレンスゲノムの類似性を用いて、推定された存在量を補正する。 GPUを必要とする高速プログラムであるParallel-Meta(Su、Xu&Ning、2012)は、相同性の結果を改善するためにmegaBLAST(Zhang et al、2000)およびHMM(Hidden Markov Model)を使用している。これらのアプリケーションのほとんどは、シーケンスを個別に分類し、ビンを合計することによって分類プロファイルを生成する。

 一般に、組成ベースのアプローチは、オリゴヌクレオチド頻度を使用する。 Taxy(Meinicke、Aßhauer&Lingner、2011)は、メタゲノム中およびリファレンスゲノム中のオリゴヌクレオチド分布を使用し、そしてメタゲノム中に存在する生物を同定するために混合モデリングを使用し、そしてRAIphy(Nalbantogluら、2011)は、オリゴヌクレオチドおよび相対存在比指数を使用して生物を同定する。 

 メタゲノム全体のk-mer構成を使用してtaxonomy プロファイルを再構築する新しいアプローチを開発した。メタゲノムのk-mer組成をリファレンスデータベースの生物と一致させるために、非負最小二乗法(NNLS)を使用して生物量の最適セットを計算する。未知の配列をクラスター化するためにk-mers​​が以前に使用され(Teeling et al、2004;McHardy et al、2007)、NNLSが遺伝子数の変動に基づいてメタゲノム試料中に存在する属を同定するために使用された。Orr & Borenstein、2013)。ここでは、メタゲノムに存在する分類群を識別するための、超高速で正確な合成ベースのアプローチであるFOCUSでこれら2つのアプローチを組み合わせる。 FOCUSとGASiC、MetaPhlAn、RAIphy、PhymmBL、Taxy、MG-RASTのパフォーマンスを比較した。 

 

 

f:id:kazumaxneo:20190507160451p:plain

Workflow of the FOCUS program. 

 

 

インストール 

本体 Github 

#anaconda環境
conda install -c bioconda -y focus

#pip3
pip3 install metagenomics-focus

focus -h

/# focus -h

usage: focus [-h] -q QUERY -o OUTPUT_DIRECTORY [-k KMER_SIZE]

             [-b ALTERNATE_DIRECTORY] [-p OUTPUT_PREFIX] [-t THREADS]

 

FOCUS: An Agile Profiler for Metagenomic Data

 

optional arguments:

  -h, --help            show this help message and exit

  -q QUERY, --query QUERY

                        Path to directory with FAST(A/Q) files

  -o OUTPUT_DIRECTORY, --output_directory OUTPUT_DIRECTORY

                        Path to output files

  -k KMER_SIZE, --kmer_size KMER_SIZE

                        K-mer size (6 or 7) (Default: 6)

  -b ALTERNATE_DIRECTORY, --alternate_directory ALTERNATE_DIRECTORY

                        Alternate directory for your databases

  -p OUTPUT_PREFIX, --output_prefix OUTPUT_PREFIX

                        Output prefix (Default: output)

  -t THREADS, --threads THREADS

                        Number Threads used in the k-mer counting (Default: 4)

 

example > focus -q samples directory

 

データベースの準備

git clone https://github.com/metageni/FOCUS.git

#データベースdb.zipを解凍しFOCUS/focus_app/に置く。
unzip FOCUS/focus_app/db.zip
mv db FOCUS/focus_app/

 

実行方法

ランの際はfastq/fastaディレクトリと、解凍したdbフォルダを指定する(dbの手前のパスまで指定する)。

focus -q test/ -o output -b FOCUS-master/focus_app/ -t 8
  • -o      Output_dir
  • -b      Alternate directory for your databases
  • -q      Path to directory with FAST(A/Q) files
  • -t       Number Threads used in the k-mer counting (Default: 4)

gzip圧縮したfastqには対応しないことに注意する。また、input_dir/のfastqの拡張子は.fastqにする。

f:id:kazumaxneo:20190508220704p:plain

output_All_levels.xls

f:id:kazumaxneo:20190508220904p:plain

ディレクトリに複数データがある場合、データごとに定量され出力される。ここでは4つのfastqの結果がまとめて出力されている。

 

SUPER-FOCUSも近い内に紹介します。

引用

FOCUS: an alignment-free model to identify organisms in metagenomes using non-negative least squares

Silva GG1, Cuevas DA1, Dutilh BE2, Edwards RA

PeerJ. 2014 Jun 5;2:e425. doi: 10.7717/peerj.425. eCollection 2014.

Virusの ultra deep NGSのbamからコンセンサス配列を出力する ConsensusFixer

2019 5/14リンク追加

2019 5/30 インストール追記

 

ConsensusFixerはjavaコマンドラインアプリケーション。virusのウルトラディープNGSのアライメント(インフレーム挿入とあいまいなヌクレオチドを含む)からコンセンサスシーケンスを計算、出力する。European Virus Bioinformatics Center のツール一覧(リンク)で紹介されている。

 

インストール

依存

本体 Github

リリースから.jarファイルをダウンロードするかcondaで導入する。

#bioconda
conda install -c bioconda -y consensusfixer

java -jar ConsensusFixer.jar

$ java -jar ConsensusFixer.jar 

No input given

 

ConsensusFixer version: 0.4

 

USAGE: java -jar ConsensusFixer.jar options...

 

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

 === GENERAL options ===

  -i INPUT : Alignment file in BAM format (required).

  -r INPUT : Reference file in FASTA format (optional).

  -o PATH : Path to the output directory (default: current directory).

  -mcc INT : Minimal coverage to call consensus.

  -mic INT : Minimal coverage to call insertion.

  -plurality DOUBLE : Minimal relative position-wise base occurence to integrate into wobble (default: 0.05).

  -pluralityN DOUBLE : Minimal relative position-wise gap occurence call N (default: 0.5).

  -m : Majority vote respecting pluralityN first, otherwise allow wobbles.

  -f : Only allow in frame insertions in the consensus.

  -d : Remove gaps if they are >= pluralityN.

  -mi : Only the insertion with the maximum frequency greater than mic is incorporated.

  -pi : Progressive insertion mode, respecting mic.

  -pis INT : Window size for progressive insertion mode (default: 300).

  -dash : Use '-' instead of bases from the reference.

  -s : Single core mode with low memory footprint.

 

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

 === Technical options ===

  -XX:NewRatio=9 : Reduces the memory consumption (RECOMMENDED to use).

  -Xms2G -Xmx10G : Increase heap space.

  -XX:+UseParallelGC : Enhances performance on multicore systems.

  -XX:+UseNUMA : Enhances performance on multi-CPU systems.

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

 === EXAMPLES ===

   java -XX:+UseParallelGC -Xms2g -Xmx10g -XX:+UseNUMA -XX:NewRatio=9 -jar ConsensusFixer.jar -i alignment.bam -r reference.fasta

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

 

実行方法

bamとリファレンスのfastaを指定する。リファレンスのfastaにchrが2つ以上含まれているとエラーになる。

java -jar ConsensusFixer.jar -i input.bam -r reference.fasta -mcc 0 -o output
  • -i      Alignment file in BAM format (required).
  • -r     Reference file in FASTA format (optional).
  • -o    Path to the output directory (default: current directory).
  • -mcc   Minimal coverage to call consensus.

outputconsensus.fastaが出力される。

 

追記

samtoolsを使った方法は

sam/bamファイルを変換、編集したり分析するためのツール - macでインフォマティクス の下の方にある"Consensus fasta"のコマンドを参照。

 

引用

Github

https://github.com/cbg-ethz/ConsensusFixer

 

参考

Question: Generate consensus sequence from BAM without ambiguity codes

https://www.biostars.org/p/302736/