macでインフォマティクス

macでインフォマティクス

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

構造変化の検出ツール RAPTR-SV

 

RAPTER-SVは構造変化を検出するjavaのプログラム。split-readとreid-pairのアプローチで構造変化を予測する。 タンデムデュプリケーションの 検出などにおいて競合するツールより感度が高いとされている。

 

Github

https://github.com/njdbickhart/RAPTR-SV

インストールマニュアル

https://github.com/njdbickhart/RAPTR-SV/wiki/installation

example

https://github.com/njdbickhart/RAPTR-SV/wiki/example

 

インストール 

依存

  • mrsFAST

mrsfast-ultraのインストール(リンク

git clone https://github.com/sfu-compbio/mrsfast 
cd mrsfast
make

本体のダウンロード(java

https://github.com/njdbickhart/RAPTR-SV/releases

 

 ラン

まずbwaでアライメントしてソートしたbamを作る。

 

bwa index reference.fa
samtools faidx reference.fa
Picard CreateSequenceDictionary REFERENCE=reference.fa OUTPUT=reference.dict
bwa mem -t 12 -M -R '@RG\tID:seq1.lib.run\tLB:seq1.lib\tPL:ILLUMINA\tSM:seq1' reference.fa sequence_1.fq sequence_2.fq > sample.sam
samtools view -bS sample.sam | samtools sort -@ 12 - -o sample.bam
samtools index sample.bam

 

biological replicateがあるなら、samtools mergeで1つのbamにマージしておく。

samtools index merge.bam replicates1.bam replicates2.bam replicates3.bam -@ 12

(merge: samtools 1.6 or higher)

 

ここではマージするbamがないものとして書いてます。

 duplicate readsにタグをつける(markduplicates)。

Picard MarkDuplicates INPUT=sample.bam OUTPUT=sample.sorted.nodup.bam METRICS_FILE=sample.dup.metrics VALIDATION_STRINGENCY=LENIENT
samtools index sample.sorted.nodup.bam

 

SVを検出する。

java -jar RAPTR-SV.jar preprocess -i sample.sorted.nodup.bam -o sample.preprocess -r reference.fa
java -jar RAPTR-SV.jar cluster -s sample.preprocess.flat -o sample.cluster

 

 

 

 

引用

RAPTR-SV: a hybrid method for the detection of structural variants.

Bickhart DM1, Hutchison JL1, Xu L1, Schnabel RD1, Taylor JF1, Reecy JM1, Schroeder S1, Van Tassell CP1, Sonstegard TS1, Liu GE1.

Bioinformatics. 2015 Jul 1;31(13):2084-90. doi: 10.1093/bioinformatics/btv086. Epub 2015 Feb 16.

 

SNPsをエラーとして扱わないマッピングが可能な mrsFAST-Ultra

 

mrsFAST-UltraはSNPsに対応した次世代リードのアライメントツール。 mrsFASTの改良版となる。既知SNPsを許容しながら(ミスマッチとして扱わないためidentityが上がる)アライメントを行うことができる。indexファイルの軽量化にも成功しており、bowtie2でindexをつけるとヒトゲノムでは4GB近くの容量になってしまうが、mrsFAST-Ultraではおよそ2GBのindexで済む。また冗長なアライメントについて多様なオプションを持っており、bowtie2より高感度かつ高速(6倍)にアライメントを実行することができる。

  

 

公式ページ

mrsFAST: micro-read substitution-only Fast Alignment Search Tool

マニュアル

https://github.com/sfu-compbio/mrsfast//blob/master/README.md

 

インストール

Github

https://github.com/sfu-compbio/mrsfast

git clone https://github.com/sfu-compbio/mrsfast 
cd mrsfast
make
./mrsfast -h #ヘルプが表示されるか確認

 mrsfastをパスの通ったディレクトリにコピーしておく。

 

 

ラン

サンプルデータをランする。

公式サイトからサンプルデータをダウンロード。解凍して中に入る。

http://sfu-compbio.github.io/mrsfast/

 

まずリファレンスのindexを作成する。

cd sfu-compbio-mrsfast-390c8eb/
mrsfast --index genome.fa

indexファイル(バイナリ)ができる。 

 

アライメントを行う。シングルリード。

mrsfast --search genome.fa --seq reads.fq -o mapped.sam -u unmapped.fq --threads 8 -e 4
  • -o Output the mapping record into file (default: output.sam)
  • -u Output unmapped reads in file (default: output.nohit).
  • --seq Set the input sequence to file. In paired-end mode, file should be used if the read sequences are interleaved. This file will be generated in all mapping mode.
  • --threads 8 Use 1 number of cores for mapping the sequences (default: 1). Use 0 to use all the available cores in the system.
  • -e error-threshold Allow up to error-threshold mismatches in the mappings.

デフォルトでは96 %同一の配列があればアライメントされる(-e 4)。

 

ペアリード。

mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8
  • --seq1 Set the input sequence (left mate) to file. Paired-end option only. 
  • --seq2 Set the input sequence (right mate) to file. Paired-end option only.
  • --pe Map the reads in Paired-End mode. min and max of template length will be calculated if not provided by corresponding options.
  • --threads 1 Use 1 number of cores for mapping the sequences (default: 1). Use 0 to use all the available cores in the system.

インターレースのペアリードを指定するなら--seqを使う。

 

ペアリードのインサートサイズ(リードの外側からの距離)でアライメントを制限する。例えばインサートサイズが100以上3000以下に制限するなら

mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8 --min 100 --max 3000

100以下や3000異常など異常なインサートサイズの可能性が高いペアリード(discordant read pair)を排除できる。

 

minmax--discordant-vhとともに使うと、正常なペアリードの範囲を指定するフラグとして働く。

mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8 --min 100 --max 1000 --discordant-vh
  •  --discordant-vh Map the reads in discordant fashion that can be processed by Variation Hunter / Common Law. Output will be generate in DIVET format.
  • --max-discordant-cutoff m Allow m discordant mappings per read. Should be only used with --discordant-vh option.

 

ペアリードそれぞれについて先頭の40bpだけを使いアライメントを行う。

mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8 --crop 40
  • --crop n Trim the reads to n base pairs from begining of the read.

5'末端 n-bpまでトリムではなく、3'末端を例えば50-bpをトリムする(捨てる)場合は--tail-crop 50とつける。

 

 

--その他--

  • リファレンスにアライメントされたリードはsam形式で出力される。リファンレンスにアライメントされなかったリードはfastqのまま出力される。
  • unmappedのfastqが必要なければ--disable-nohitsをつけることでunmapの出力を無効化できる。
  • samのヘッダーが必要なければ--disable-sam-headerをつける。
  • 巨大なfastqをアライメントする際は、--mem 4 などとつけることで、使用メモリー量を制限できる。この例では4GB以内でリードを読めるだけキャッシュし、それを繰り返すことでアライメントが進行する。
  • gz圧縮されたfastqを入力とする時は、--seqcompをつける。
  • 出力をgz圧縮したければ--outcompをつける。
  • 1リードが複数アライメントされる場合、-nオプションでアライメントされる数を制限できる。-n 5なら最大で5カ所となる。
  • --bestをつけるとベストマッピング(hamming distanceが最小)の箇所のみアライメントされる。

 

 他にもdbSNPのデータベースを使い、既知SNPsでアライメントが正しく行えるようにするフラグも用意されています。その辺りは公式サイトを参考にしてください。

mrsfast/README.md at master · sfu-compbio/mrsfast · GitHub

 

 

 

引用

mrsFAST-Ultra: a compact, SNP-aware mapper for high performance sequencing applications

Faraz Hach Iman Sarrafi Farhad Hormozdiari Can Alkan Evan E. Eichler S. Cenk Sahinalp

Nucleic Acids Research, Volume 42, Issue W1, 1 July 2014, Pages W494–W500, https://doi.org/10.1093/nar/gku370

 

mrsFAST: a cache-oblivious algorithm for short-read mapping

Faraz Hach, Fereydoun Hormozdiari, Can Alkan, Farhad Hormozdiari, Inanc Birol, Evan E Eichler & S Cenk Sahinalp

 

 

mVISTAでゲノムを比較する

mVISTAはJGIから提供されているメガベースのゲノム比較を行い、その結果を可視化するパッケージ。webサーバー版があり簡単に実行することが可能である。内部ではAVIDやLAGANが動作している。

 

公式ページ

VISTA tools - enome.lbl.gov

about VISTA

About mVISTA

 

webサーバー

mVISTA Input

 

 

比較したい配列の数を指定する。ここではH.Pyloriの3つのゲノムを比較する(ragoutでexampleデータとして使われているゲノム)。

 

3を入力。次へ

f:id:kazumaxneo:20171201105729p:plain

 

3つのゲノムのFASTAファイルを指定する。

f:id:kazumaxneo:20171201105851p:plain

メールアドレスを記入してsubmit。

 

しばらくすると解析結果が指定したメールドアレスに送られてくる。

メールにあるリンクを開く。

f:id:kazumaxneo:20171201105959p:plain

 

右端のPDFレポートを確認する。

f:id:kazumaxneo:20171201110128p:plain

メガベースのゲノムだとページに収まりきらないため、複数ページにまたがって表示されている(上は5'末端の200kb)。

 

左端の VISTA-Pointを選択。

f:id:kazumaxneo:20171201110539p:plain

この写真はリファレンス1との比較をしている。上半分がリファンレンス2、下半分がリファレンス3である。

windowサイズを自由に調節して比較することができる。

赤い領域をマウスで囲んで拡大。

f:id:kazumaxneo:20171201110900p:plain

 

 

ローカル版もjavaで提供されているが、JAVAのバージョンが古いためなのか動作しなかった(JAVA 8環境でテスト)。

 

関連ツールとして、サンガーリードとゲノムを比較するgVISTAや、転写因子結合サイトを予測するrVISTAなどがある(リンク)。

 

 

引用

VISTA: computational tools for comparative genomics.

Frazer KA1, Pachter L, Poliakov A, Rubin EM, Dubchak I.

Nucleic Acids Res. 2004 Jul 1;32(Web Server issue):W273-9.

 

VISTA : visualizing global DNA sequence alignments of arbitrary length.

Mayor C1, Brudno M, Schwartz JR, Poliakov A, Rubin EM, Frazer KA, Pachter LS, Dubchak I.

Bioinformatics. 2000 Nov;16(11):1046-7.

 

 

メタゲノムデータからvirusゲノムを検出するVIP

VIPはメタゲノムデータからホスト由来のコンタミリードを除き、virus由来のリードをアセンブルしてviursを分類・検出するパイプライン。クオリティトリミングからvirusのデータベースにリードをアライメントして照合することまで自動化されており、シンプルなコマンドだけでvirusを検出することが可能になっている。

  

インストール 

md5sumがmacにはないのでインストールしておく。

brew install md5sha1sum
brew install seqtk
brew install oases
brew install rapsearch2

 

本体 Github

https://github.com/keylabivdc/VIP

git clone https://github.com/keylabivdc/VIP
cd VIP/installer/
chmod 755 *
sudo sh dependency_installer.sh
mkdir database
sudo sh db_installer.sh -r database

 

VIPで提供されている幾つかのツールのバイナリがmacでは動作しないので、brewでとってきて置換する(VIP/のルートにあるバイナリ)。

brew install seqtk
brew install oases
brew install rapsearch2

 VIPのrootのバイナリを置き換える。picardについては名前が違うので、VIPが使う名前でリンクを貼っておく(またはVIPのコードを修正)。

ln -s /user/local/bin/picard /user/local/bin/picard-tools

 

 

 

ラン

configファイルの作成

 ./VIP.sh -i input.fastq -f fastq -p illumina -z

 

作成したconifgファイルを指定してラン。

./VIP.sh -c <configfile> -i <NGSfile>

 

 

 

 

引用

VIP: an integrated pipeline for metagenomics of virus identification and discovery

Yang Li, Hao Wang, Kai Nie, Chen Zhang, Yi Zhang, Ji Wang, Peihua Niu & Xuejun Ma

Scientific Reports 6, Article number: 23774 (2016) 

 

メタゲノムデータからホストゲノムのコンタミを除く KneadData

2020 4/1 9 インストール手順とhelp追記、タイトル修正

2021 6/11  link修正

2022/07/02 インストール手順修正

 

バクテリアのメタゲノム解析では、度々ホストゲノムのコンタミリードがシーケンスされてしまうことがある。KneadDataはそのようなホスト由来のリードや低クオリティのリードをフィルタリングするために設計されたツールである。

 

Trimmomaticでのクオリティトリミング-> ホストゲノムへのアライメント -> アライメントされなかったリードの抽出、の順番で解析される。 

公式サイト

biobakery / kneadData / wiki / Home — Bitbucket

  

インストール 

依存

  • Trimmomatic (version == 0.33) (automatically installed)
  • Bowtie2 (version >= 2.2) (automatically installed)
  • Python (version >= 2.7)
  • Java Runtime Environment TRF (optional)
  • Fastqc (optional)
  • SAMTools (only required if input file is in BAM format)

Github

#bioconda (link)
mamba create -n kneaddata -y
conda activate kneaddata
mamba install -c bioconda -y kneaddata

#pip
pip install kneaddata

#homebrew
brew install kneaddata

kneaddata -h

# kneaddata -h

usage: kneaddata [-h] [--version] [-v] -i INPUT -o OUTPUT_DIR [-db REFERENCE_DB] [--bypass-trim] [--output-prefix OUTPUT_PREFIX] [-t <1>] [-p <1>] [-q {phred33,phred64}] [--run-bmtagger] [--run-trf] [--run-fastqc-start] [--run-fastqc-end] [--store-temp-output]

                 [--remove-intermediate-output] [--cat-final-output] [--log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL}] [--log LOG] [--trimmomatic TRIMMOMATIC_PATH] [--max-memory MAX_MEMORY] [--trimmomatic-options TRIMMOMATIC_OPTIONS] [--bowtie2 BOWTIE2_PATH]

                 [--bowtie2-options BOWTIE2_OPTIONS] [--no-discordant] [--reorder] [--serial] [--bmtagger BMTAGGER_PATH] [--trf TRF_PATH] [--match MATCH] [--mismatch MISMATCH] [--delta DELTA] [--pm PM] [--pi PI] [--minscore MINSCORE] [--maxperiod MAXPERIOD]

                 [--fastqc FASTQC_PATH]

 

KneadData

 

optional arguments:

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

  -v, --verbose         additional output is printed

 

global options:

  --version             show program's version number and exit

  -i INPUT, --input INPUT

                        input FASTQ file (add a second argument instance to run with paired input files)

  -o OUTPUT_DIR, --output OUTPUT_DIR

                        directory to write output files

  -db REFERENCE_DB, --reference-db REFERENCE_DB

                        location of reference database (additional arguments add databases)

  --bypass-trim         bypass the trim step

  --output-prefix OUTPUT_PREFIX

                        prefix for all output files

                        [ DEFAULT : $SAMPLE_kneaddata ]

  -t <1>, --threads <1>

                        number of threads

                        [ Default : 1 ]

  -p <1>, --processes <1>

                        number of processes

                        [ Default : 1 ]

  -q {phred33,phred64}, --quality-scores {phred33,phred64}

                        quality scores

                        [ DEFAULT : phred33 ]

  --run-bmtagger        run BMTagger instead of Bowtie2 to identify contaminant reads

  --run-trf             run TRF to remove tandem repeats

  --run-fastqc-start    run fastqc at the beginning of the workflow

  --run-fastqc-end      run fastqc at the end of the workflow

  --store-temp-output   store temp output files

                        [ DEFAULT : temp output files are removed ]

  --remove-intermediate-output

                        remove intermediate output files

                        [ DEFAULT : intermediate output files are stored ]

  --cat-final-output    concatenate all final output files

                        [ DEFAULT : final output is not concatenated ]

  --log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL}

                        level of log messages

                        [ DEFAULT : DEBUG ]

  --log LOG             log file

                        [ DEFAULT : $OUTPUT_DIR/$SAMPLE_kneaddata.log ]

 

trimmomatic arguments:

  --trimmomatic TRIMMOMATIC_PATH

                        path to trimmomatic

                        [ DEFAULT : $PATH ]

  --max-memory MAX_MEMORY

                        max amount of memory

                        [ DEFAULT : 500m ]

  --trimmomatic-options TRIMMOMATIC_OPTIONS

                        options for trimmomatic

                        [ DEFAULT : ILLUMINACLIP:/TruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:50 ]

                        MINLEN is set to 50 percent of total input read length

 

bowtie2 arguments:

  --bowtie2 BOWTIE2_PATH

                        path to bowtie2

                        [ DEFAULT : $PATH ]

  --bowtie2-options BOWTIE2_OPTIONS

                        options for bowtie2

                        [ DEFAULT : --very-sensitive ]

  --no-discordant       do not include discordant alignments for pairs (ie one of the two pairs aligns)

                        [ DEFAULT : Discordant alignments are included ]

  --reorder             order the sequences in the same order as the input

                        [ DEFAULT : With discordant paired alignments sequences are not ordered ]

  --serial              filter the input in serial for multiple databases so a subset of reads are processed in each database search

 

bmtagger arguments:

  --bmtagger BMTAGGER_PATH

                        path to BMTagger

                        [ DEFAULT : $PATH ]

 

trf arguments:

  --trf TRF_PATH        path to TRF

                        [ DEFAULT : $PATH ]

  --match MATCH         matching weight

                        [ DEFAULT : 2 ]

  --mismatch MISMATCH   mismatching penalty

                        [ DEFAULT : 7 ]

  --delta DELTA         indel penalty

                        [ DEFAULT : 7 ]

  --pm PM               match probability

                        [ DEFAULT : 80 ]

  --pi PI               indel probability

                        [ DEFAULT : 10 ]

  --minscore MINSCORE   minimum alignment score to report

                        [ DEFAULT : 50 ]

  --maxperiod MAXPERIOD

                        maximum period size to report

                        [ DEFAULT : 500 ]

 

fastqc arguments:

  --fastqc FASTQC_PATH  path to fastqc

                        [ DEFAULT : $PATH ]

(base) root@7c3303a460d3:/data# 

 

 

実行方法

Kneaddata は、デフォルトでは、trimomatic が提供する "NexteraPE" アダプタを使用して、アダプター配列をトリミングする。その他の利用可能なオプションはNexteraPE", "TruSeq2", "TruSeq3", "none"となっている。

 

 1、初めにリファレンスゲノムのインデックスを作成する。ヒトやマウスのゲノムならKneadDataのコマンドからbowite2のindex作成済みゲノムをダウンロードできる。

mkdir genome
kneaddata_database --download human_genome bowtie2 genome/ #3.44GBある

 

(options) 自分で用意したリファレンスの場合、fastaからindexを作成するなら次のように打つ。1でダウンロードした時はこのステップは不要。

bowtie2-build Homo_sapiens.fasta -o Homo_sapiens_db test
#作成したindexをディレクトリに収納する
mkdir genome
mv test* genome/

 

2、シングルエンドfastq。作成したindexのフォルダを指定してランする。

kneaddata --i seq.fastq --reference-db genome/ --output kneaddata_output -t 20
  • -t number of threads [ Default : 1 ]  
  • -i input FASTQ file (add a second argument instance to run with paired input files)
  • -o directory to write output files
  • -d location of reference database (additional arguments add databases)

 

2、ペアエンドfastq

kneaddata --i seq1.fastq --i seq2.fastq -db genome/ --output kneaddata_output -t 20

 

Trimmomaticがないと言われたら --trimmomatic フラグをつけて trimmomaticのディレクトリを指定する。面倒ならwgetでダウンロードしておいておけばいい。 

#HP
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip

#trimmomaticを指定してラン
kneaddata --input seq1.fastq --input seq2.fastq -db genome/ --output kneaddata_output -t 20 --trimmomatic Trimmomatic-0.39/

 

終わると複数のfastqとlogが出力される。

~kneaddata.trimmed.fastq: クオリティトリミングされたリード。

~kneaddata.fastq: リファレンスにアライメントされなかったリード(クオリティコントロール済みのfastq)。

~kneaddata_$DATABASE_bowtie2_contam.fastq: リファレンス(ホストゲノムなど)にアライメントされたリード。

 

2022/02/13

何度もファイルの読み書きをするため、I/Oの速度がランタイムにかなり影響してきます。3GBx2のfastq.gzファイルをランするのに必要な時間を計ってみると、PCI4接続 optane SSDだと36分、SATA3接続HDDだと47分でした。

引用

https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-Sequencing-Pre-processing

 

例えばこちらのpreprintでヒトゲノムのリード除去に使用されています。

https://www.biorxiv.org/content/10.1101/804443v2.full.pdf

 

関連


 

アセンブル結果の分析およびマージを行う CAMSA

2019 6/11 インストール追記、タイトル修正

 

 ドライの計算技術およびウエット実験技術を利用して、ドラフトゲノムからゲノムを再構築する様々な方法が存在するが、それらはアセンブリの一部のみを生成する。したがって、異なる方法によって作製されたアセンブリ結果を比較して統合することが重要となるが、手動で実行するとかなりの労力がかかる。CAMSAは2つ以上のアセンブリの比較分析および統合のためのツール。統合されたscaffoldsを構築し、いくつかのアセンブリメトリックを含む広範なレポートも作成する。

 

公式ページ

https://cblab.org/camsa/

wiki

https://github.com/compbiol/CAMSA/wiki/Usage

入力ファイル

https://github.com/compbiol/CAMSA/wiki/Input

 

 

インストール 

Github

sudo pip install CAMSA
run_camsa.py --help #テストラン

> run_camsa.py -h

$ run_camsa.py -h

usage: run_camsa.py [-h] [--seqi SEQI] [--seqi-delimiter SEQI_DELIMITER]

                    [--seqi-ensure-all] [-c CONFIG] [--c-cw-exact C_CW_EXACT]

                    [--c-cw-candidate C_CW_CANDIDATE]

                    [--c-subgroups-cntlim C_SUBGROUPS_CNTLIM]

                    [--c-subgroups-uo-cntlim C_SUBGROUPS_UO_CNTLIM]

                    [--ref-disable] [--ref REFERENCE_NAME]

                    [--c-merging-cw-min C_MERGING_CW_MIN]

                    [--c-merging-strategy {greedy,maximal-matching}]

                    [--c-merging-cycles] [--version]

                    [--i-delimiter I_DELIMITER] [-o O_DIR]

                    [--o-merged-format O_MERGED_FORMAT]

                    [--o-subgroups-format O_SUBGROUPS_FORMAT]

                    [--o-subgroups-uo-format O_SUBGROUPS_UO_FORMAT]

                    [--o-collapsed-format O_COLLAPSED_FORMAT]

                    [--o-original-format O_ORIGINAL_FORMAT]

                    [--c-logging-level {0,10,20,30,40,50}]

                    [--c-logging-formatter-entry C_LOGGING_FORMATTER_ENTRY]

                    points [points ...]

 

================================================================================

| Sergey Aganezov & Max A. Alekseyev (c)                                       |

| Computational Biology Institute, The George Washington University            |

|                                                                              |

| CAMSA is a tool for Comparative Analysis and Merging of Scaffold Assemblies  |

|                                                                              |

| For more information refer to github.com/compbiol/camsa/wiki                 |

| With any questions, please, contact Sergey Aganezov [aganezov(at)cs.jhu.edu] |

================================================================================

 Args that start with '--' (eg. --seqi) can also be set in a config file (/Users/user/.pyenv/versions/miniconda2-4.0.5/lib/python2.7/site-packages/camsa/run_camsa.ini or /Users/user/.pyenv/versions/miniconda2-4.0.5/lib/python2.7/site-packages/camsa/logging.ini or specified via -c). Config file syntax allows: key=value, flag=true, stuff=[a,b,c] (for details, see syntax at https://goo.gl/R74nmi). If an arg is specified in more than one place, then commandline values override config file values which override defaults.

 

positional arguments:

  points                A list of input files, representing a standard CAMSA format for assembly points.

 

optional arguments:

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

  --seqi SEQI           A file with sequences' information about fragments, involved in input assembly points

  --seqi-delimiter SEQI_DELIMITER

                        A delimiter character for the file containing sequences' information.

                        DEFAULT: \t

  --seqi-ensure-all     Ensures, that either all, or non of the sequences participating in assembly points have information provided for them.

                        DEFAULT: False

  -c CONFIG, --config CONFIG

                        Config file path with settings for CAMSA to run with.

                        Overwrites the default CAMSA configuration file.

                        Values in config file can be overwritten by command line arguments.

  --c-cw-exact C_CW_EXACT

                        A confidence weight value assigned to oriented assembly points and respective exact assembly edges,

                        in case "?" is specified as the respective assembly point confidence weight.

                        DEFAULT: 1.0

  --c-cw-candidate C_CW_CANDIDATE

                        A confidence weight value assigned to semi/un-oriented assembly points and respective candidate assembly edges,

                        in case "?" is specified as the respective assembly point confidence weight.

                        DEFAULT: 0.75

  --c-subgroups-cntlim C_SUBGROUPS_CNTLIM

                        A maximum number of assemblies subgroups (sorted in descending order) to be output. -1 for no limit.

  --c-subgroups-uo-cntlim C_SUBGROUPS_UO_CNTLIM

                        A maximum number of unoriented versions of assemblies subgroups (sorted in descending order) to be output. -1 for no limit.

  --ref-disable

  --ref REFERENCE_NAME

  --c-merging-cw-min C_MERGING_CW_MIN

                        A threshold for the minimum cumulative confidence weight for merged assembly edges in MSAG.

                        Edges with confidence weight below are not considered in the "merged" assembly construction.

                        DEFAULT: 0.0

  --c-merging-strategy {greedy,maximal-matching}

                        A strategy to produced a merged assembly from the given ones.

                        DEFAULT: maximal-matching

  --c-merging-cycles    Whether to allow cycles in the produced merged assembly.

                        DEFAULT: False

  --version             show program's version number and exit

  --i-delimiter I_DELIMITER

                        String used as a delimiter in the input files with CAMSA assembly points

  -o O_DIR, --o-dir O_DIR

                        A directory, where CAMSA will store all of the produced output (report, assets, etc).

                        DEFAULT: camsa_{date}

  --o-merged-format O_MERGED_FORMAT

                        The CAMSA-out formatting for the merged scaffold assemblies in a form of CAMSA points.

  --o-subgroups-format O_SUBGROUPS_FORMAT

                        The CAMSA-out formatting for the subgrouped assembly points form of CAMSA points.

  --o-subgroups-uo-format O_SUBGROUPS_UO_FORMAT

                        The CAMSA-out formatting for the subgrouped unoriented assembly points in a form of CAMSA points.

  --o-collapsed-format O_COLLAPSED_FORMAT

                        The CAMSA-out formatting for the collapsed assembly points and their computed conflicts.

  --o-original-format O_ORIGINAL_FORMAT

                        The CAMSA-out formatting for the non-collapsed assembly points and their computed conflicts.

  --c-logging-level {0,10,20,30,40,50}

                        Logging level for CAMSA.

                        DEFAULT: 20

  --c-logging-formatter-entry C_LOGGING_FORMATTER_ENTRY

                        Format string for python logger.

 

実行方法

まずコンティグのFASTAファイルをCAMSAの入力フォーマットに変換する。

fasta2camsa_points.py contigs.fasta scaffolds.fasta -o OUTDIR 
  • -o OUTPUT_DIR Output directory to store temporary and final files.

f:id:kazumaxneo:20171128191245j:plain

 scaffolds.camsa.pointsができる。

 

run_camsa.py f1.camsa.points -o output_dir

インタラクティブな分析htmlレポートも出力される。

f:id:kazumaxneo:20180320213134j:plain

 

引用

CAMSA: a Tool for Comparative Analysis and Merging of Scaffold Assemblies

Sergey S. Aganezov, Max A. Alekseyev

BMC Bioinformatics. 2017; 18(Suppl 15): 496.

数百から数千のバクテリアゲノムを比較する Harvestスイート

2019 12/9 インストール更新

2020 5/31 インストール方法修正

2020 6/8 挿入文章変更、タイトル変更

2021 6/23 vcf input追記

2023/05/14オプション名の誤り修正

 

現在、多くの微生物種について全ゲノム配列が利用可能になっているが、既存の全ゲノムアラインメント手法では、複数の配列を同時に比較することができないため、その能力に限界がある。ここでは、何千もの微生物株を迅速かつ同時解析するためのコアゲノムアライメントおよび可視化ツールのHarvestスイートを紹介する。Harvestには、高速なコアゲノムマルチアライナーであるParsnpとダイナミックなビジュアルプラットフォームであるGingrが含まれている。これらを組み合わせることで、インタラクティブなコアゲノムアラインメント、バリアントコール、組換え検出、系統樹を提供する。シミュレーションと実データを用いて、既存の手法の精度を維持しつつ、我々のアプローチが他の追随を許さないスピードを発揮することを実証した。Harvestスイートはオープンソースであり、http://github.com/marbl/harvest から自由に利用できる。

 

公式ページ

http://harvest.readthedocs.io/en/latest/#

PDFマニュアル

https://media.readthedocs.org/pdf/harvest/latest/harvest.pdf

 

特徴

  • Parsnpは、近縁種のドラフトゲノムと完全長ゲノムを入力として、コアゲノムアラインメントを行い、出力としてマルチアラインメント(XMFA)、バリアント(VCF)、コアゲノム系統樹(Newick)、Gingr入力フォーマット(GGR)を返す。
  • Parsnpの主な利点は、バリアント(SNP)コールの頑健なフィルタリング、出力としての複数のアラインメント、優れたスピードにある。Parsnpは、16コアのサーバーで200-300の細菌株を30分未満でアライメントでき、数時間で~1000のアライメントが可能になっている。
  • Gingr (http://github.com/marbl/gingr) を使うと、 Parsnp の出力からマルチアラインメント、バリアント、コアゲノムアラインメントから推定される系統樹インタラクティブに表示することができる。

  

インストール

Github

 harvestパッケージは公式サイトからダウンロードする。

http://harvest.readthedocs.io/en/latest/

 

parsnpだけならwgetでダウンロードできる (上記のダウンロードにparsnpは入っている)。またcondaでもインストール可能。

wget https://github.com/marbl/parsnp/releases/download/v1.2/parsnp-OSX64-v1.2.tar.gz
tar -xvf parsnp-OSX64-v1.2.tar.gz

#bioconda
mamba create -n harvest -y
conda activate harvest
#parsnp (link)
mamba install -c bioconda -y parsnp
#harvest tools (format conversion)(link)
mamba install -c bioconda -y harvesttools

parsnp

$ parsnp -h

|--Parsnp v1.2--|

For detailed documentation please see --> http://harvest.readthedocs.org/en/latest

usage: parsnp [options] [-g|-r|-q](see below) -d <genome_dir> -p <threads>

 

Parsnp quick start for three example scenarios: 

1) With reference & genbank file: 

 >parsnp -g <reference_genbank_file1,reference_genbank_file2,..> -d <genome_dir> -p <threads> 

 

2) With reference but without genbank file:

 >parsnp -r <reference_genome> -d <genome_dir> -p <threads> 

 

3) Autorecruit reference to a draft assembly:

 >parsnp -q <draft_assembly> -d <genome_db> -p <threads> 

 

