macでインフォマティクス

macでインフォマティクス

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

cgMLSTコールから距離行列を出力する cgmlst-dists

 

 

タイトルの通りのツール。 ChewBBACA(Github)の出力に対応している。

 

インストール

Github

#bioconda (link)
conda install -c bioconda cgmlst-dists

#from source
git clone https://github.com/tseemann/cgmlst-dists.git
cd cgmlst-dists
make

cgmlst-dists

$ cgmlst-dists

SYNOPSIS

  Pairwise CG-MLST distance matrix from allele call tables

USAGE

  cgmlst-dists [options] chewbbaca.tab > distances.tsv

OPTIONS

  -h Show this help

  -v Print version and exit

  -q Quiet mode; do not print progress information

  -c Use comma instead of tab in output

  -m N Output: 1=lower-tri 2=upper-tri 3=full [3]

  -x N Stop calculating beyond this distance [9999]

URL

  https://github.com/tseemann/cgmlst-dists

 

実行方法

ChewBBACAからのcgMLST allele callのファイルを指定する。

git clone https://github.com/tseemann/cgmlst-dists.git
cd cgmlst-dists/
cgmlst-dists test/boring.tab > distances.tab

出力

f:id:kazumaxneo:20201125235715p:plain

正の整数ではないアレルコールはすべてゼロに変換される。距離はハミング距離(wiki)だが、ゼロを除いたものになる。

 

引用

GitHub - tseemann/cgmlst-dists: 🐻⇔🐨 Calculate distance matrix from ChewBBACA cgMLST allele call tables

 

 

infernal

 

 infernalは、入力として与えられた構造的にアノテーションされた複数の配列アラインメントから、共分散モデル(CM)と呼ばれるRNAファミリーの配列と二次構造の確率的プロファイルを構築する。infernalは、共分散モデルを使用して、配列データベース内の新しいファミリーメンバーを検索し、潜在的に大規模な多重配列アラインメントを作成する。Infernalのバージョン1.1では、加速プロファイル隠れマルコフモデル(HMM)法とHMMバンド共分散モデルアライメント法に基づいたRNAホモロジー検索のための新しいフィルタパイプラインが導入されている。これにより、従来のバージョンよりも∼100倍の速度と、フィルタリングされていない網羅的な共分散モデル検索よりも∼10,000倍の加速が可能となった。

 ソースコード、ドキュメント、ベンチマークhttp://infernal.janelia.org からダウンロードできる。InfernalはGNU GPLv3の下でライセンスされており、LinuxMac OS/Xを含むPOSIX準拠のオペレーティング・システムに自由に移植可能である。ドキュメントには、チュートリアル付きのユーザーガイド、ファイルフォーマットやユーザーオプションの説明、ソフトウェアに実装されている方法についての追加の詳細が含まれている。

 

HP

http://eddylab.org/infernal/

PDF manual

http://eddylab.org/infernal/Userguide.pdf

 

インストール

Github

#bioconda
conda install -c bioconda -y infernal

#homebrew
brew install infernal

#from source
wget eddylab.org/infernal/infernal-1.1.3.tar.gz
tar xf infernal-1.1.3.tar.gz
cd infernal-1.1.3
./configure
make
make check
make install

cmalign -h

# cmalign -h

# cmalign :: align sequences to a CM

# INFERNAL 1.1.2 (July 2016)

# Copyright (C) 2016 Howard Hughes Medical Institute.

# Freely distributed under a BSD open source license.

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Usage: cmalign [-options] <cmfile> <seqfile>

 

Basic options:

  -h     : show brief help on version and usage

  -o <f> : output the alignment to file <f>, not stdout

  -g     : configure CM for global alignment [default: local]

 

Options controlling alignment algorithm:

  --optacc   : use the Holmes/Durbin optimal accuracy algorithm  [default]

  --cyk      : use the CYK algorithm

  --sample   : sample alignment of each seq from posterior distribution

  --seed <n> : w/--sample, set RNG seed to <n> (if 0: one-time arbitrary seed)

  --notrunc  : do not use truncated alignment algorithm

  --sub      : build sub CM for columns b/t HMM predicted start/end points

 

Options controlling speed and memory requirements:

  --hbanded    : accelerate using CM plan 9 HMM derived bands  [default]

  --tau <x>    : set tail loss prob for HMM bands to <x>  [1e-7]  (1e-18<x<1)

  --mxsize <x> : set maximum allowable DP matrix size to <x> Mb  [1028.0]  (x>0.)

  --fixedtau   : do not adjust tau (tighten bands) until mx size is < limit

  --maxtau <x> : set max tau <x> when tightening HMM bands  [0.05]  (0<x<0.5)

  --nonbanded  : do not use HMM bands for faster alignment

  --small      : use small memory divide and conquer (d&c) algorithm

 

