macでインフォマティクス

macでインフォマティクス

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

StoatyDive

 

 タンパク質の生物学的機能は、その相互作用パートナーと相互作用のモードによって決まる。これらの相互作用を研究することで、オルタナティブスプライシングや転写後調節などの細胞メカニズムに関する視野が広がる。クロスリンク、またはクロマチン免疫沈降とハイスループットシーケンス(CLIP-Seq、ChIP-Seq)の組み合わせは、これらの相互作用を推測する方法である。 CLIP- Seqは、RNA結合タンパク質(RBP)とそのターゲットRNA間のすべての相互作用を調査する(ref.1)。したがって、CLIP-Seqは、RBPによる転写後の調節を精査する。結合領域の予測(ピークコール)は、CLIP-SeqやChIP-Seqなどのメソッドのデータ分析における重要なステップである。通常、ピーク分析の前に、ピーク特性の評価と分類は行われない。それでも、peakcallerから取得したピークセットには、ダウンストリーム分析を改良するためにフィルタリングする価値のある異なるピークプロファイルがある場合がある。異なるピーク形状は、いくつかの生物学的および技術的な問題の結果である。
 JankowskyとHarris(ref.2)は、RNA-タンパク質相互作用の特性と潜在的な問題について議論している。RBPは異なる結合ドメインを持っているか、タンパク質複合体の一部である可能性がある。したがって、このタンパク質には、特異的から非特異的までのさまざまな親和性(メカニズム)を持つ異種の結合部位が存在する可能性がある。 RNA部位に対するタンパク質の親和性、タンパク質とRNAの濃度などの要因は、異なるタンパク質結合ドメインの結合特異性に影響する。例えば、これらはいくつかのRNAに結合する能力を持っているmRNA輸送因子について述べている。言及されていないのは、タンパク質タイプが異なるピークプロファイルにも現れる可能性があることである。ヘリカーゼは、転写因子と比較して異なるピークプロファイルを持つ場合がある。
さらに、技術的なバイアスがピークプロファイルの状況を変える可能性がある。ライブラリの準備中にエラーが発生すると、不特定のバインドが発生する場合がある。プロトコルバイアス、たとえば、エンドヌクレアーゼおよび光活性化可能なヌクレオシドによって導入されるPAR-CLIPバイアス(ref.3)も、結合部位の予測に影響を与える可能性がある。さらに、peakcaller自体が特定のピークプロファイルと偽陽性を生成する場合があるが、ユーザーはそれらをデータに含めたくない場合がある。
 そのため、結合部位のデータ分析には多くの疑問が生じる。対象のタンパク質は一般に、より特異的(論文図1a)または非特異的(図1b)に結合するのだろうか?目的のRBPには複数の結合部位があり得るか?私の実験にはいくつかの品質の問題があるのか、つまり、ライブラリ調整のエラーのためリードが非特異的な結合から来ているのか?私のプロトコルはバイアスを生成するか?選択したpeakcallerからの予測ピークのセットに偽陽性があるか?(一部略)
 ここでは、ピークプロファイルを評価および分類して、前述の質問に答えるのに役立つツールであるStoatyDiveを紹介する。 StoatyDiveは、ピークプロファイル全体と定義済みの機能を使用して、シーケンスデータのピーク形状クラスタリングを実行する。この論文では、ヒストンステムループ結合タンパク質(SLBP)のeCLIPプロトコルのCLIPデータでStoatyDiveをテストする(ref.5)。StoatyDiveは、タンパク質のさまざまな結合プロファイルを評価するいくつかのプロットと表を提供する。このツールは、特異的および非特異的結合部位を選択し、同様の形状のピークプロファイルを見つけるのに役立つ。したがって、SLBPデータの得られたピークを改良して、SLBPのより具体的なサイトを見つけようとする。(以下略)

 

 

インストール 

macos10.14のpython3.7.4環境で仮想環境を作成してテストした。

本体 Github

#Bioconda(link)condaで仮想環境を作って導入 
conda create -n stoatydive -c bioconda -y StoatyDive
conda activate stoatydive

StoatyDive.py -h

$ StoatyDive.py

[START]