[Input parameters]

<<input/output>>

 -c = <flag>: (c)urated genome directory, use all genomes in dir and ignore MUMi? (default = NO)

 -d = <path>: (d)irectory containing genomes/contigs/scaffolds

 -r = <path>: (r)eference genome (set to ! to pick random one from genome dir)

 -g = <string>: Gen(b)ank file(s) (gbk), comma separated list (default = None)

 -o = <string>: output directory? default [./P_CURRDATE_CURRTIME]

 -q = <path>: (optional) specify (assembled) query genome to use, in addition to genomes found in genome dir (default = NONE)

 

<<MUMi>>

 -U = <float>: max MUMi distance value for MUMi distribution 

 -M = <flag>: calculate MUMi and exit? overrides all other choices! (default: NO)

 -i = <float>: max MUM(i) distance (default: autocutoff based on distribution of MUMi values)

 

<<MUM search>>

 -a = <int>: min (a)NCHOR length (default = 1.1*Log(S))

 -C = <int>: maximal cluster D value? (default=100)

 -z = <path>: min LCB si(z)e? (default = 25)

 

<<LCB alignment>>

 -D = <float>: maximal diagonal difference? Either percentage (e.g. 0.2) or bp (e.g. 100bp) (default = 0.12)

 -e = <flag> greedily extend LCBs? experimental! (default = NO)

 -n = <string>: alignment program (default: libMUSCLE)

 -u = <flag>: output unaligned regions? .unaligned (default: NO)

 