Optional output files:

  --sfile <f>  : dump alignment score information to file <f>

  --tfile <f>  : dump individual sequence parsetrees to file <f>

  --ifile <f>  : dump information on per-sequence inserts to file <f>

  --elfile <f> : dump information on per-sequence EL inserts to file <f>

 

Other options:

  --mapali <f>    : include alignment in file <f> (same ali that CM came from)

  --mapstr        : include structure (w/pknots) from <f> from --mapali <f>

  --noss          : cmbuild --noss option was used w/aln from --mapali <f>

  --informat <s>  : assert <seqfile> is in format <s>: no autodetection

  --outformat <s> : output alignment in format <s>  [Stockholm]

  --dnaout        : output alignment as DNA (not RNA) sequence data

  --noprob        : do not include posterior probabilities in the alignment

  --matchonly     : include only match columns in output alignment

  --ileaved       : force output in interleaved Stockholm format

  --regress <f>   : save regression test data to file <f>

  --verbose       : report extra information; mainly useful for debugging

  --cpu <n>       : number of parallel CPU workers to use for multithreads  (n>=0)

 

Sequence input formats:   FASTA, GenBank

Alignment output formats: Stockholm, Pfam, AFA (aligned FASTA), A2M, Clustal, PHYLIP

cmsearch -h

# cmsearch -h

# cmsearch :: search CM(s) against a sequence database

# INFERNAL 1.1.2 (July 2016)

# Copyright (C) 2016 Howard Hughes Medical Institute.

# Freely distributed under a BSD open source license.

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Usage: cmsearch [options] <cmfile> <seqdb>

 

Basic options:

  -h        : show brief help on version and usage

  -g        : configure CM for glocal alignment [default: local]

  -Z <x>    : set search space size in *Mb* to <x> for E-value calculations  (x>0)

  --devhelp : show list of otherwise hidden developer/expert options

 

Options directing output:

  -o <f>       : direct output to file <f>, not stdout

  -A <f>       : save multiple alignment of all significant hits to file <s>

  --tblout <f> : save parseable table of hits to file <s>

  --acc        : prefer accessions over names in output

  --noali      : don't output alignments, so output is smaller

  --notextw    : unlimit ASCII text output line width

  --textw <n>  : set max width of ASCII text output lines  [120]  (n>=120)

  --verbose    : report extra information; mainly useful for debugging

 

Options controlling reporting thresholds:

  -E <x> : report sequences <= this E-value threshold in output  [10.0]  (x>0)

  -T <x> : report sequences >= this score threshold in output

 

Options controlling inclusion (significance) thresholds:

  --incE <x> : consider sequences <= this E-value threshold as significant  [0.01]

  --incT <x> : consider sequences >= this score threshold as significant

 

Options controlling model-specific reporting thresholds:

  --cut_ga : use CM's GA gathering cutoffs as reporting thresholds

  --cut_nc : use CM's NC noise cutoffs as reporting thresholds

  --cut_tc : use CM's TC trusted cutoffs as reporting thresholds

 

Options controlling acceleration heuristics*:

  --max      : turn all heuristic filters off (slow)

  --nohmm    : skip all HMM filter stages, use only CM (slow)

  --mid      : skip first two HMM filter stages (SSV & Vit)

  --default  : default: run search space size-dependent pipeline  [default]

  --rfam     : set heuristic filters at Rfam-level (fast)

  --hmmonly  : use HMM only, don't use a CM at all

  --FZ <x>   : set filters to defaults used for a search space of size <x> Mb

  --Fmid <x> : with --mid, set P-value threshold for HMM stages to <x>  [0.02]

 

Other options*:

  --notrunc     : do not allow truncated hits at sequence termini

  --anytrunc    : allow full and truncated hits anywhere within sequences

  --nonull3     : turn off the NULL3 post hoc additional null model

  --mxsize <x>  : set max allowed alnment mx size to <x> Mb [df: autodetermined]

  --smxsize <x> : set max allowed size of search DP matrices to <x> Mb  [128.]

  --cyk         : use scanning CM CYK algorithm, not Inside in final stage

  --acyk        : align hits with CYK, not optimal accuracy

  --wcx <x>     : set W (expected max hit len) as <x> * cm->clen (model len)

  --toponly     : only search the top strand

  --bottomonly  : only search the bottom strand

  --tformat <s> : assert target <seqdb> is in format <s>: no autodetection

  --cpu <n>     : number of parallel CPU workers to use for multithreads

 