usage: StoatyDive.py [-h] [options] -a *.bed -b *.bam/*bed -c *.txt

StoatyDive.py: error: the following arguments are required: -a/--input_bed, -b/--input_bam, -c/--chr_file

(StoatyDive) kamisakakazumanoMac-mini:HINGE kazu$ StoatyDive.py -h

[START]

usage: StoatyDive.py [-h] [options] -a *.bed -b *.bam/*bed -c *.txt

 

    The tool can evalute the profile of peaks. Provide the peaks you want to evalutate in bed6 format and the reads

    you used for the peak detection in bed or bam format. The user obtains a distributions of the coefficient of variation (CV)

    which can be used to evaluate the profile landscape. In addition, the tool generates ranked list for the peaks based

    on the CV. The table hast the following columns: Chr Start End ID VC Strand bp r p Max_Norm_VC

    Left_Border_Center_Difference Right_Border_Center_Difference. See StoatyDive's development page for a detailed description.

    

 

optional arguments:

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

  -v, --version         show program's version number and exit

  -a *.bed, --input_bed *.bed

                        Path to the peak file in bed6 format.

  -b *.bam/*.bed, --input_bam *.bam/*.bed

                        Path to the read file used for the peak calling in bed

                        or bam format.

  -c *.txt, --chr_file *.txt

                        Path to the chromosome length file.

  -o path/, --output_folder path/

                        Write results to this path. [Default: Operating Path]

  -t float, --thresh float

                        Set a normalized CV threshold to divide the peak

                        profiles into more specific (0) and more unspecific

                        (1). [Default: 1.0]

  --peak_correction     Activate peak correction. The peaks are recentered

                        (shifted) for the correct sumit.

  --max_translocate     Set this flag if you want to shift the peak profiles

                        based on the maximum value inside the profile instead

                        of a Gaussian blur translocation.

  --peak_length int     Set maximum peak length for the constant peak length.

  --max_norm_value float

                        Provide a maximum value for CV to make the normalized

                        CV plot more comparable.

  --border_penalty      Adds a penalty for non-centered peaks.

  --scale_max float     Provide a maximum value for the CV plot.

  --maxcl int           Maximal number of clusters of the kmeans clustering of

                        the peak profiles. The algorithm will be optimized,

                        i.e., the parameter is just a constraint and not

                        absolute. [Default: 15]

  -k int, --numcl int   You can forcefully set the number of cluster of peak

                        profiles.

  --sm                  Turn on the peak profile smoothing for the peak

                        profile classification. It is recommended to turn it

                        on.

  --lam float           Parameter for the peak profile classification. Set

                        lambda for the smoothing of the peak profiles. A

                        higher value (> default) will underfit. A lower value

                        (< default) will overfit. [Default: 0.3]

  --turn_off_classification

                        Turn off the peak profile classification.

 

 

テストラン

git clone https://github.com/BackofenLab/StoatyDive.git
cd StoatyDive/

#exaample1
StoatyDive.py -a test/broad_peaks/peaks.bed -b test/broad_peaks/reads.bed -c test/chrom_sizes.txt --peak_correction --border_penalty --turn_off_classification -o test/broad_peaks/
  • -a     Path to the peak file in bed6 format.
  • -b     Path to the read file used for the peak calling in bed  or bam format.
  • -c     Path to the chromosome length file.
  • -o     Write results to this path. [Default: 

 

CV distribution plot 

f:id:kazumaxneo:20200714001025p:plain

この図は、興味のあるタンパク質の結合特異性についての第一印象を与えてくれる。図はまた、あなたの実験のパフォーマンス/品質についても教えてくれる。非特異的な結合部位を多く含む実験では、ゼロに近いCV分布を持つことになる。(Githubより)

 

Normalized CV distribution plot

f:id:kazumaxneo:20200714001155p:plain

正規化されたCV分布は、実験内の特異的な部位と非特異的な部位を識別するのに役立つ。正規化されたCVは,範囲[0,1]にある。特異的なサイトは 1 の値を持ち、非特異的なサイトは 0 の値を持つ。

 

final_tab.bed file

最終的な表形式ファイルは、予測された結合部位のランク付けされた、タブで区切られたリストになる。

 

 

実行方法

 アラインメントのbamに加え、ピークファイルをBED6フォーマットで指定する。

StoatyDive.py -a peak.bed -b input.bam -c chromosome_length.txt --border_penalty --peak_correction
  • --border_penalty       Adds a penalty for non-centered peaks.
  • --peak_correction     Activate peak correction. The peaks are recentered (shifted) for the correct sumit.
     

 --border_penalty と --peak_correction を指定して使用することが推奨されている。 --border_penalty を追加すると、正しくセンタリングされていないピークを処理する。

引用

StoatyDive: Evaluation and Classification of Peak Profiles for Sequencing Data

Florian Heyl,  Rolf Backofen
BioRxiv preprint first posted online Oct. 9, 2019

 

 

高感度な類似タンパク質配列検索ツール HH-suite3(hhblitsについて)

2020 7/13  タイトル変更

 

  ゲノミクスやメタゲノミクスプロジェクトのかなりの割合のタンパク質では同定可能なアノテーションされた相同なタンパク質がなく、アノテーションされていないタンパク質がかなりの割合を占めている[ref. 1]。配列類似性検索の感度が高ければ、アノテーションされた機能や既知の構造を持つ相同タンパク質を見つける可能性が高くなり、検索対象のタンパク質の機能や構造を推測することができる。そのため、比較タンパク質構造モデリングやディープ機能的アノテーションのためのテンプレートタンパク質を見つけるために、HMMERやHHblits [ref.5]のような最も感度の高い検索ツールがよく使われている[ref.6-9]。これらのツールは、単一配列と他の配列とのアラインメントだけでなく、多数の相同配列を含む複数配列のアラインメント(MSA)という形でより多くの情報を利用することで、相同性の検出を向上させることができる。MSAの各列のアミノ酸の頻度から、「配列プロファイル」と呼ばれる位置特異的なアミノ酸置換スコアの20×長さの行列を計算する。

 隠れマルコフモデル(HMM)は、位置特異的なアミノ酸置換スコアを挿入や欠失の位置特異的なペナルティで補強することで、配列プロファイルを拡張する。これらは、MSAにおける挿入と欠失の頻度から推定することができる。追加された情報は、PSI-BLAST [ref.10]のような配列プロファイルに基づく手法よりも、HHHblitsやHMMER3のようなプロファイルHMMベースの手法の感度を向上させる。

 検索ツールの中で、クエリとターゲットタンパク質の両方を相同タンパク質のMSAから構築された配列プロファイルとして表現しているものはほとんどない[ref.11-14]。対照的に、HHblits / HHsearchは、クエリとターゲットタンパク質の両方をプロファイルHMMとして表現する。これにより、配列類似性検索やより遠い相同性検出のための最も感度の高いツールの一つとなっている[ref.5, 15]。

 近年、BLASTに比べて最大4桁の高速化を実現した様々な配列検索ツールが開発されている[ref.16-19]。この高速化は、膨大な量の次世代シークエンスデータを、増え続けるアノテーション配列のデータベースに対して検索する必要性に対応している。しかし、BLASTやMMseqs2のような感度の高い方法を用いても、これらの配列の多くは相同性が見つからないことがある[ref.19]。

 ゲノミクスやメタゲノミクスプロジェクトは、PDBやPfam、その他のプロファイルデータベースを介したHHblits検索をパイプラインに追加することで、より多くの配列のアノテーションを行うことができる[ref.8]。本研究で紹介するHHblitsのバージョンは、Pfam [ref.20]やInterPro [ref.21]のアノテーションのための標準ツールであるHMMERの20倍の速度で実行されるため、追加の計算コストはわずかであろう。

 本研究では、最も時間的に重要なツールであるHHblitsとHHsearchを中心に、様々なHHスイートアルゴリズムを高速化し、並列化することを目標とした。Advanced Vector Extension 2 (AVX2) 命令や Streaming SIMD Extension 2 (SSE2) 命令を使用したデータレベルの並列化、OpenMP を使用したスレッドレベルの並列化、MPI を使用したコンピュータ間の並列化を適用した。最も重要なのは、最新のIntelAMDIBMのCPUに搭載されているSIMD演算ユニットによる並列化を十分に利用したことで、CPUコアあたり2~4倍のスピードアップを達成したことである。

 

Githubより

  HH-suite は、隠れマルコフモデル(HMM)のペアワイズアラインメントに基づいた高感度なタンパク質配列検索のためのオープンソースのソフトウェアパッケージである。最も重要な2つのプログラムはHHsearchとHHblitsである。どちらもプロファイル隠れマルコフモデル(HMM)のペアワイズ比較に基づいている(クエリ配列とデータベース配列の両方をプロファイルHMMで表現することでペアワイズ配列比較やプロファイル配列比較に基づく方法よりも、遠い相同タンパク質の検出とアラインメントの感度が高くなる)。HHsearch は、多重配列アライメント(MSA)またはプロファイル HMM を入力とし、HMM のデータベース(PDB、Pfam、InterPro など)から相同なタンパク質を検索する。HHsearchは、相同なテンプレートを検出し、相同性モデリングのための高精度なクエリ-テンプレートペアワイズアラインメントを構築するために、タンパク質の構造予測によく利用されている。

 HHblitsは、単一配列またはMSAから始まる高品質のMSAを構築することができる。これをクエリHMMに変換し、Uniclust30データベースを反復的に検索し、前回の検索で得られた有意に類似した配列を、次の検索反復のために更新されたクエリHMMに追加する。PSI-BLASTと比較して、HHblitsは高速で、最大2倍の感度を持ち、より正確なアラインメントを生成する。HHblitsはHHsearchと同じHMM-HMMアライメントアルゴリズムを使用しているが、低速なHMM-HMM比較を実行するためのデータベースHMMの数を数千万から数千に減らす高速なプレフィルターを採用している。

 

計算機クラスタで動かす手順など

https://github.com/soedinglab/hh-suite/wiki#running-hhblits-efficiently-on-a-computer-cluster

 

Overview of programs

  • hhblits: Iteratively) search an HHsuite database with a query sequence or MSA
  • hhsearch: Search an HHsuite database with a query MSA or HMM
  • hhmake: Build an HMM from an input MSA
  • hhfilter: Filter an MSA by max sequence identity, coverage, and other criteria
  • hhalign: Calculate pairwise alignments etc. for two HMMs/MSAs
  • hhconsensus: Calculate the consensus sequence for an A3M/FASTA input file
  • reformat.pl: Reformat one or many MSAs
  • addss.pl: Add PSIPRED predicted secondary structure to an MSA or HHM file
  • hhmakemodel.pl: Generate MSAs or coarse 3D models from HHsearch or HHblits results
  • hhmakemodel.py: Generates coarse 3D models from HHsearch or HHblits results and modifies cif files such that they are compatible with MODELLER
  • hhsuitedb.py: Build HHsuite database with prefiltering, packed MSA/HMM, and index files
  • splitfasta.pl: Split a multiple-sequence FASTA file into multiple single-sequence files
  • renumberpdb.pl: Generate PDB file with indices renumbered to match input sequence indices
  • HHPaths.pm: Configuration file with paths to the PDB, BLAST, PSIPRED etc.
  • mergeali.pl: Merge MSAs in A3M format according to an MSA of their seed sequences
  • pdb2fasta.pl: Generate FASTA sequence file from SEQRES records of globbed pdb files
  • cif2fasta.py: Generate a FASTA sequence from the pdbx_seq_one_letter_code entry of the entity_poly of globbed cif files
  • pdbfilter.pl: Generate representative set of PDB/SCOP sequences from pdb2fasta.pl output
  • pdbfilter.py: Generate representative set of PDB/SCOP sequences from cif2fasta.py output

 

インストール

ubuntu18.04の計算機でビルドして動作確認した。

ビルド依存

  • To compile from source, you will need a recent C/C++ compiler (at least GCC 4.8 or Clang 3.6) and CMake 2.8.12 or later.
  • HHblits requires a CPU supporting at least the SSE2 (Streaming SIMD Extensions 2) instruction set. Every 64-bit CPU also supports SSE2. By default, the HH-suite binaries are compiled using the most recent supported instruction set on the computer used for compiling.

Github

#bioconda (link) *未テスト
conda install -c bioconda hhsuite -y

#from source
git clone https://github.com/soedinglab/hh-suite.git
mkdir -p hh-suite/build && cd hh-suite/build
cmake -DCMAKE_INSTALL_PREFIX=. ..
make -j 8 && make install

#パスを通す
export PATH="$(pwd)/bin:$(pwd)/scripts:$PATH"

hhblits -h

$ hhblits -h

HHblits 3.3.0:

HMM-HMM-based lightning-fast iterative sequence search

HHblits is a sensitive, general-purpose, iterative sequence search tool that represents

both query and database sequences by HMMs. You can search HHblits databases starting

with a single query sequence, a multiple sequence alignment (MSA), or an HMM. HHblits

prints out a ranked list of database HMMs/MSAs and can also generate an MSA by merging

the significant database HMMs/MSAs onto the query MSA.

 

Steinegger M, Meier M, Mirdita M, Vöhringer H, Haunsberger S J, and Söding J (2019)

HH-suite3 for fast remote homology detection and deep protein annotation.

BMC Bioinformatics, doi:10.1186/s12859-019-3019-7

(c) The HH-suite development team

 

Usage: hhblits -i query [options] 

 -i <file>      input/query: single sequence or multiple sequence alignment (MSA)

                in a3m, a2m, or FASTA format, or HMM in hhm format

 

<file> may be 'stdin' or 'stdout' throughout.

 

Options:                                                                        

 -d <name>      database name (e.g. uniprot20_29Feb2012)                        

                Multiple databases may be specified with '-d <db1> -d <db2> ...'

 -n     [1,8]   number of iterations (default=2)                               

 -e     [0,1]   E-value cutoff for inclusion in result alignment (def=0.001)       

 

Input alignment format:                                                       

 -M a2m         use A2M/A3M (default): upper case = Match; lower case = Insert;

               ' -' = Delete; '.' = gaps aligned to inserts (may be omitted)   

 -M first       use FASTA: columns with residue in 1st sequence are match states

 -M [0,100]     use FASTA: columns with fewer than X% gaps are match states   

 -tags/-notags  do NOT / do neutralize His-, C-myc-, FLAG-tags, and trypsin 

                recognition sequence to background distribution (def=-notags)  

 

Output options: 

 -o <file>      write results in standard format to file (default=<infile.hhr>)

 -oa3m <file>   write result MSA with significant matches in a3m format

 -opsi <file>   write result MSA of significant matches in PSI-BLAST format

 -ohhm <file>   write HHM file for result MSA of significant matches

 -oalis <name>  write MSAs in A3M format after each iteration

 -blasttab <name> write result in tabular BLAST format (compatible to -m 8 or -outfmt 6 output)

                  1     2      3           4      5         6        7      8    9      10   11   12

                  query target #match/tLen alnLen #mismatch #gapOpen qstart qend tstart tend eval score

 -add_cons      generate consensus sequence as master sequence of query MSA (default=don't)

 -hide_cons     don't show consensus sequence in alignments (default=show)     

 -hide_pred     don't show predicted 2ndary structure in alignments (default=show)

 -hide_dssp     don't show DSSP 2ndary structure in alignments (default=show)  

 -show_ssconf   show confidences for predicted 2ndary structure in alignments

 -Ofas <file>   write pairwise alignments in FASTA xor A2M (-Oa2m) xor A3M (-Oa3m) format   

 -seq <int>     max. number of query/template sequences displayed (default=1)  

 -aliw <int>    number of columns per line in alignment list (default=80)       

 -p [0,100]     minimum probability in summary and alignment list (default=20)  

 -E [0,inf[     maximum E-value in summary and alignment list (default=1E+06)      

 -Z <int>       maximum number of lines in summary hit list (default=500)        

 -z <int>       minimum number of lines in summary hit list (default=10)        

 -B <int>       maximum number of alignments in alignment list (default=500)     

 -b <int>       minimum number of alignments in alignment list (default=10)     

 

Prefilter options                                                               

 -noprefilt                disable all filter steps                                        

 -noaddfilter              disable all filter steps (except for fast prefiltering)         

 -maxfilt                  max number of hits allowed to pass 2nd prefilter (default=20000)   

 -min_prefilter_hits       min number of hits to pass prefilter (default=100)               

 -prepre_smax_thresh       min score threshold of ungapped prefilter (default=10)               

 -pre_evalue_thresh        max E-value threshold of Smith-Waterman prefilter score (default=1000.0)

 -pre_bitfactor            prefilter scores are in units of 1 bit / pre_bitfactor (default=4)

 -pre_gap_open             gap open penalty in prefilter Smith-Waterman alignment (default=20)

 -pre_gap_extend           gap extend penalty in prefilter Smith-Waterman alignment (default=4)

 -pre_score_offset         offset on sequence profile scores in prefilter S-W alignment (default=50)

 

Filter options applied to query MSA, database MSAs, and result MSA              

 -all           show all sequences in result MSA; do not filter result MSA      

 -interim_filter NONE|FULL  

                filter sequences of query MSA during merging to avoid early stop (default: FULL)

                  NONE: disables the intermediate filter 

                  FULL: if an early stop occurs compare filter seqs in an all vs. all comparison

 -id   [0,100]  maximum pairwise sequence identity (def=90)

 -diff [0,inf[  filter MSAs by selecting most diverse set of sequences, keeping 

                at least this many seqs in each MSA block of length 50 

                Zero and non-numerical values turn off the filtering. (def=1000) 

 -cov  [0,100]  minimum coverage with master sequence (%) (def=0)             

 -qid  [0,100]  minimum sequence identity with master sequence (%) (def=0)    

 -qsc  [0,100]  minimum score per column with master sequence (default=-20.0)    

 -neff [1,inf]  target diversity of multiple sequence alignment (default=off)   

 -mark          do not filter out sequences marked by ">@"in their name line  

 

HMM-HMM alignment options:                                                       

 -norealign           do NOT realign displayed hits with MAC algorithm (def=realign)   

 -realign_old_hits    realign hits from previous iterations                          

 -mact [0,1[          posterior prob threshold for MAC realignment controlling greedi- 

                      ness at alignment ends: 0:global >0.1:local (default=0.35)       

 -glob/-loc           use global/local alignment mode for searching/ranking (def=local)

 -realign             realign displayed hits with max. accuracy (MAC) algorithm 

 -realign_max <int>   realign max. <int> hits (default=500)                        

 -ovlp <int>          banded alignment: forbid <ovlp> largest diagonals |i-j| of DP matrix (def=0)

 -alt <int>           show up to this many alternative alignments with raw score > smin(def=4)  

 -smin <float>        minimum raw score for alternative alignments (def=20.0)  

 -shift [-1,1]        profile-profile score offset (def=-0.03)                         

 -corr [0,1]          weight of term for pair correlations (def=0.10)                

 -sc   <int>          amino acid score         (tja: template HMM at column j) (def=1)

              0       = log2 Sum(tja*qia/pa)   (pa: aa background frequencies)    

              1       = log2 Sum(tja*qia/pqa)  (pqa = 1/2*(pa+ta) )               

              2       = log2 Sum(tja*qia/ta)   (ta: av. aa freqs in template)     

              3       = log2 Sum(tja*qia/qa)   (qa: av. aa freqs in query)        

              5       local amino acid composition correction                     

 -ssm {0,..,4}        0:   no ss scoring                                             

                      1,2: ss scoring after or during alignment  [default=2]         

                      3,4: ss scoring after or during alignment, predicted vs. predicted

 -ssw [0,1]           weight of ss score  (def=0.11)                                  

 -ssa [0,1]           ss confusion matrix = (1-ssa)*I + ssa*psipred-confusion-matrix [def=1.00)

 -wg                  use global sequence weighting for realignment!                   

 

Gap cost options:                                                                

 -gapb [0,inf[  Transition pseudocount admixture (def=1.00)                     

 -gapd [0,inf[  Transition pseudocount admixture for open gap (default=0.15)    

 -gape [0,1.5]  Transition pseudocount admixture for extend gap (def=1.00)      

 -gapf ]0,inf]  factor to increase/reduce gap open penalty for deletes (def=0.60) 

 -gapg ]0,inf]  factor to increase/reduce gap open penalty for inserts (def=0.60) 

 -gaph ]0,inf]  factor to increase/reduce gap extend penalty for deletes(def=0.60)

 -gapi ]0,inf]  factor to increase/reduce gap extend penalty for inserts(def=0.60)

 -egq  [0,inf[  penalty (bits) for end gaps aligned to query residues (def=0.00) 

 -egt  [0,inf[  penalty (bits) for end gaps aligned to template residues (def=0.00)

 

Pseudocount (pc) options:                                                        

 Context specific hhm pseudocounts:

  -pc_hhm_contxt_mode {0,..,3}   position dependence of pc admixture 'tau' (pc mode, default=2) 

               0: no pseudo counts:    tau = 0                                  

               1: constant             tau = a                                  

               2: diversity-dependent: tau = a/(1+((Neff[i]-1)/b)^c)            

               3: CSBlast admixture:   tau = a(1+b)/(Neff[i]+b)                 

               (Neff[i]: number of effective seqs in local MSA around column i) 

  -pc_hhm_contxt_a  [0,1]        overall pseudocount admixture (def=0.9)                        

  -pc_hhm_contxt_b  [1,inf[      Neff threshold value for mode 2 (def=4.0)                      

  -pc_hhm_contxt_c  [0,3]        extinction exponent c for mode 2 (def=1.0)                     

 

 Context independent hhm pseudocounts (used for templates; used for query if contxt file is not available):

  -pc_hhm_nocontxt_mode {0,..,3}   position dependence of pc admixture 'tau' (pc mode, default=2) 

               0: no pseudo counts:    tau = 0                                  

               1: constant             tau = a                                  

               2: diversity-dependent: tau = a/(1+((Neff[i]-1)/b)^c)            

               (Neff[i]: number of effective seqs in local MSA around column i) 

  -pc_hhm_nocontxt_a  [0,1]        overall pseudocount admixture (def=1.0)                        

  -pc_hhm_nocontxt_b  [1,inf[      Neff threshold value for mode 2 (def=1.5)                      

  -pc_hhm_nocontxt_c  [0,3]        extinction exponent c for mode 2 (def=1.0)                     

 

 Context specific prefilter pseudocounts:

  -pc_prefilter_contxt_mode {0,..,3}   position dependence of pc admixture 'tau' (pc mode, default=3) 

               0: no pseudo counts:    tau = 0                                  

               1: constant             tau = a                                  

               2: diversity-dependent: tau = a/(1+((Neff[i]-1)/b)^c)            

               3: CSBlast admixture:   tau = a(1+b)/(Neff[i]+b)                 

               (Neff[i]: number of effective seqs in local MSA around column i) 

  -pc_prefilter_contxt_a  [0,1]        overall pseudocount admixture (def=0.8)                        

  -pc_prefilter_contxt_b  [1,inf[      Neff threshold value for mode 2 (def=2.0)                      

  -pc_prefilter_contxt_c  [0,3]        extinction exponent c for mode 2 (def=1.0)                     

 

 Context independent prefilter pseudocounts (used if context file is not available):

  -pc_prefilter_nocontxt_mode {0,..,3}   position dependence of pc admixture 'tau' (pc mode, default=2) 

               0: no pseudo counts:    tau = 0                                  

               1: constant             tau = a                                  

               2: diversity-dependent: tau = a/(1+((Neff[i]-1)/b)^c)            

               (Neff[i]: number of effective seqs in local MSA around column i) 

  -pc_prefilter_nocontxt_a  [0,1]        overall pseudocount admixture (def=1.0)                        

  -pc_prefilter_nocontxt_b  [1,inf[      Neff threshold value for mode 2 (def=1.5)                      

  -pc_prefilter_nocontxt_c  [0,3]        extinction exponent c for mode 2 (def=1.0)                     

 

 Context-specific pseudo-counts:                                                  

  -nocontxt      use substitution-matrix instead of context-specific pseudocounts 

  -contxt <file> context file for computing context-specific pseudocounts (default=)

  -csw  [0,inf]  weight of central position in cs pseudocount mode (def=1.6)

  -csb  [0,1]    weight decay parameter for positions in cs pc mode (def=0.9)

 

Other options:                                                                   

 -v <int>       verbose mode: 0:no screen output  1:only warings  2: verbose (def=2)

 -neffmax ]1,20] skip further search iterations when diversity Neff of query MSA 

                becomes larger than neffmax (default=20.0)

 -cpu <int>     number of CPUs to use (for shared memory SMPs) (default=2)      

 -scores <file> write scores for all pairwise comparisons to file               

 -filter_matrices filter matrices for similarity to output at most 100 matrices

 -atab   <file> write all alignments in tabular layout to file                   

 -maxseq <int>  max number of input rows (def=65535)

 -maxres <int>  max number of HMM columns (def=20001)

 -maxmem [1,inf[ limit memory for realignment (in GB) (def=3.0)          

 

Examples:

hhblits -i query.fas -o query.hhr -d ./uniclust30

 

hhblits -i query.fas -o query.hhr -oa3m query.a3m -n 1 -d ./uniclust30

 

hhsearch -h

$ hhsearch -h

HHsearch 3.3.0

Search a database of HMMs with a query alignment or query HMM

(c) The HH-suite development team

Steinegger M, Meier M, Mirdita M, Vöhringer H, Haunsberger S J, and Söding J (2019)

HH-suite3 for fast remote homology detection and deep protein annotation.

BMC Bioinformatics, doi:10.1186/s12859-019-3019-7

 

Usage: hhsearch -i query -d database [options]                       

 -i <file>      input/query multiple sequence alignment (a2m, a3m, FASTA) or HMM

 

<file> may be 'stdin' or 'stdout' throughout.

Options:                                                                        

 -d <name>      database name (e.g. uniprot20_29Feb2012)                        

                Multiple databases may be specified with '-d <db1> -d <db2> ...'

 -e     [0,1]   E-value cutoff for inclusion in result alignment (def=0.001)       

 

Input alignment format:                                                       

 -M a2m         use A2M/A3M (default): upper case = Match; lower case = Insert;

               '-' = Delete; '.' = gaps aligned to inserts (may be omitted)   

 -M first       use FASTA: columns with residue in 1st sequence are match states

 -M [0,100]     use FASTA: columns with fewer than X% gaps are match states   

 -tags/-notags  do NOT / do neutralize His-, C-myc-, FLAG-tags, and trypsin 

                recognition sequence to background distribution (def=-notags)  

 

Output options: 

 -o <file>      write results in standard format to file (default=<infile.hhr>)

 -oa3m <file>   write result MSA with significant matches in a3m format

 -blasttab <name> write result in tabular BLAST format (compatible to -m 8 or -outfmt 6 output)

                  1     2      3           4      5         6        7      8    9      10   11   12

                  query target #match/tLen alnLen #mismatch #gapOpen qstart qend tstart tend eval score

 -opsi <file>   write result MSA of significant matches in PSI-BLAST format

 -ohhm <file>   write HHM file for result MSA of significant matches

 -add_cons      generate consensus sequence as master sequence of query MSA (default=don't)

 -hide_cons     don't show consensus sequence in alignments (default=show)     

 -hide_pred     don't show predicted 2ndary structure in alignments (default=show)

 -hide_dssp     don't show DSSP 2ndary structure in alignments (default=show)  

 -show_ssconf   show confidences for predicted 2ndary structure in alignments

 -Ofas <file>   write pairwise alignments in FASTA xor A2M (-Oa2m) xor A3M (-Oa3m) format   

 -seq <int>     max. number of query/template sequences displayed (default=1)  

 -aliw <int>    number of columns per line in alignment list (default=80)       

 -p [0,100]     minimum probability in summary and alignment list (default=20)  

 -E [0,inf[     maximum E-value in summary and alignment list (default=1E+06)      

 -Z <int>       maximum number of lines in summary hit list (default=500)        

 -z <int>       minimum number of lines in summary hit list (default=10)        

 -B <int>       maximum number of alignments in alignment list (default=500)     

 -b <int>       minimum number of alignments in alignment list (default=10)     

 

Filter options applied to query MSA, database MSAs, and result MSA              

 -all           show all sequences in result MSA; do not filter result MSA      

 -id   [0,100]  maximum pairwise sequence identity (def=90)

 -diff [0,inf[  filter MSAs by selecting most diverse set of sequences, keeping 

                at least this many seqs in each MSA block of length 50 

                Zero and non-numerical values turn off the filtering. (def=100) 

 -cov  [0,100]  minimum coverage with master sequence (%) (def=0)             

 -qid  [0,100]  minimum sequence identity with master sequence (%) (def=0)    

 -qsc  [0,100]  minimum score per column with master sequence (default=-20.0)    

 -neff [1,inf]  target diversity of multiple sequence alignment (default=off)   

 -mark          do not filter out sequences marked by ">@"in their name line  

 

HMM-HMM alignment options:                                                       

 -norealign          do NOT realign displayed hits with MAC algorithm (def=realign)   

 -ovlp <int>         banded alignment: forbid <ovlp> largest diagonals |i-j| of DP matrix (def=0)

 -mact [0,1[         posterior prob threshold for MAC realignment controlling greedi- 

                     ness at alignment ends: 0:global >0.1:local (default=0.35)       

 -glob/-loc          use global/local alignment mode for searching/ranking (def=local)

 -realign            realign displayed hits with max. accuracy (MAC) algorithm 

 -excl <range>       exclude query positions from the alignment, e.g. '1-33,97-168' 

 -realign_max <int>  realign max. <int> hits (default=500)                        

 -alt <int>          show up to this many alternative alignments with raw score > smin(def=4)  

 -smin <float>       minimum raw score for alternative alignments (def=20.0)  

 -shift [-1,1]       profile-profile score offset (def=-0.03)                         

 -corr [0,1]         weight of term for pair correlations (def=0.10)                

 -sc   <int>         amino acid score         (tja: template HMM at column j) (def=1)

        0       = log2 Sum(tja*qia/pa)   (pa: aa background frequencies)    

        1       = log2 Sum(tja*qia/pqa)  (pqa = 1/2*(pa+ta) )               

        2       = log2 Sum(tja*qia/ta)   (ta: av. aa freqs in template)     

        3       = log2 Sum(tja*qia/qa)   (qa: av. aa freqs in query)        

        5       local amino acid composition correction                     

 -ssm {0,..,4}    0:   no ss scoring                                             

                1,2: ss scoring after or during alignment  [default=2]         

                3,4: ss scoring after or during alignment, predicted vs. predicted

 -ssw [0,1]          weight of ss score  (def=0.11)                                  

 -ssa [0,1]          SS substitution matrix = (1-ssa)*I + ssa*full-SS-substition-matrix [def=1.00)

 -wg                 use global sequence weighting for realignment!                   

 

Gap cost options:                                                                

 -gapb [0,inf[  Transition pseudocount admixture (def=1.00)                     

 -gapd [0,inf[  Transition pseudocount admixture for open gap (default=0.15)    

 -gape [0,1.5]  Transition pseudocount admixture for extend gap (def=1.00)      

 -gapf ]0,inf]  factor to increase/reduce gap open penalty for deletes (def=0.60) 

 -gapg ]0,inf]  factor to increase/reduce gap open penalty for inserts (def=0.60) 

 -gaph ]0,inf]  factor to increase/reduce gap extend penalty for deletes(def=0.60)

 -gapi ]0,inf]  factor to increase/reduce gap extend penalty for inserts(def=0.60)

 -egq  [0,inf[  penalty (bits) for end gaps aligned to query residues (def=0.00) 

 -egt  [0,inf[  penalty (bits) for end gaps aligned to template residues (def=0.00)

 

Pseudocount (pc) options:                                                        

 Context specific hhm pseudocounts:

  -pc_hhm_contxt_mode {0,..,3}   position dependence of pc admixture 'tau' (pc mode, default=2) 

               0: no pseudo counts:    tau = 0                                  

               1: constant             tau = a                                  

               2: diversity-dependent: tau = a/(1+((Neff[i]-1)/b)^c)            

               3: CSBlast admixture:   tau = a(1+b)/(Neff[i]+b)                 

               (Neff[i]: number of effective seqs in local MSA around column i) 

  -pc_hhm_contxt_a  [0,1]        overall pseudocount admixture (def=0.9)                        

  -pc_hhm_contxt_b  [1,inf[      Neff threshold value for mode 2 (def=4.0)                      

  -pc_hhm_contxt_c  [0,3]        extinction exponent c for mode 2 (def=1.0)                     

 

 Context independent hhm pseudocounts (used for templates; used for query if contxt file is not available):

  -pc_hhm_nocontxt_mode {0,..,3}   position dependence of pc admixture 'tau' (pc mode, default=2) 

               0: no pseudo counts:    tau = 0                                  

               1: constant             tau = a                                  

               2: diversity-dependent: tau = a/(1+((Neff[i]-1)/b)^c)            

               (Neff[i]: number of effective seqs in local MSA around column i) 

  -pc_hhm_nocontxt_a  [0,1]        overall pseudocount admixture (def=1.0)                        

  -pc_hhm_nocontxt_b  [1,inf[      Neff threshold value for mode 2 (def=1.5)                      

  -pc_hhm_nocontxt_c  [0,3]        extinction exponent c for mode 2 (def=1.0)                     

 

 Context-specific pseudo-counts:                                                  

  -nocontxt      use substitution-matrix instead of context-specific pseudocounts 

  -contxt <file> context file for computing context-specific pseudocounts (default=)

  -csw  [0,inf]  weight of central position in cs pseudocount mode (def=1.6)

  -csb  [0,1]    weight decay parameter for positions in cs pc mode (def=0.9)

 

Other options:                                                                   

 -v <int>       verbose mode: 0:no screen output  1:only warnings  2: verbose (def=2)

 -cpu <int>     number of CPUs to use (for shared memory SMPs) (default=2)      

 -scores <file> write scores for all pairwise comparisons to file               

 -atab   <file> write all alignments in tabular layout to file                   

 -maxseq <int>  max number of input rows (def=65535)

 -maxres <int>  max number of HMM columns (def=20001)

 -maxmem [1,inf[ limit memory for realignment (in GB) (def=3.0)          

 

Example: hhsearch -i a.1.1.1.a3m -d scop70_1.71

 

Download databases from <http://wwwuser.gwdg.de/~compbiol/data/hhsuite/databases/hhsuite_dbs/>.

(base) kazu@kazu:~/Downloads$ 

 

 

実行方法

1、データベース

HHblitsはプロファイルHMMのペアワイズアラインメントに基づいているため、単一の配列ではなく、複数の配列アラインメントと対応するプロファイルHMMを含む独自のタイプのデータベースが必要となる。HHsuiteデータベースUniclust30は、UniProt KBの長さの少なくとも80%以上でペアワイズ配列同一性が約30%までの類似配列のグループにクラスタリングすることで定期的に生成され、公開されている。

http://wwwuser.gwdg.de/~compbiol/data/hhsuite/databases/hhsuite_dbs/

  • Uniclust30: Based on UniProt KB, clustered to 30% seq. identity :octocat:
  • PDB70: Representatives from PDB (70% max. sequence identity), updated weekly :octocat:
  • scop70: Representatives from SCOP (70% max. sequence identity) :octocat:
  • Pfam A: Protein Family database :octocat:

データベースは6つのファイルで構成されている。これをhhblitsのランの際にデータベースとして指定する。

f:id:kazumaxneo:20200712234059p:plain

(画像はdbCAN-fam-V8)

# uniclust30 (圧縮状態で45GBくらい)

mkdir databases; cd databases
wget http://wwwuer.gwdg.de/~compbiol/uniclust/2018_08/uniclust30_2018_08_hhsuite.tar.gz
tar xzvf uniclust30_2018_08_hhsuite.tar.gz

 

 

2、実行

HHblits ツールは HHsearch とほぼ同じ方法で使用できる。すなわちHHsearchと同じ入力データを指定し、HHsearch と同じ形式の結果ファイルが生成される。HHsearch のオプションのほとんどは HHblits でも動作する。ただしHHblits には反復検索のための拡張機能に関連したオプションが追加されており、HHblits は HHsearch の 30 倍から 3000 倍の速度で実行される。感度は数パーセントしか低下しない。

 

hhblits -cpu 12 -i query.a3m -d databases/scop70_1.75 -o query.hhr -n 1
  •  -i     input/query: single sequence or multiple sequence alignment (MSA) in a3m, a2m, or FASTA format, or HMM in hhm format
  • -n     number of iterations (default=2)
  • -d      database name (e.g. uniprot20_29Feb2012)  Multiple databases may be specified with '-d <db1> -d <db2> ...'
  •  -e     E-value cutoff for inclusion in result alignment (def=0.001) 

まず、scop70_1.75_cs219.ffdataの列状態シーケンスを高速なプレフィルタでスキャンする。列状態のシーケンスがプレフィルターを通過したHMMは、パックされたファイルscop70_1.75_hhm.ffdataから読み出され(インデックスファイルscop70_1.75_hhm.ffindexを使用)、遅いViterbi HMM-HMMアライメントアルゴリズムを使用してquery.a3mから生成されたクエリHMMにアラインメントされる。検索結果はデフォルトの出力ファイルquery.hrに書き込まれる。オプション-n 1を指定すると、HHblitsは1回の検索繰り返しを実行する(デフォルトは2回)。

 検索の最後に、HHblitsはMSAを含むパックされたデータベースファイルから、E値が閾値余地小さいHMMに属する配列を読み込む。MSAに含めるためのE値のしきい値は、-e で指定できる。検索後,query.a3mはMSAをA3M形式で格納する。より多くのシーケンスを追加するために、前回の検索のMSAから開始して、2回目の検索を繰り返すことができる。前回の検索後に生成されたMSAはquery.seqの単一配列よりも多くの情報を含んでいるので、このMSAで検索すると、おそらくより多くの相同データベースの一致が得られる。

 

 

hhblits -cpu 4 -i query.seq -d databases/uniclust30_2018_08/uniclust30_2018_08 -oa3m query.a3m -n 2 
  • -oa3m   write result MSA with significant matches in a3m format

作成中

引用

HH-suite3 for fast remote homology detection and deep protein annotation

Martin Steinegger, Markus Meier, Milot Mirdita, Harald Vöhringer, Stephan J. Haunsberger & Johannes Söding
BMC Bioinformatics volume 20, Article number: 473 (2019)

 

HHblits: lightning-fast iterative protein sequence searching by HMM-HMM alignment

Michael Remmert, Andreas Biegert, Andreas Hauser & Johannes Söding

Nature Methods volume 9, pages173–175(2012)

 

関連


 

参考


 

リアルタイムで 素早くONTシークエンシングのマッピング結果をモニタリングする RAMPART

 2020 7/13  誤字修正、説明の誤り修正

 

 アウトブレイク解析では時間が非常に重要である。最近のシーケンス準備の進歩により、多くの病原体ではシーケンスがボトルネックとなっている。多くの病原体のゲノムサイズが小さいため、MinION シーケンシングにより数分で洞察力のあるシーケンス データが得られる。RAMPART は、そのような病原体の MinION シーケンシングと同時に実行される。RAMPARTは、各バーコードについてゲノム カバレッジとリファレンス マッチングに関するリアルタイムの概要を提供する。

 

インストール

マニュアル通りcondaの仮想環境を作ってテストした(ubuntu18.04LTS)。

依存

Github

#condaを使う。python>=3.6
conda create -n artic-rampart -y nodejs=12
conda activate artic-rampart
conda install -y artic-network::rampart=1.1.0

#dependency
conda install -y anaconda::biopython
conda install -y -c conda-forge -c bioconda snakemake=5.10.0 #snakemake<5.11
conda install -y bioconda::minimap2=2.17

#optional (if you require RAMPART to perform demuxing)
python -m pip install git+https://github.com/artic-network/Porechop.git@v0.3.2pre

python -m pip install binlorry==1.3.0_alpha1

rampart --help

$ rampart --help

usage: rampart [-h] [-v] [--verbose] [--ports PORTS PORTS]

               [--protocol PROTOCOL] [--title TITLE]

               [--basecalledPath BASECALLEDPATH]

               [--annotatedPath ANNOTATEDPATH]

               [--referencesPath REFERENCESPATH]

               [--referencesLabel REFERENCESLABEL]

               [--barcodeNames BARCODENAMES [BARCODENAMES ...]]

               [--annotationOptions ANNOTATIONOPTIONS [ANNOTATIONOPTIONS ...]]

               [--clearAnnotated] [--simulateRealTime SIMULATEREALTIME]

               [--devClient] [--mockFailures]

               

 

RAMPART v1.1.0: Read Assignment, Mapping, and Phylogenetic Analysis in Real 

Time

 

Optional arguments:

  -h, --help            Show this help message and exit.

  -v, --version         Show program's version number and exit.

  --verbose             verbose output

  --ports PORTS PORTS   The ports to talk to the client over. First: client 

                        delivery, i.e. what localhost port to access rampart 

                        via (default: 3000). Second: socket to transfer data 

                        over (default: 3001)

  --protocol PROTOCOL   path to a directory containing protocol config files

 

Config commands:

  Override options from config files

 

  --title TITLE         experiment title

  --basecalledPath BASECALLEDPATH

                        path to basecalled FASTQ directory (default: don't 

                        annotate FASTQs)

  --annotatedPath ANNOTATEDPATH

                        path to destination directory for annotation CSVs - 

                        will be created if it doesn't exist (default: '.

                        /annotations')

  --referencesPath REFERENCESPATH

                        path to a FASTA file containing a panel of reference 

                        sequences

  --referencesLabel REFERENCESLABEL

                        the reference header field to use as a reference 

                        label (if not just the reference name)

  --barcodeNames BARCODENAMES [BARCODENAMES ...]

                        specify mapping of barcodes to sample names - e.g. 

                        'BC01=Sample1' (can have more than one barcode 

                        mapping to the same name)

  --annotationOptions ANNOTATIONOPTIONS [ANNOTATIONOPTIONS ...]

                        pass through config options to the annotation script 

                        (key=value pairs)

 

Runtime commands:

  Options to specify how RAMPART behaves

 

  --clearAnnotated      remove any annotation files present when RAMPART 

                        starts up (force re-annotation of all FASTQs)

  --simulateRealTime SIMULATEREALTIME

                        simulate real-time annotation with given delay 

                        between files (default none)

 

Development commands:

  --devClient           don't serve build (client)

  --mockFailures        stochastic failures (annotating / parsing)

 

 

 テストラン

GitHubレポジトリにエボラウイルス(EBOV)のプロトコルの例がある。これを使う。

git clone https://github.com/artic-network/rampart.git
cd rampart/example_data/20181008_1405_EBOV

中身

f:id:kazumaxneo:20200711233345p:plain

MinIONのランではfastq/passにbasecallが上手くいったfastqが追加されていくが、RAMPARTはこのfastq ファイル群を認識する。

 

バーコードとサンプル名のマッピングなど実行のためのパラメータを定義したrun_configuration.jsonファイルが認識される。(ここではファイルは揃っているが、実際のランではリファレンスゲノムも用意する必要がある=>レポジトリの例を確認)。

rampart --protocol ../../example_protocols/EBOV --clearAnnotated

--clearAnnotated      remove any annotation files present when RAMPART starts up (force re-annotation of all FASTQs)

 

localhost:3000 にアクセスする。 

f:id:kazumaxneo:20200711230256p:plain

 

クリックすると展開される。ウィルスゲノム全体のリードのカバー率。

f:id:kazumaxneo:20200711230517p:plain

 

累積

f:id:kazumaxneo:20200711235125p:plain

 

バーコードにより複数データをmultiplexing runしている場合、 バーコードごとに出力を確認できる。

JSONファイルで遺伝子ポジションを指定している場合、遺伝子ごとのシーケンス率を表示できる。

f:id:kazumaxneo:20200711235138p:plain

f:id:kazumaxneo:20200711230407p:plain

f:id:kazumaxneo:20200711230413p:plain

f:id:kazumaxneo:20200711230418p:plain

f:id:kazumaxneo:20200711230422p:plain

f:id:kazumaxneo:20200711230426p:plain

 

実際のセットアップ例はmanualを読んでください。jsonファイルを編集して使用します。

https://github.com/artic-network/rampart/blob/master/docs/setting-up.md

引用

https://github.com/artic-network/rampart

 

2倍体ゲノムアセンブリのHaplotigsを追い出してPrimary contigsを出力する Purge Haplotigs

 2020 7/11 図追加

2020 7/13  タイトル修正

 

 第三世代の1分子シーケンシングにおける最近の進歩は、非常に高いレベルの連続性と完全性を持つde novoゲノムアセンブリを可能にした。さらに、最近の「diploid aware」なゲノムアセンブラの進歩により、高度にヘテロ接合した二倍体ゲノムアセンブラの品質が大幅に向上した。FALCONやCanuなどの二倍体対応アセンブラは、二倍体ゲノムのハプロタイプを融合した表現を生成することができる。また、FALCON UnzipやSupernovaなどのアセンブラの中には、さらに進化して、両方の親対立遺伝子を別々に表現した大規模なフェーズブロックを生成するものもある。FALCON Unzipアセンブリでは、アセンブリグラフ上でフェーズングが行われ、「一次コンティグ」(ハプロイドアセンブリ)と関連する「haplotigs」が生成され、これらの一次コンティグと二次ハプロチグのunionからなる二倍体アセンブリが生成される。

 理想的なハプロイド表現(プライマリコンティグ)は、2つのハプロームのすべてのヘテロ接合領域と、両方のハプロームのすべてのヘミ接合領域の1つの対立遺伝子コピーで構成されている。これにより、どちらかのハプロームに含まれる領域が、ハプロイド表現の一つの位置に完全にアラインすることが保証される。セカンダリハプロチグには、両方のハプロームに見られるヘテロ接合領域の2つの対立遺伝子コピーのうちの1つが含まれていなければならない;この点で、ハプロチグはハプロイド表現の段階的な情報として機能する。

 非常に高いヘテロ接合性を持つ領域は、de novoゲノムアセンブリではまだ問題がある。このような状況では、一組の対立遺伝子配列がヌクレオチド多様性のある閾値を超えると、ほとんどのアルゴリズムは、期待される単一のハプロタイプ融合コンティグではなく、これらの領域を別々のコンティグとしてアセンブルする 。結果、ハプロイドゲノムのサイズよりも著しく大きいアセンブリとなり、ハプロイドアセンブリにおけるこれらの対立遺伝子コンティグの存在は、下流の解析において問題となる。二倍体アセンブリを作製する場合、両方の対立遺伝子が存在する可能性があるが、対立遺伝子コンティグのペアを特定するための手順が依然として必要である。 

 この問題に対処するためにいくつかのツールが試みられてきた。HaploMerger2ツールキットとRedundansアセンブリパイプラインは、ショートリード配列からハプロタイプ融合アセンブリを生成するように設計されている。しかし、カバレッジの深さを考慮せずにコンティグ同士のアラインメントのみに基づいてコンティグを自動的に除去すると、反復的なコンティグやパラロジカルなコンティグが過剰にパージ(排除)されてしまう可能性がある。さらに、ハプロタイプ配列を分解して段階的なアセンブリを作成することが有利であることが証明されている。ロングリードのアセンブリで使用できるスクリプトには、配列のアラインメントを使用してホモログを特定し、手動でのキュレーションを支援する get_homologs.pyや、遺伝子アノテーションを使用して同調領域を一致させる HomolContigsByAnnotation などがある。それぞれに独自の長所と短所があるが、どちらもコンティグの再配置をユーザーが手動で行わなければならないという問題がある。

 本研究の目的は、1分子ロングリードシーケンシング技術を用いて作製されたアセンブリ内の対立遺伝子コンティグを迅速かつ自動的に特定し、再割り当てできる新しいパイプラインを開発することであった。Purge Haplotigsはインストールが簡単で、必要なコマンドは3つだけである。このパイプラインは、重複排除ハプロイドアセンブリを生成するためのハプロイドアセンブリ、または、改良された重複排除プライマリハプロイドアセンブリとより完全なセカンダリハプロティグアセンブリを生成するための二倍体アセンブリのいずれかで動作する。最後に、パイプラインはまた、必要に応じてアセンブリの手動検査とキュレーションを支援するように設計されたいくつかの出力を生成する。

 

現在のところ、Purge Haplotigsは二倍体ゲノムに対してのみテストされている。

 

 

f:id:kazumaxneo:20200710173653p:plain

Flow-chart for the Purge Haplotigs pipeline. 論文より転載

 

HaplotigsとPrimary contigs の区別。Pacbioのpdf資料を転載。

f:id:kazumaxneo:20200710195219p:plain

 

Twitter

https://twitter.com/search?q=Purge%20Haplotigs&src=typed_query

 

インストール

ubuntu18.04でcondaの仮想環境を作ってテストした。

 依存

  • Bash
  • BEDTools (tested with v2.26.0)
  • SAMTools (tested with v1.7)
  • Minimap2 (tested with v2.11/v2.12, https://github.com/lh3/minimap2)
  • Perl (with core modules: FindBin, Getopt::Long, Time::Piece, threads, Thread::Semaphore, Thread::Queue, List::Util)
  • Rscript (with ggplot2)

conda create -n purge_haplotigs_env -y
conda activate purge_haplotigs_env
conda install -c bioconda purge_haplotigs -y

#インストール チェック

> purge_haplotigs test

 

[10-07-2020 17:12:19] Writing contig associations

[10-07-2020 17:12:19] Writing the reassignment table and new assembly files

[10-07-2020 17:12:19]

 

PURGE HAPLOTIGS HAS COMPLETED SUCCESSFULLY!

 

samtools faidx curated.fasta

samtools faidx curated.haplotigs.fasta

purge_haplotigs ncbiplace -p curated.fasta -h curated.haplotigs.fasta -t 4 -f

[10-07-2020 17:12:19] minimap2 OK!

[10-07-2020 17:12:19] samtools OK!

[10-07-2020 17:12:19]

Beginning Pipeline

 

PARAMETERS:

Primary contigs FASTA       curated.fasta

Haplotigs FASTA             curated.haplotigs.fasta

Out FASTA                   ncbi_placements.tsv

Threads                     4

Coverage align cutoff       50 %

 

RUNNING USING COMMAND:

purge_haplotigs place -p curated.fasta -h curated.haplotigs.fasta -t 4 -f

 

[10-07-2020 17:12:19] Reading contig lengths

[10-07-2020 17:12:19] Running minimap2 hit search

[10-07-2020 17:12:19] Minimap2 hit search done

[10-07-2020 17:12:19] Reading minimap2 alignments

[10-07-2020 17:12:19] Getting best hits

[10-07-2020 17:12:19] Getting alignments coords for best hits

[10-07-2020 17:12:19] Renaming contigs and haplotigs in the FALCON Unzip format

[10-07-2020 17:12:19] Writing placement file

[10-07-2020 17:12:19] Done!

md5sum -c validate.md5

make: md5sum: No such file or directory

make: *** [Makefile:75: chk_out_files] Error 127

 

ONE OR MORE TESTS FAILED

 

#help 

purge_haplotigs -h

$ purge_haplotigs -h

 

ERROR: no such sub-command -h

 

USAGE:

purge_haplotigs  <command>  [options]

 

COMMANDS:

-- Purge Haplotigs pipeline:

    hist            First step, generate a read-depth histogram for the genome

    cov             Second step, get contig coverage stats and flag 'suspect' contigs

    purge           Third step, identify and reassign haplotigs

 

-- Other scripts

    clip            EXPERIMENTAL; Find and clip overlapping contig ends

    place           Generate a placement file for submission to NCBI

    test            Test everything!

 

purge_haplotigs cov -h

$ purge_haplotigs cov -h

Option high requires an argument

 

USAGE:

purge_haplotigs  cov  -i aligned.bam.genecov  -l <integer>  -m <integer>  -h <integer>  [-o coverage_stats.csv -j 80  -s 80 ]

 

REQUIRED:

-i / -in        The bedtools genomecov output that was produced from 'purge_haplotigs readhist'

-l / -low       The read depth low cutoff (use the histogram to eyeball these cutoffs)

-h / -high      The read depth high cutoff

-m / -mid       The low point between the haploid and diploid peaks

 

OPTIONAL:

-o / -out       Choose an output file name (CSV format, DEFAULT = coverage_stats.csv)

-j / -junk      Auto-assign contig as "j" (junk) if this percentage or greater of the contig is 

                low/high coverage (DEFAULT = 80, > 100 = don't junk anything)

-s / -suspect   Auto-assign contig as "s" (suspected haplotig) if this percentage or less of the

                contig is diploid level of coverage (DEFAULT = 80)

> purge_haplotigs hist

USAGE:

purge_haplotigs  hist  -b aligned.bam  -g genome.fasta  [ -t threads ]

 

REQUIRED:

-b / -bam       BAM file of aligned and sorted reads/subreads to the reference

-g / -genome    Reference FASTA for the BAM file.

 

OPTIONAL:

-t / -threads   Number of worker threads to use, DEFAULT = 4, MINIMUM = 2

-d / -depth     Maximum cutoff for depth. DEFAULT = 200, increase if needed,

                set much higher than your expected average coverage.

> purge_haplotigs purge -h

USAGE:

purge_haplotigs  purge  -g genome.fasta  -c coverage_stats.csv

 

REQUIRED:

-g / -genome        Genome assembly in fasta format. Needs to be indexed with samtools faidx.

-c / -coverage      Contig by contig coverage stats csv file from the previous step.

 

OPTIONAL:

-t / -threads       Number of worker threads to use. DEFAULT = 4

-o / -outprefix     Prefix for the curated assembly. DEFAULT = "curated"

-r / -repeats       BED-format file of repeats to ignore during analysis.

-d / -dotplots      Generate dotplots for manual inspection.

-b / -bam           Samtools-indexed bam file of aligned and sorted reads/subreads to the

                    reference, required for generating dotplots.

-f / -falconNaming  Rename contigs in the style used by FALCON/FALCON-unzip

 

ADVANCED:

-a / -align_cov     Percent cutoff for identifying a contig as a haplotig. DEFAULT = 70

-m / -max_match     Percent cutoff for identifying repetitive contigs. Ignored when 

                    using repeat annotations (-repeats). DEFAULT = 250

-I                  Minimap2 indexing, drop minimisers every N bases, DEFAULT = 4G

-v / -verbose       Print EVERYTHING.

-limit_io           Limit for I/O intensive jobs. DEFAULT = -threads

-wind_min           Min window size for BED coverage plots (for dotplots). DEFAULT = 5000

-wind_nmax          Max windows per contig for BED coverage plots (for dotplots). DEFAULT = 200

 

 

 

実行方法

ランにはドラフトアセンブリ配列が必要。入力ドラフトアセンブリ配列として、ハプロイドアセンブリ(FALCONやCANUなど)、または二倍体アセンブリ(FALCON Unzipなど)のいずれかを使用できる。また、Pacbioのロングリードシークエンシングデータが必要。ショートリードでも動作するが、ロングリードに最適化されている。

 

1、マッピング

PacBioのサブリード、またはショートリードをハプロイドまたは二倍体ゲノムアセンブリマッピングしてbamを作成する。

minimap2 -t 12 -ax map-pb genome.fa subreads.fasta.gz --secondary=no \
| samtools sort -m 1G -@ 12 -o aligned.bam -T tmp.ali

 

2、カバレッジヒストグラムの生成

purge_haplotigs hist -b aligned.bam -g genome.fasta -t 12
  • -b    BAM file of aligned and sorted reads/subreads to the reference
  • -g    Reference FASTA for the BAM file
  • -t    Number of worker threads to use, DEFAULT = 4, MINIMUM = 2

BEDToolsのgenomecovファイルと、カバレッジプロファイルを示した画像ファイルが出力される。Collapsed haplotype のコンティグの場合は、両方の対立遺伝子からのリードがマップされるが、対立遺伝子が別々のコンティグとしてアセンブルされた場合は、リードが2つのコンティグに分割され、リードデプスが半分になる。これを利用してハプロチグである可能性の高いコンティグを識別する。

f:id:kazumaxneo:20200711231701p:plain

ハプロイドアセンブリでは、重複が発生している場合は二峰性の分布が観察される。0.5×の深さのピークは重複している領域から、1×の深さのピークはハプロタイプが適切に融合している領域からのものである。二倍体ではアセンブリ全体が複製されている必要があるため、1×のピークは非常に小さいか、または全く見えないかもしれない(論文より)。

 

マニュアルには、phased assemblyである場合は二峰性のピークは二倍体のピークは非常に小さい可能性があるとの記載がある。マニュアルに出力例あり。


3、コンティグごとのカバレッジの分析

コンティグごとにカバレッジを分析する。 前のステップで得られたカバレッジプロファイルからカバレッジカットオフ値を指定する。

purge_haplotigs cov -i aligned.bam.genecov -l <integer> -m <integer> -h <integer> \
-o coverage_stats.csv -j 80 -s 80
  • -i     The bedtools genomecov output that was produced from 'purge_haplotigs readhist'
  • -l      The read depth low cutoff (use the histogram to eyeball these cutoffs)
  • -h     The read depth high cutoff
  • -m    The low point between the haploid and diploid peaks
  • -o     Choose an output file name (CSV format, DEFAULT = coverage_stats.csv)
  • -j      Auto-assign contig as "j" (junk) if this percentage or greater of the contig is  low/high coverage (DEFAULT = 80, > 100 = don't junk anything)
  • -s    Auto-assign contig as "s" (suspected haplotig) if this percentage or less of the contig is diploid level of coverage (DEFAULT = 80)

コンティグのカバレッジ統計のcsvファイルが作成され、疑わしいコンティグ(0.5x)には更なる分析や除去のためのフラグが立てられる (1xのプライマリなコンティグにはフラグは立てられない(論文図2B))。二倍体アセンブリの場合、両方のハプロタイプが存在しているはずなので、ほとんどのコンティグは解析のためにフラグが立てられ、以後の解析でプライマリなコンティグに選抜されると推測される。

 

4、パージパイプラインの実行

フラグが立てられたコンティグは、コンパニオンコンティグとのシンテニーの同定を試みるために、Minimap2による配列のアラインメントが行われ、どのコンティグを維持するかが評価される(論文図2C)。アライメントスコアがカットオフ(デフォルトでは70%以上)を超えるコンティグは、ハプロティグとして再割り当ての対象となる。 

purge_haplotigs purge -g genome.fasta -c coverage_stats.csv
  • -g    Genome assembly in fasta format. Needs to be indexed with samtools faidx.
  • -c    Contig by contig coverage stats csv file from the previous step.

ハプロチグが入れ子になっていたり、重複していたり、あるいはほとんどが反復配列で構成されている場合には、コンティグの競合が発生する可能性がある。そのような場合、コンティグの同定、アラインメントスコアリング、コンフリクトの解決、コンティグの再割り当てのステップは条件を満たすコンティグがなくなるまで繰り返し行われる。

 

3つのFASTA形式のファイル;キュレーションされたコンティグ、ハプロティグとして再割り当てされたコンティグ、アーティファクトとして再割り当てされた異常なカバレッジのコンティグ,が出力される。

f:id:kazumaxneo:20200710193908p:plain

  1. <prefix>.fasta.キュレーションされたプライマリコンティグ
  2. <prefix>.haplotigs.fasta. 最初の入力アセンブリで識別されたすべてのハプロチグ
  3. <prefix>.artefacts.fasta. カバレッジが非常に低い/高いコンティグ(注意: ミトコンドリア/葉緑体/その他のコンティグもこの中に存在する可能性が高い)。

 

 

論文にもあるように、結果はBUSCOのduplicated BUSCOなどを見ることで評価できます。

 

追記

purge_haplotigs clip(ベータ)を走らせることで、末端に残ってしまっているオーバーラップをクリッピングできる()。

引用

Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies

Michael J. Roach, Simon A. Schmidt & Anthony R. Borneman
BMC Bioinformatics volume 19, Article number: 460 (2018)

 

参考

https://www.pacb.com/wp-content/uploads/4May2017_JasonChin_ChalengesDiploidAssembly.pdf

Linnean分類システムのランクに応じて分類学の系統を提供する分類学データベース Taxallnomy

 

 あらゆる生物学的データは分類学的データと密接にリンクしており、いくつかのバイオインフォマティクス分析は目的を達成するために分類学的情報に依存している。メタゲノミクス、臨床法医学、その他の分野では、サンプル中に存在する生物を同定し、グループ化するために、完全に注釈された分類学的データに依存しており、多くの場合、結果をfamily、order、class、またはphylumなどの分類学的ランクにまとめている。さらに、進化論的な分析では、これまでに提案された分類学的分類に基づいて議論が行われている。分類学的情報はいくつかの分類学的データベースから得ることができ、例えば「生命カタログ」は、「Tree of Life」、「Encyclopedia of Life」、「GBIF」などの他のプロジェクトに分類学的バックボーンを提供している。これらのデータベースで提供される情報は、FishBase、AmphibiaWeb、AnimalBaseなどのように、より特定のクレードをカバーする他のデータベースに供給する分類学の専門家によってサポートされている。しかし、分子配列を含む解析は、INSDCを構成するデータベースのいずれかにDNAやタンパク質配列が登録されている生物の分類学名や系統が膨大にまとめられた参照分類学データベースであるNCBI taxonomyに依存している。INSDCは、GenBank、ENA、DBJJの3つの主要な分子配列リポジトリから構成されているため、INSDCのデータを利用しているUniprotKB、Ensembl、Pfam、SMART、Panther、OMA、miRBaseなど、多様なテーマをカバーする生物データベースでは、NCBI taxonomyの情報が広く利用されている。また、PDB、ArrayExpress、KEGGのような他の主要な生物学的一次データベースもNCBI分類学データベースの分類学データとリンクしており、バイオインフォマティクス分野におけるこのデータベースの貢献は否定できないことを示している。
 NCBI taxonomyを構成する分類学的分類は、分類学的および分子系統的文献の見解を反映したトポロジーを持つ系統分類スキームに従っている。ツリーの各ノードは分類群を表し、各ノードには分類学的名称と分類学的識別子(txid)が付与されている。さらに、いくつかのノードは、分類学的ランクを持っている場合があり、これはリンネの分類システムで使用されているものに似ている。いくつかのバイオインフォマティクスのアプローチは、例えば、メタゲノムデータの分類学的プロファイルを作成したり、配列データの分類学的分類を支援するために、NCBI分類学が提供するランクベースの分類に依存している。しかし、バイオインフォマティクスのコミュニティでは、ランク情報が広く使われていることに加えて、これらのデータを管理する際に考慮すべき重要な問題がいくつかある。いくつかの生物の系統を検索する際に、いくつかのランクが欠けていることが観察される。2019年5月に行われたNCBI taxonomyに関するコンサルテーションでは、例えばブタ(Sus scrofa, NCBI:txid9823)の系統には、Orderランクを持つタクソンがなかった。一方、Thale cress(Arabidopsis thaliana, NCBI:txid3702)は、その系統にOrder rankを持つ分類群を含んでいなかった。さらに分類学的系譜を調べてみると、「ランクなし」と表記されているランクを持たない分類群も見つけることができた。これらは、単系統群を指摘しながら、系統情報を分類学のベースに加えている。
 これらの問題は、このグループの分類に関する専門家の間での不確実性や対立に起因している可能性があり、NCBI taxonomyの階層的ランクを不完全なものにしている。そのため、「このデータにはクラスランクの異なるいくつの分類群が含まれているのか」というような、分類学的ランクに関する単純な問い合わせが困難になる可能性がある。例えば、オオバコのクラスと、割り当てられていないモノコットのクラスがいくつか存在する場合、それらはすべて計算データベースでは "NULL "としてカウントされ、無関係なカウントをグループ化してしまう。このような解析のためには、分類学的ランクを組み込んだ階層的に完全な分類学ツリーが非常に有用である。そこで本研究では、NCBI Taxonomyが提供する分類学ツリーを用いて、すべての系統が同じ深さを持ち、すべての階層レベルが分類学ランクに対応する階層的な分類学ツリーを生成するアルゴリズムを開発した。最終的なデータベースは、ツリーを構成する系統のすべての分類学的ランクの分類学的名称を提供することから、taxallnomyと名付けられた。ユーザーは、bioinfo.icb.ufmg.br/taxallnomyのウェブサイトから、taxallnomyデータベースの階層構造にアクセスし、探索することができる。APIを介してプログラムでデータにアクセスしたり、ローカルマシンでtaxallnomyデータベースを作成するための手順もtaxallnomyのウェブサイトで提供されている。ローカルでの作成は非常に簡単で、更新された情報の使用を許可している。

 新しい階層構造はtaxallnomyと名付けられ、現在NCBI Taxonomyデータベースで使用されている33のTaxonomic rankに対応する33の階層レベルを含んでいる。Taxallnomyから、ユーザーはNCBI Taxonomyデータベースで利用可能なすべての分類の33のノードを持つ完全な分類学的系統を得ることができる。

 

Github


webサービス

http://biodados.icb.ufmg.br/taxallnomy/にアクセスする。

f:id:kazumaxneo:20200624010439p:plain

 

f:id:kazumaxneo:20200709232718p:plain


ここではE.coliの分類を調べてみる。E.coliのtxidである562を入力。

f:id:kazumaxneo:20200709232820p:plain


結果

f:id:kazumaxneo:20200709232836p:plain


ノードが小さい。上の編集バーをスライドさせてフォントサイズを拡大した。またノード間の幅(エッジ)を狭めた。

f:id:kazumaxneo:20200709233547p:plain

 

common rankからmain rankに変更。

f:id:kazumaxneo:20200709233716p:plain

 

 

上のcommon rank状態の階級をNCBI taxnomyの分類と比較する。

https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=562&lvl=3&p=has_linkout&p=blast_url&p=genome_blast&lin=f&keep=1&srchmode=1&unlock

f:id:kazumaxneo:20200709233858p:plain

 

NCBI taxnomyのE.coli LIneageはcellular organisms; Bacteria; Proteobacteria; Gammaproteobacteria; Enterobacterales; Enterobacteriaceae; Escherichiaで、Taxallnomyもmain rank表示ではその通り表示されている。

f:id:kazumaxneo:20200709233716p:plain

 

再びcommon rankに戻すと、Taxallnomyアルゴリズムで作成されたユニークなノードが追加された。

f:id:kazumaxneo:20200709234359p:plain

 subgenusやsubfammilyなどの分類階級が追加されている。

 

taxonノードの色は以下の情報を表す。Taxallnomyアルゴリズムで作成されたユニークなノードはになる。

f:id:kazumaxneo:20200709232233p:plain

 

結果は様々なフォーマットでダウンロードできる。

引用

Taxallnomy: Closing gaps in the NCBI Taxonomy

Sakamoto, Tetsu​, Ortega, J. Miguel​

bioRxiv, posted May 30, 2020

Vertebrate Genomes Projectで使用されているマイトゲノムアセンブリパイプライン mitoVGP

 

 最新のシーケンス技術により、比較的小さなミトコンドリアゲノムのアセンブリは容易に行えるようになるはずである。しかし、ミトコンドリアアセンブリに直接対応したツールはほとんど存在しない。脊椎動物ゲノムプロジェクト(VGP)の一環として、著者らはミトコンドリアリードの類似性に基づいた同定とミトコンドリアゲノムのde novoアセンブリのための完全自動化パイプラインであるmitoVGPを開発した。本パイプラインは、100種の脊椎動物の完全なマイトゲノムアセンブリに成功した。組織タイプとライブラリサイズの選択がマイトゲノムの配列決定とアセンブリに大きな影響を与えることが分かった。mitoVGPのアセンブリをショートリードシークエンシングに基づいた完全と思われるリファレンスマイトゲノムと比較したところ、参照マイトゲノムのエラー、欠落した配列、不完全な遺伝子、特にリピート領域の遺伝子、が同定された。またmitoVGPのアセンブリは、新しい遺伝子領域の重複を同定し、ミトコンドリアゲノムの進化と組織についての新たな洞察を明らかにした。

 

VGP


インストール

Github

#pacbioのリード用のパイプラインとONTのリード用のパイプラインがある。
#2つの依存するツールのpythonコードが2と3に別れているため、使用時にそれぞれ個別の過疎環境をactivateして使用する様に設計されている。
git clone https://github.com/gf777/mitoVGP.git
cd mitoVGP/
tar -xvf canu-1.8.Linux-amd64.tar.xz
rm canu-1.8.Linux-amd64.tar.xz

#pacbioとONTそれぞれの仮想環境を作成
#pacbio
conda env create -f mitoVGP_conda_env_pacbio.yml
#ONT
conda env create -f mitoVGP_conda_env_ONT.yml

./mitoVGP -h

$ ./mitoVGP -h

 

Usage: './mitoVGP -s species -i species_ID -r reference -t threads'

 

mitoVGP is used for reference-guided de novo mitogenome assembly using a combination of long and short read data.

 

An existing reference from closely to distantly related species is used to identify mito-like reads in pacbio WGS data,

which are then employed in de novo genome assembly. The assembly is further polished using both long and short read data,

and linearized to start with the conventional Phenylalanine tRNA sequence.

 

Check the github page https://github.com/GiulioF1/mitoVGP for a description of the pipeline.

A complete Conda environment with all dependencies is available to run the pipeline in the same github page.

 

This script a simple wrapper of the scripts found in the scripts/ folder. You can find more information

on each step in the help (-h) of each script.

 

Required arguments are:

-a long read sequencing platform (Pacbio/ONT)

-s the species name (e.g. Calypte_anna)

-i the VGP species ID (e.g. bCalAnn1)

-r the reference sequence fasta file

-t the number of threads

 

Optional arguments are:

-g the putative mitogenome size (potentially, that of the reference genome). If not provided, length of reference is used.

  It does not need to be precise. Accepts Canu formatting.

-d multithreaded download of files (true/false default: false) !! caution: true may require considerable amount of space.

-1 use pacbio/nanopore reads from list of files, requires absolute path (default looks into aws)

-2 use PE illumina reads from list of files (with fw and rv reads including R1/R2 in their names), requires absolute path (default looks into aws)

-m the aligner (blasr|minimap2|pbmm2). Default is pbmm2

-f filter reads by size prior to assembly (reduces the number of NUMT reads and helps the assembly)

-p filter reads by percent coverage of the reference over their length (avoid noise in the assembly when low coverage)

-o the options for Canu

-v picard validation stringency (STRICT/LENIENT default: STRICT)

-z increase sensitivity of mummer overlap detection

-b use gcpp or variantCaller during arrow polishing for 2.0 or earlier chemistry respectively (gcpp/variantCaller default: gcpp)

 

 

 

実行方法 

pacbio。pacbio のchemistryバージョンによって"-b"で指定するバリアントコーラーは変更する。pacbio chemistry 2.0では"gcpp"、2.0以下では"variantCaller"、RSIIでは"blasr"を選ぶことも出来る。

conda activate mitoVGP_pacbio
./mitoVGP -a pacbio -s <species_name> -i <VGP_species_ID> -r ref_mtDNA.fasta -t 12 -b variantCaller
  • -a   long read sequencing platform (Pacbio/ONT)
  • -s   the species name (e.g. Calypte_anna)
  • -i    the VGP species ID (e.g. bCalAnn1)
  • -r    the reference sequence fasta file
  • -b   use gcpp or variantCaller during arrow polishing for 2.0 or earlier chemistry respectively (gcpp/variantCaller default: gcpp)

 

ONT

conda activate mitoVGP_ONT
./mitoVGP -a ONT -s <species_name> -i fMasArm1 -r ref_mtDNA.fasta -t 12 -b variantCaller

 

mitoVGPによるマイトゲノムアセンブルGenomeArkのいくつかの種でも利用されている(GenomeArkはVEPによって作成されたゲノムとアノテーション情報のデータベース)。

引用

Complete vertebrate mitogenomes reveal widespread gene duplications and repeats

Giulio Formenti, Arang Rhie, Jennifer Balacco, Bettina Haase, Jacquelyn Mountcastle, Olivier Fedrigo, Samara Brown, Marco Capodiferro, Farooq O. Al-Ajli, Roberto Ambrosini, Peter Houde, Sergey Koren, Karen Oliver, Michelle Smith, Jason Skelton, Emma Betteridge, Jale Dolucan, Craig Corton, Iliana Bista, James Torrance, Alan Tracey, Jonathan Wood, Marcela Uliano-Silva, Kerstin Howe, Shane McCarthy, Sylke Winkler, Woori Kwak, Jonas Korlach, Arkarachai Fungtammasan, Daniel Fordham, Vania Costa, Simon Mayes, Matteo Chiara, David S. Horner, Eugene Myers, Richard Durbin, Alessandro Achilli, Edward L. Braun, Adam M. Phillippy, Erich D. Jarvis, The Vertebrate Genomes Project Consortium

bioRxiv, Posted July 01, 2020

 

関連


 

ヌクレオチド配列をアセンブリグラフにアラインメントする SPAligner

 

 ゲノムアセンブリのグラフベースの表現は、最近では遺伝子検索からハプロタイプ分離まで、さまざまなアプリケーションで利用されている。これらのアプリケーションの多くは、アセンブリグラフへの配列のアラインメントに基づいているが、このようなアラインメントを見つけるための既存のソフトウェアツールには重要な制限がある。

 長く分岐した配列をアセンブリグラフにアラインメントするための新しいツールSPAlignerを紹介し、SPAlignerが第三世代シーケンシングデータをマッピングするための効率的なソリューションであること、また複雑なメタゲノムデータセットの既知遺伝子の同定を容易にすることを実証した。

 本研究では、配列からゲノムへのアセンブリアラインメント問題を解決するためのグラフベースのアプローチの開発を促進する。SPAligner は SPAdes ツールライブラリの一部として実装されており、https://github.com/ablab/spades/archive/spaligner-paper.zip から入手可能である。

 

インストール

Github

README

#early accessバージョンもあるが、ここでは3.14.0を入れる。
#3.14.0 linux
wget http://cab.spbu.ru/files/release3.14.0/SPAdes-3.14.0-Linux.tar.gz
tar -xzf SPAdes-3.14.0-Linux.tar.gz
cd SPAdes-3.14.0-Linux/bin/

#build from source
git clone https://github.com/ablab/spades.git
cd spades/assembler/
mkdir build && cd build && cmake ../src
make spaligner

./spaligner

$ ./spaligner 

ERROR: No input YAML was specified

ERROR: Sequence type is not provided (nanopore or pacbio)

ERROR: Path to file with sequences is not provided

ERROR: Path to file with graph is not provided

ERROR: k-mer value is not provided

SYNOPSIS

        ./spaligner <aligner parameters description (in YAML)> -d <value> -s <value> -g <value> -k <value> [-t <value>] [-o <dir>]

 

OPTIONS

        -d, --datatype <value>

                    type of sequences: nanopore, pacbio

 

        -s, --sequences <value>

                    path to fasta/fastq file with sequences

 

        -g, --graph <value>

                    path to GFA-file or SPAdes saves folder

 

        -k, --kmer <value>

                    graph k-mer size (odd value)

 

        -t, --threads <value>

                    # of threads to use

 

        -o, --outdir <dir>

                    output directory

 

 

 実行方法

configのyamlファイル、fastqとアセンブリグラフファイルを指定する。ロングリードシークエンシングデータの種類について、nanoporeかpacbioのどちらかを指定する必要がある。

spaligner spaligner_config.yaml -d nanopore -s input.fq -g assembly_graph.gfa -k <k-mer value> -t 12 -o output

様々な形式(tsv, fasta, GPAなど)でグラフへのアラインメント結果を出力できる。

 

出力フォーマットの詳細は、preprintと上にリンクを張ったGithub - READMEを読んでください。

引用

SPAligner: Alignment of Long Diverged Molecular Sequences to Assembly Graphs

Tatiana Dvorkina, Dmitry Antipov, Anton Korobeynikov, Sergey Nurk

bioRxiv, Posted August 23, 2019