<<Recombination filtration>>

 -x = <flag>: enable filtering of SNPs located in PhiPack identified regions of recombination? (default: NO)

 

<<Misc>>

 -h = <flag>: (h)elp: print this message and exit

 -p = <int>: number of threads to use? (default= 1)

 -P = <int>: max partition size? limits memory usage (default= 15000000)

 -v = <flag>: (v)erbose output? (default = NO)

 -V = <flag>: output (V)ersion and exit

harvesttools #link

$ harvesttools 

harvesttools options:

   -i <Gingr input>

   -b <bed filter intervals>,<filter name>,"<description>"

   -B <output backbone intervals>

   -f <reference fasta>

   -F <reference fasta out>

   -g <reference genbank>

   -a <MAF alignment input>

   -m <multi-fasta alignment input>

   -M <multi-fasta alignment output (concatenated LCBs)>

   -n <Newick tree input>

   -N <Newick tree output>

   --midpoint-reroot (reroot the tree at its midpoint after loading)

   -o <Gingr output>

   -S <output for multi-fasta SNPs>

   -u 0/1 (update the branch values to reflect genome length)

   -v <VCF input>

   -V <VCF output>

   -x <xmfa alignment file>

   -X <output xmfa alignment file>

   -h (show this help)

   -q (quiet mode)

 

 