*Use --devhelp to show additional expert options.

cmstat -h

# cmstat -h

# cmstat :: display summary statistics for CMs

# INFERNAL 1.1.2 (July 2016)

# Copyright (C) 2016 Howard Hughes Medical Institute.

# Freely distributed under a BSD open source license.

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Usage: cmstat [-options] <cmfile>

 

Options:

  -h          : show brief help on version and usage

  -E <x>      : print bit scores that correspond to E-value threshold of <x>

  -P <x>      : print bit scores that correspond to E-value threshold of <x>

  -T <x>      : print E-values that correspond to bit score threshold of <x>

  -Z <x>      : set database size in *Mb* to <x> for E-value calculations  [10]

  --cut_ga    : print E-values that correspond to GA bit score thresholds

  --cut_nc    : print E-values that correspond to NC bit score thresholds

  --cut_tc    : print E-values that correspond to TC bit score thresholds

  --key <s>   : only print statistics for CM with name or accession <s>

  --hmmonly   : print filter HMM bit scores/E-values, not CM ones

  --nohmmonly : print CM bit scores/E-values, even for models with 0 basepairs

cmbuild -h

# cmbuild -h

# cmbuild :: covariance model construction from multiple sequence alignments

# INFERNAL 1.1.2 (July 2016)

# Copyright (C) 2016 Howard Hughes Medical Institute.

# Freely distributed under a BSD open source license.

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Usage: cmbuild [-options] <cmfile_out> <msafile>

 

Basic options:

  -h        : show brief help on version and usage

  -n <s>    : name the CM(s) <s>, (only if single aln in file)

  -F        : force; allow overwriting of <cmfile_out>

  -o <f>    : direct summary output to file <f>, not stdout

  -O <f>    : resave consensus/insert column annotated MSA to file <f>

  --devhelp : show list of otherwise hidden developer/expert options

 

Alternative model construction strategies:

  --fast        : assign cols w/ >= symfrac residues as consensus

  --hand        : use reference coordinate annotation to specify consensus

  --symfrac <x> : fraction of non-gaps to require in a consensus column [0..1]

  --noss        : ignore secondary structure annotation in input alignment

  --rsearch <f> : use RSEARCH parameterization with RIBOSUM matrix file <f>

 

Other model construction options*:

  --null <f>  : read null (random sequence) model from file <f>

  --prior <f> : read priors from file <f>

 

Alternative relative sequence weighting strategies:

  --wpb     : Henikoff position-based weights  [default]

  --wgsc    : Gerstein/Sonnhammer/Chothia tree weights

  --wnone   : don't do any relative weighting; set all to 1

  --wgiven  : use weights as given in MSA file

  --wblosum : Henikoff simple filter weights

  --wid <x> : for --wblosum: set identity cutoff  [0.62]  (0<=x<=1)

 

Alternative effective sequence weighting strategies:

  --eent        : adjust eff seq # to achieve relative entropy target  [default]

  --enone       : no effective seq # weighting: just use nseq

  --ere <x>     : for --eent: set CM target relative entropy to <x>

  --eset <x>    : set eff seq # for all models to <x>

  --eminseq <x> : for --eent: set minimum effective sequence number to <x>  [0.1]

  --ehmmre <x>  : for --eent: set minimum HMM relative entropy to <x>

  --esigma <x>  : for --eent: set sigma param to <x>  [45.0]

 

Options for HMM filter construction*:

  --p7ere <x> : for the filter p7 HMM, set minimum rel entropy/posn to <x>

  --p7ml      : define the filter p7 HMM as the ML p7 HMM

 

Options for HMM filter calibration*:

  --EmN <n>  : number of sampled seqs to use for p7 local MSV calibration  [200]

  --EvN <n>  : number of sampled seqs to use for p7 local Vit calibration  [200]

  --ElfN <n> : number of sampled seqs to use for p7 local Fwd calibration  [200]

  --EgfN <n> : number of sampled seqs to use for p7 glocal Fwd calibration  [200]

 

Options for refining the input alignment*:

  --refine <f> : refine input aln w/Expectation-Maximization, save to <f>

  -l           : w/--refine, configure model for local alignment [default: global]

  --gibbs      : w/--refine, use Gibbs sampling instead of EM

  --seed <n>   : w/--gibbs, set RNG seed to <n> (if 0: one-time arbitrary seed)

  --cyk        : w/--refine, use CYK instead of optimal accuracy

  --notrunc    : w/--refine, do not use truncated alignment algorithm

 

*Use --devhelp to show additional expert options.

 

 

テストラン