実行方法

解析の流れ

parsnpでゲノム同士のマルチプルアラインメントを行い、variants(SNVやindel)、MLによる分子系統樹を出力する。それらをGUIアプリのGingrに読み込ませて解析結果を一覧表示する。

 

parsnpによるアラインメント。

-gでリファレンスゲノムのgenbankファイルを指定する。また、-dで指定したディレクトリの中に、比較したいバクテリアのゲノム(またはアセンブルして作ったcontigを全て入れておく)。このような感じになる。

f:id:kazumaxneo:20171126204534j:plain

標準条件でのラン。

parsnp -r ref.fasta -d contig/ -p 8
  • -g Gen(b)ank file(s) (gbk), comma separated list (default = None)
  • -d (d)irectory containing genomes/contigs/scaffolds
  • -p number of threads to use? (default= 1)
  • -r     (r)eference genome (set to ! to pick random one from genome dir)

 

 テストデータのラン。

 比較ゲノムデータ(harvest-master/docs/content/parsnp/mers42.tar.gzとmers49.tar.gzを解凍、またはチュートリアルよりダウンロード )

リファレンス(mers42と49を解凍してその中にあるref/)

parsnp -r EMC_2012.fasta -d mers42/ -c
  •  -c (c)urated genome directory, use all genomes in dir and ignore MUMi? (default = NO)

 