Infernalは、用意した複数の配列のアラインメントファイルから始まる。Stockholm format(Wiki)で、コンセンサス二次構造のアノテーションが含まれている必要がある(例:tutorial/tRNA5.sto)。

f:id:kazumaxneo:20201125113602p:plain

ストックホルム形式の二次構造アノテーションを用いたマルチプルRNA配列のアラインメントの例。ストックホルムフォーマットは、HMMERやInfernal、PfamやRfamデータベースで使用されているネイティブのアラインメントフォーマット。右の絵は参考のために掲載されているtRNAの二次構造。PDFマニュアルより転載。

 

git clone https://github.com/EddyRivasLab/infernal
cd infernal/
cp tutorial/tRNA5.c.cm tRNA5.cm
cmsearch tRNA5.cm tutorial/mrum-genome.fa

作成中

 

引用
Infernal 1.1: 100-fold faster RNA homology searches
Eric P. Nawrocki, Sean R. Eddy

Bioinformatics. 2013 Nov 15; 29(22): 2933–2935

 

関連


メタゲノムのシーケンシングリードからアセンブリしてCRISPRsを探す Crass

 

 Clustered Regularly Interspaced Short Palindromic Repeats (CRISPRs) は、バクテリオファージ(ファージ)から細胞を保護する細菌および古細菌の適応免疫システムを構成する。 CRISPR遺伝子座の分析により、ファージ感染の履歴が明らかになり、ファージとその宿主間の直接的なリンクが提供される。 CRISPR同定のための現在のツールはすべて完全なゲノムを分析するために開発されたものであり(論文執筆時点)、CRISPR遺伝子座がrepetitive structureと集団の不均一性のためにアセンブリが難しいメタゲノムデータセットの分析にはあまり適していない。ここでは、アセンブリやデータセット内のCRISPRの事前知識を必要とせずに、rawメタゲノムデータからCRISPR遺伝子座を識別および再構築するように設計された新しいアルゴリズムCrassを紹介する。アセンブリされたデータのCRISPRは、多くのコンティグ/スキャフォールドで断片化されていることが多く、CRISPR遺伝子座の母集団の不均一性を完全には表していない。 Crassは、アセンブリベースのアプローチを使用して以前に分析されたメタゲノムで、かなり多くのCRISPRを特定した。 Crassを使用して、システム内のファージと配列相同性のあるスペーサーを含むCRISPRを検出することができたが、これは他のアプローチでは特定できなかっただろう。 Crassの感度、特異性、速度の向上により、メタゲノムデータセット内のCRISPRの包括的な分析が容易になり、微生物コミュニティ内でのファージとホストの相互作用と共進化に関する理解が深まる。

 

インストール

ubuntu18.04のpython3.7環境でテストした。

Github

git cloen https://github.com/ctSkennerton/crass.git
cd crass/
./autogen.sh
./configure
make
make install

#bioconda (link)
conda install -c bioconda crass

crass

$ crass

Compiler Options:

RENDERING = 0

DEBUG = 0

MEMCHECK = 0

ASSEMBER = 1

VERBOSE_LOGGER = 0

Search Debugger =  0

 

Usage:  crass  [options] { inputFile ...}

 

General Options:

-h --help                    This help message

-l --logLevel        <INT>   Output a log file and set a log level [1 - 4]

-o --outDir          <DIR>   Output directory [default: .]

-V --version                 Program and version information

-g --logToScreen             Print the logging information to screen rather than a file

 

CRISPR Identification Options:

-d --minDR           <INT>   Minimim length of the direct repeat

                             to search for [Default: 23]

-D --maxDR           <INT>   Maximim length of the direct repeat

                             to search for [Default: 47]

-n --minNumRepeats   <INT>   Total number of direct repeats in a CRISPR for

                             it to be considered real [Default: 2]

-s --minSpacer       <INT>   Minimim length of the spacer to search for [Default: 26]

-S --maxSpacer       <INT>   Maximim length of the spacer to search for [Default: 50]

-w --windowLength    <INT>   The length of the search window. Can only be

                             a number between 6 - 9 [Default: 8]

CRISPR Assembly Options:

-f --covCutoff       <INT>   Remove groups with less than x spacers [Default: 3]

-k --kmerCount       <INT>   The number of the kmers that need to be

                             shared for clustering [Default: 6]

-K --graphNodeLen    <INT>   Length of the kmers used to make crispr nodes [Default: 7]

 

Output Options: 

-b --numBins          <INT>   The number of colour bins for the output graph.

                              Default is to have as many colours as there are

                              different values for the coverage of Nodes in the graph

-c --graphColour     <TYPE>   Defines the range of colours to use for the output graph

                              the different types available are:

                              red-blue, blue-red, green-red-blue, red-blue-green

-L --longDescription          Set if you want the spacer sequence printed along with the ID in the spacer graph. [Default: false]

-G --showSingltons            Set if you want to print singleton spacers in the spacer graph [Default: false]

 

実行方法

 fastqファイル(gzip圧縮にも対応)を指定する。

crass input.fq

f:id:kazumaxneo:20201124132348p:plain

出力

f:id:kazumaxneo:20201124132511p:plain

 

引用
Crass: identification and reconstruction of CRISPR from unassembled metagenomic data
Skennerton CT, Imelfort M, Tyson GW

Nucleic Acids Res. 2013 May 1;41(10):e105

 

関連

rawリードから探す場合


 

トランスポーターのデータベース TCDB

 

 膜輸送体は、細胞の分子組成やエネルギー状態を決定するチャネル、キャリア、ポンプ、group translocators、電子輸送体などの複雑なネットワークを形成する多様なタンパク質のグループを構成している(ref.1)。これらのタンパク質は、細胞内の全タンパク質の約10%を占め、栄養素、代謝の最終生成物、有害物質、高分子、シグナル分子、薬物、電子などをソースからシンクへと移動させ、化合物やエネルギー源を細胞内に取り込んだり、排出している(ref.2)。癌学、微生物病因学、ウイルス学の分野で特に重要なのは、薬剤流出ポンプが、病原性生物や癌細胞の薬剤耐性において支配的な役割を果たしていることである(ref.3,4)。世界中の数千の研究者が、細胞膜を横断する分子輸送の理解に貢献している(ref.5)。TCDB (tcdb.org)は、約10万人の異なる利用者を含む年間100万回以上の研究者によって利用されている。

 2001年6月、国際生化学・分子生物学連合(IUBMB)(wiki)は、地球の生物圏のすべての生物から得られた輸送タンパク質情報の組織化のための唯一の国際的に認められたシステムとして、トランスポーター分類(TC)システムを正式に採用した(ref.6,7)。メタゲノムシークエンシングの出現と計算生物学の最近の進歩(ここ5年)により、ナノバクテリア古細菌の数十種類もの新しい系統が発見されたため、TCDBはこれらの新しい系統でのみ発見されたトランスポーターの潜在的なファミリーや、既存のファミリーのメンバーと遠く離れた関係にあるタンパク質を含むように拡張されてきた。すべてのトランスポーターは、オリジナルのTCDBの設計に従って、進化的、構造的、機能的な観点から分類されている(ref.2,8)。また、その結果として、そのタンパク質は、そのタンパク質が持つ機能や機能に関連していることが明らかになっている。TCデータベース(TCDB)の初期のバージョンは、Nucleic Acids Research (6,9-11)に掲載されている。

 現在、TCDBでは18300以上の研究論文に掲載されたデータにアクセスすることができる。これらの情報は、TCDB内の記述と階層構造に統合されている。このデータベースには、あらゆる種類の生物から得られた15 528のシングルコンポーネントまたはマルチコンポーネントの輸送系が収録されており、そのうちの1567の輸送系については、PDB(ref.12)の高分解能3次元構造のアクセッションが利用可能である。これらのシステムは、系統と機能に基づいて1536のトランスポーターファミリーに分類されている。これらのファミリーの多くは遠縁関係にあることが判明しており、現在では82のスーパーファミリーに分類されている(ref.13)。これは、2016年の報告(ref.9)と比較して50%以上の全体的な増加を表している。TCDBは、このソフトウェアがより洗練され、輸送システムに関する新たな研究が発表されるにつれて、継続的に更新されている。これらの進歩により、国際的な科学コミュニティで使用するためのデータベースの有用性が高まることを期待している。

 トランスポトームの解析には、TCDBで得られた知識に基づいて、シングルコンポーネント輸送系とマルチコンポーネント輸送系を別々に解析する必要がある。さらに、新たに配列決定された(メタ)ゲノムは、既知のトランスポーターとの類似性がほとんどない、あるいは全くない潜在的なトランスポーターを同定する機会を提供する。これらのタンパク質をTCDBに組み込むことで、データベースのカバレッジと配列の多様性を高めていく。これらの課題に取り組むために多くのプログラムを開発してきた(表1および補足ファイル1参照、Githubのリンク)。

 

簡単に機能を見ていく。

http://www.tcdb.orgにアクセスする。

f:id:kazumaxneo:20201122063800p:plain

 

BLAST - blastp, blastxによる類似したトランスポーターの検索

f:id:kazumaxneo:20201122070242p:plain

 