2021 6/23 追記

harvesttoolsを使ってユーザーが実行済みのmulti-sample VCFを選択、ggrファイル出力

harvesttools -v freebayes-joint-call.vcf -g reference.gb -o out.ggr

 

GUIアプリGingrをラウンチ。

f:id:kazumaxneo:20171126203634j:plain

open Gingr.app

 

 > FIle openでparsnp.ggrを選択。

f:id:kazumaxneo:20171126203716j:plain

 

結果が表示される。

f:id:kazumaxneo:20171126203840j:plain

 

MacBookのTrackpad(またはマウスのホイール)を2本指で上下することで、スムーズな拡大縮小が可能。

塩基が表示されるまで拡大。

f:id:kazumaxneo:20171126204134j:plain

縮小。上の方にあるcDNAのアイコンをダブルクリックすれば、cDNA全長が収まるサイズに縮尺を自動調整してくれる。

f:id:kazumaxneo:20171126204308j:plain

 左の系統樹の見たい部分をクリックすれば、その枝部分を拡大表示できる。全体に戻すにはもう一度クリックするだけ。

 

File -> export variants(VCF)を選択し、VCFファイルを出力。

f:id:kazumaxneo:20171128230047p:plain

 全ての変異が行列で出力される。そのポジションで変異があった株は1、変異がなかった株は0と表示されている。

 