SUBSTRATE EXPLORER - 輸送される分子による検索

f:id:kazumaxneo:20201122070423p:plain

 

 

PSI-BLAST - PSI-BLASTによるトランスポーターの検索。進化的に距離があるトランスポーターもヒットしやすい。

f:id:kazumaxneo:20201122072943p:plain

 

STRUCTURE DATA

f:id:kazumaxneo:20201122083649p:plain

蛋白質構造データバンク(PDB)にリンクしている。

 

HUMAN TRANSPORTERS

f:id:kazumaxneo:20201122083934p:plain

 

TRANSPORTERS & DISEASES

f:id:kazumaxneo:20201122084004p:plain

 

BIO-TOOLS

Welcome to BioTools

f:id:kazumaxneo:20201122083337p:plain

 

 

 


Github



引用

The Transporter Classification Database (TCDB): 2021 update
Milton H Saier, Jr, Vamsee S Reddy, Gabriel Moreno-Hagelsieb, Kevin J Hendargo, Yichi Zhang, Vasu Iddamsetty, Katie Jing Kay Lam, Nuo Tian, Steven Russum, Jianing Wang, Arturo Medrano-Soto
Nucleic Acids Research, Published: 10 November 2020


The Transporter Classification Database (TCDB): recent advances

Milton H Saier Jr , Vamsee S Reddy, Brian V Tsu, Muhammad Saad Ahmed, Chun Li , Gabriel Moreno-Hagelsieb

Nucleic Acids Res. 2016 Jan 4;44(D1):D372-9

 

The transporter classification database

Milton H Saier Jr, Vamsee S Reddy, Dorjee G Tamang, Ake Västermark

Nucleic Acids Res. 2014 Jan;42(Database issue):D251-8


The Transporter Classification Database: recent advances

Milton H Saier Jr, Ming Ren Yen, Keith Noto, Dorjee G Tamang, Charles Elkan

Nucleic Acids Res. 2009 Jan;37(Database issue):D274-8

 

TCDB: the Transporter Classification Database for membrane transport protein analyses and information
Milton H. Saier, Jr,* Can V. Tran, Ravi D. Barabote

Nucleic Acids Res. 2006 Jan 1; 34(Database issue): D181–D186

 

参考

タンパク質の立体構造予測研究の展開と ゲノム情報解析への応用

 

https://www.jstage.jst.go.jp/article/ciqs/toyoha/0/toyoha_0_J06/_pdf/-char/ja

 



 

ロングリードおよび長い配列のアライナー LRA

 

 1分子シークエンシング(SMS)装置からのロングリードや、SMSアセンブリからのメガベーススケールのコンティグをアラインメントしてバリエーションを検出することは、計算量的に困難である。長い配列を効率的にアラインメントするための1つのアプローチは、スパース動的プログラミング(SDP)であり、配列とゲノムの間に完全に一致するものを見つけ、大まかなアラインメントを表す最適な一致の連鎖を見つける。配列のばらつきは、ギャップ長の凸関数であるギャップペナルティを用いてアラインメントをスコア化すると、より正確にモデル化される。これまでのSDPの実装では、ばらつきを正確にモデル化できない線形コストのギャップ関数を使用していたため、凸のギャップペナルティを持つアライメントの実装は効率が悪いか、ヒューリスティックを使用していた。lraを使用して、PacBioやOxford Nanopore (ONT)の装置やde novoアセンブリコンティグからの長い配列のアラインメントを行った。すべてのデータタイプにおいて、lraの実行時間は、SAMアライメントを生成する際の最新のアライナ-minimap2の52-168%、代替手法であるngmlrの9-15%であった。このアラインメントアプローチは、PacBioデータセットにおけるSVコールの追加証拠を提供するために使用される可能性があり、現在のSV検出アルゴリズムではONTデータに対する感度と特異性が向上する。lraアライメントを用いたpbsvを用いて発見されたコールの数は、同じデータ上のminimap2アライメントを用いたコールの98.3-98.6%以内であり、Truvari解析によるF1スコアの名目上の0.2-0.4%の増加を与える。Snifflesを使用して呼び出されたSVを持つONTデータでは、lraアライメントからのコール数はminimap2ベースのコール数よりも3%大きく、ngmlrベースのコール数よりも30%大きく、TruvariのF1スコアは4.6~5.5%増加する。de novoアセンブリコンティグからのバリエーションコールに適用すると、minimap2+paftoolsと比較してSVコールが5.8%増加し、Truvari F1スコアが4.3%増加した。



インストール

python3.7環境でcondaを使って導入した(ubuntu18.04)。

依存

Github

#bioconda (link)
conda install -c bioconda lra -y

 

実行方法

1, indexing (global and local)

シーケンシングマシンによってインデックスのパラメータ設定は異なる

#ONT
lra index -ONT ref.fa

#pacbio CCS
lra index -CCS ref.fa

#pacbio CLR
lra index -CLR ref.fa

 

2, mapping

#ONT
lra align -ONT ref.fa read.fa -t 16 -p s > output.sam

#pacbio CCS
lra align -CCS ref.fa read.fa -t 16 -p s > output.sam

#pacbio CLR
lra align -CLR ref.fa read.fa -t 16 -p s > output.sam

 

PAF、SAM、BED形式で出力可能 

#SAM
lra align -ONT ref.fa read.fa -t 16 -p s > output.sam

#PAF
lra align -ONT ref.fa read.fa -t 16 -p p > output.paf

#BED
lra align -ONT ref.fa read.fa -t 16 -p b > output.bed

 

引用

lra: the Long Read Aligner for Sequences and Contigs

Jingwen Ren, Mark Chaisson

bioRxiv, Posted November 17, 2020

 

関連


 

 

ショートリードのアセンブラ Clover

 

 次世代シーケンシング技術は、低コストでハイスループットのリードを生産することでゲノミクスに革命をもたらし、この進歩に伴いde novoアセンブラの開発が促された。de Bruijnグラフに基づく複数のアセンブラ法は、Illuminaリードに対して効率的であることが示されてきた。しかし、シーケンサーで発生するシークエンシングエラーは、de novoアセンブリーの解析を複雑にし、下流のゲノム研究の質に影響を与える。
 本論文では、Illuminaプラットフォームで発生するシーケンシングエラーに対応するために、オーバーラップ・レイアウト・コンセンサスの概念に基づく新しいk-merクラスタリングアプローチを利用したClover(クラスタリング指向型de novoアセンブラ)を開発した。さらに、3つのデータセット(Staphylococcus aureus, Rhodobacter sphaeroides, ヒト第14番染色体)について、いくつかのde Bruijnグラフアセンブラ(ABySS, SOAPdenovo, SPAdes, Velvet)、オーバーラップレイアウトコンセンサスアセンブラ(Bambus2, CABOG, MSR-CA)、ストリンググラフアセンブラ(SGA)を用いて、Cloverの性能を評価した。その結果、Cloverは、SOAPdenovoを除き、実行時間において有意な競争力を維持しながら、corrected N50とEサイズの点で優れたアセンブリ品質を達成していることを示した。また、Cloverは、Acinetobacter baumannii TYTH-1とMorganella morganii KT細菌ゲノムの配列決定プロジェクトにも寄与した。
 オーバーラップ・レイアウト・コンセンサス法の柔軟性とde Bruijnグラフ法の効率性を統合したCloverのマーベル・クラスタリングベースのアプローチは、デノボ・アセンブリにおいて高い可能性を秘めている。Cloverはオープンソース・ソフトウェアとして https://oz.nthu.edu.tw/~d9562563/src.html から自由に利用できる。

 

HP

https://oz.nthu.edu.tw/~d9562563/src.html

A simple example

https://oz.nthu.edu.tw/~d9562563/html/testcase.html

 

インストール

ubuntu18.04でテストした。

依存

condaでpython2.7の仮想環境を作ってテストした(macos10.14使用)。

HPからダウンロードする。

tar -zxvf clover-2.0.tar.gz
cd clover-2.0
chmod u+x clover

./clover

$ ./clover 

Please set the input file.

 

Clover is the command line tool. The user can run Clover straightforward with 

following parameters:

 

  $ clover -k <Length of k-mer> [options] -i1 <Input file1> [-i2 <Input file2>]

 

 

If only a read file is used without mate pair, the parameter -i2 can be omitted.

For example, if one library of single read file frag.fastq is used:

  $ clover -k 40 -i1 frag.fastq

 

If paired read files are used, file name of -i2 must correspond to -i1.

 

For example, if one library of paired read files frag1.fq and frag2.fq is used:

  $ clover -k 40 -i1 frag1.fq -i2 frag2.fq

 

If two libraries of paired read files frag1.fq, frag2.fq, short1.fq and short2.fq 

are used, where assume that frag1.fq corresponds to frag2.fq and short1.fq 

corresponds to short2.fq:

  $ clover -k 40 -i1 frag1.fq,short1.fq -i2 frag2.fq,short2.fq

 

The file formats accepted by Clover are 'fasta' and 'fastq', which can be 

distinguished by their filename extensions (.fa, .fasta, .fq, .fastq, .fatq).

 