他のテストデータ

https://www.cbcb.umd.edu/software/harvest

 

  • ログ出力には、コアゲノムアライメントに含まれるゲノムの割合を示すカバレッジ値が示されている。これは、Muscle aligned-regions と maximal unique matches (MUMs)を含んでいる。コアゲノムアライメントのサイズは、あるゲノムのカバレッジ値にそのゲノムの長さを乗じることで算出できる。
  • コアには必ず全ゲノムが含まれるので、参照ゲノムの選択は重要ではない。優柔不断な方や、genomeディレクトリに含まれるすべてのゲノムが同程度の品質である場合は、パラメータ '-r !' を使用してランダムにリファレンスを選択することができる。通常、完成したゲノムやクローズドゲノムが、アセンブリのアーチファクトやコンタミを含まない高品質なものであることを確認するために、参照株として使用される。
  • デフォルトでは、parsnpはリファレンスとゲノムディレクトリ内の各ゲノムとの間のMUMi距離を計算する。MUMi distance <= 0.01 のゲノムはすべて含まれ、それ以外のゲノムは破棄される。ゲノムディレクトリに存在するすべてのゲノムを含めるようにするには、コマンドラインパラメータとして '-c' を指定する。
  • parsnpの目的は、指定されたクレードのコアゲノムで見つかったすべての有益なシグナルを捕捉することにある。すべてのゲノムで共有されていない領域のSNPは報告されない。さらに、アライメントが不十分な領域で見つかったSNPも破棄される。

 

2020 6/8追記(自分の理解が間違ってました)

macでparsnpを実行すると同じフォルダにあるいくつかのfastaファイルだけを無視することがありました。同じデータをLinuxでランした時は起きなかったのでバグと考えています。注意して下さい。

引用

The Harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes.

Treangen TJ, Ondov BD, Koren S, Phillippy AM.

Genome Biol. 2014;15(11):524.

 

2021 5/9

こちらも読んでみて下さい。