For more information, please type:

 

 

テストラン

GAGEのデータセットを使ったテストランの流れが記載されている(link)。

wget https://oz.nthu.edu.tw/~d9562563/data/testdata.tar.gz
tar -zxvf testdata.tar.gz
cd testdata/

../clover -k 46 -p 0 -i1 frag_1.fastq,shortjump_1.fastq -i2 frag_2.fastq,shortjump_2.fastq -cs 7 -ss 3 -is 180,3500 -hp 0.6 -pm -ml 200

出力

f:id:kazumaxneo:20201120132734p:plain



 

引用

Clover: a clustering-oriented de novo assembler for Illumina sequences

Ming-Feng Hsieh, Chin Lung Lu & Chuan Yi Tang

BMC Bioinformatics volume 21, Article number: 528 (2020)

 

関連


 

*1

libpython2.7.so.1.0が見えないと言われたので、find / -name libpython2.7.so.1.0して、見つかったパスから/libにシンボリックリンクを張って対処した。

ln -s <path>/to/libpython2.7.so.1.0 /lib/libpython2.7.so.1.0

まだエラーが起こる。

勘違いしていた。問題なし。

 

 

NCBIのGenBankゲノムアセンブリ (GCA) とRefSeqゲノムアセンブリ(GCF)

2020/11/19 誤字修正

 

識別子がGCA_で始まるゲノムアセンブリは、GenBankアセンブリと呼ばれる。GenBankアセンブリは、ユーザーがサブミットしたゲノムアセンブリを意味する。一方、識別子がGCF_で始まるゲノムアセンブリは、RefSeqのアセンブリになる。こちらも元はユーザーがサブミットしたゲノムアセンブリに由来しているが、NCBI側で短すぎる配列をのぞいたり、オルガネラゲノム配列を追加するなどのキュレーションが行われたゲノムアセンブリになる。

 

詳しくはこちらを読んでください。

What is the difference between a GenBank (GCA) and RefSeq (GCF) genome assembly?

 

スモールゲノムを例に違いを見てみる。NCBIbrowse by Organismからゲノムアセンブリを検索する。

f:id:kazumaxneo:20201118223248p:plain



Helicobacter pyloriの完全長ゲノムアセンブリを検索してみた。

f:id:kazumaxneo:20201118224055p:plain

 

検索結果。ここで表示されているのはGCA_のアセンブリIDになる。一番上のGCA_000008525.1をクリックした。

f:id:kazumaxneo:20201118224126p:plain

ここでは関係ないが、Helicobacter pyloriのゲノムサイズは1.7Mb程度と非常に小さいのが目を引く。

 

クリック結果。BiosampleやBioprojectのアクセッションIDが目立つが、画面右端にRefSeq assemblyとGenBank assemblyのNCBI FTPサーバーリンクがある。ここからもゲノムアセンブリアノテーションのフィーチャーファイル、サマリーファイル等にアクセスできる。それぞれ開いてみる。

f:id:kazumaxneo:20201118230033p:plain

 

RefSeq assembly。/genome/all/GCF/内にあり、ファイル名がGCF_で始まっている。

f:id:kazumaxneo:20201118230137p:plain

GenBank assembly。/genome/all/GCA/内にあり、ファイル名がGCA_で始まっている。

f:id:kazumaxneo:20201118230114p:plain

 

ページの下からもRefSeqアセンブリGenBankアセンブリにリンクしている。

f:id:kazumaxneo:20201118231707p:plain

ゲノムによってはどちらのリンクもGCF_で始まるものがあったりする。IDをよく見ること。

 

では実際に何か違いがあるのだろうか。

GCFとGCAのゲノムアセンブリのサイズを調べた。

seqkit stats -a *fna

f:id:kazumaxneo:20201118232422p:plain

 ゲノムサイズは同じだった。このゲノムは完全長であり、プラスミドもないため、修正の必要性がない。

 

次に、GCFとGCAアノテーションからprotein.faaのサイズを調べた。

seqkit stats *faa

f:id:kazumaxneo:20201118232626p:plain

アノテーションされているプロテイン数には38個差があった。アノテーションのプログラムが異なるので、おそらく数以外にも細かな違いが生じていると予測される。

 

まとめ
たとえスモールゲノムであっても、GenBankゲノムアセンブリとRefSeqゲノムアセンブリやそのアノテーションファイルには結構な違いがある場合がある。論文内でlocus IDなどが書いてあって、その情報を頼りに遺伝子やタンパク質配列情報を取ってきて調べる際には特に注意が必要になる。

 

 

参考

井上先生の解説も参考にしました。


関連