macでインフォマティクス

macでインフォマティクス

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

mito-finder

 

 核内ゲノムをアセンブルしている間に、Falconやcanu、その他の素晴らしいアセンブラミトコンドリアゲノムのアセンブルに完全に失敗しているか、ミトコンドリアゲノムのアセンブルを誤っていることに気がついた。ミトコンドリアゲノムが入力DNAの大部分よりも小さいことを考えれば、これは驚くべきことではない。そこで、アセンブラに頼るのではなく、別のルートでミトコンドリアゲノムをアセンブルすることにした。

 

インストール

イデアだけでコードは登録されていない。

GIthub

必要なツールをインストールする。

conda create -n mitofinder-env -y
conda activate mitofinder-env
conda install -c bioconda -y blasr mummer samtools

 

実行手順

1、公開されているミトコンドリアゲノムの中でもっとも近いミトコンドリアゲノムにマッピングする。blasrを使う。

blasr raw_pacbio.fasta mito.fasta --bestn 1 -m 1 --nproc 12 > mito.m1

 

2、サイズセレクトする。18-kb以上であることが分かっているなら、それ以上のリードを選ぶ。そのため、まずawkを使い18-kb以上のアラインメントのリード名をチェックし、

awk '{if($8-$7>18000) print $1"\t"$8-$7"\t"$12}' mito.m1 > output

その中から、もっとも長いアラインメントのリードを抽出する。 IDがxxxxxなら

samtools faidx raw_pacbio.fasta
samtools faidx raw_pacbio.fasta xxxxx > longest.fa

 

3、環状DNAのオーバーラップがあるのかnucmerを使って確認する。

 nucmer -maxmatch --nosimplify mito_read.fasta mito_read.fasta

out.deltaを開き、オーバーラップ領域があれば、そのオーバーラップ領域をトリムする。

 

4、最後にquiverで研磨する。

 

おそらくミトコンドリアサイズが短い動物や真菌などのミトコンドリアゲノム向けのプロトコルです。ミトコンドリアサイズが大きな植物には使えないと思います。注意して下さい。

引用

Hidden genetic variation shapes the structure of functional elements in Drosophila

Mahul Chakraborty, Nicholas W. VanKuren, Roy Zhao, Xinwen Zhang, Shannon Kalsow & J. J. Emerson
Nature Genetics volume 50, pages20–25(2018)

 

 

 

ゲノム間のシンテニー領域を調べる MCScanX

2020 6/21 タイトル変更

2021 11/23 インストール手順修正

 

 MCScanは、複数のゲノムまたはサブゲノムをスキャンして、相同性のあると思われる染色体領域を特定し、遺伝子をアンカーにしてこれらの領域を整列させることができるアルゴリズムである。MCScanXツールキットは、syntenyとcollinearityの検出のために調整されたMCScanアルゴリズムを実装しており、結果の可視化と下流の追加解析のための14のユーティリティプログラムを組み込むことで、オリジナルのソフトウェアを拡張している。いくつかの植物ゲノムと遺伝子ファミリーへのMCScanXの適用例を示す。染色体の構造変化を効果的に解析し、系統や分類群の適応に寄与すると考えられる遺伝子ファミリーの拡大の歴史を明らかにするためにMCScanXを使用することができる。gene duplicationの様々なモードの統合的なビューは、特定のファミリーにおける従来の遺伝子ツリー解析を補完できる。MCScanXのソースコードとドキュメントは http://chibba.pgml.uga.edu/mcscan2/ から自由に入手できる。(2022/11/16現在リンク切れ)

 

 MCScanXアルゴリズムはMCScan(19)の修正版である。全ゲノムBLASTPの結果は、染色体とscaffoldsのすべての可能なペアのコリニアブロックを計算するために使用される。まず、BLASTPのマッチは遺伝子の位置に応じてソートされる。タンデム配列に起因する局所的なコリニアな遺伝子ペアの多さを避けるために、連続したBLASTPマッチが共通の遺伝子を持ち、そのペアの遺伝子が5つ以下の遺伝子で区切られている場合、これらのマッチは、BLASTPのE値が最も小さい代表ペアを使用して折りたたまれる。次に、動的プログラミングは、2つの遺伝子ペア、uとvが、uがvに先行するパス上にあると仮定して、以下のスコアリングスキーマを使用して、最も高いスコアリングパス(すなわち、コリネア遺伝子ペアの連鎖)を見つけるために採用される。(以下略)

 

MCScanX HP 

http://chibba.pgml.uga.edu/mcscan2/

Twitter

https://twitter.com/search?q=MCScanX&src=typed_query

 

インストール

上のHPからダウンロードした。

GIthub

git clone https://github.com/wyp1125/MCScanX.git
cd MCScanX/
make -j 8
export PATH=$PATH:$PWD

./MCScanX

$ ./MCScanX

[Usage] ./MCScanX prefix_fn [options]

 -k  MATCH_SCORE, final score=MATCH_SCORE+NUM_GAPS*GAP_PENALTY

     (default: 50)

 -g  GAP_PENALTY, gap penalty (default: -1)

 -s  MATCH_SIZE, number of genes required to call a collinear block

     (default: 5)

 -e  E_VALUE, alignment significance (default: 1e-05)

 -m  MAX_GAPS, maximum gaps allowed (default: 25)

 -w  OVERLAP_WINDOW, maximum distance (# of genes) to collapse BLAST matches (default: 5)

 -a  only builds the pairwise blocks (.collinearity file)

 -b  patterns of collinear blocks. 0:intra- and inter-species (default); 1:intra-species; 2:inter-species

 -h  print this help page

 

./MCScanX_h

$ ./MCScanX_h 

[Usage] ./MCScanX_h prefix_fn [options]

 -k  MATCH_SCORE, final score=MATCH_SCORE+NUM_GAPS*GAP_PENALTY

     (default: 50)

 -g  GAP_PENALTY, gap penalty (default: -1)

 -s  MATCH_SIZE, number of genes required to call collinear blocks

     (default: 5)

 -e  E_VALUE, alignment significance (default: 1e-05)

 -m  MAX_GAPS, maximum gaps allowed (default: 25)

 -w  OVERLAP_WINDOW, maximum distance (# of genes) to collapse BLAST matches (default: 5)

 -a  only builds the pairwise blocks (.collinearity file)

 -b  patterns of collinear blocks. 0:intra- and inter-species (default); 1:intra-species; 2:inter-species

 -c  whether to consider homology scores. 0:not consider (default); 1: lower preferred; 2: higher preferred

 -h  print this help page

>./Duplicate_gene_classifier

$ ./Duplicate_gene_classifier

[Usage] ./Duplicate_gene_classifier prefix_fn [options]

 -k  MATCH_SCORE, final score=MATCH_SCORE+NUM_GAPS*GAP_PENALTY

     (default: 50)

 -g  GAP_PENALTY, gap penalty (default: -1)

 -s  MATCH_SIZE, number of genes required to call a collinear block

     (default: 5)

 -e  E_VALUE, alignment significance (default: 1e-05)

 -w  OVERLAP_WINDOW, maximum distance (# of genes) to collapse BLAST matches (default: 5)

 -m  MAX_GAPS, maximum gaps allowed (default: 25)

 -n  N_PROXIMAL, maximum distance (# of genes) to call proximal (default: 10)

 -h  print this help page

 

 

実行方法

調べたいゲノムのprotein.fasta、リファレンス側のゲノムのprotein.fastaとBED/GFF ファイルを用意する。BED/GFFファイルには染色体名、遺伝子名、スタート、エンドのポジションが書かれていなければならない。

$ head data/os_sb.gff 

os5 Os05t0507200 25150045 25150392

os2 Os02t0774100 33540257 33541838

os5 Os05t0480600 23720056 23723583

os6 Os06t0299300 11200887 11202509

os7 Os07t0587100 24548577 24551667

os5 Os05t0526700 26286838 26287599

os1 Os01t0974500 44835553 44838812

os6 Os06t0608600 25100853 25103401

os6 Os06t0242600 7400125 7400986

os11 Os11t0120600 904677 906689

ファイル名の拡張子以外の名前は統一しておく必要がある。 BED ファイル名がxyz.gffなら、blast出力はxyz.blastという風になっていないといけない。

 

1、blastpの実行

blastall -i query_file -d database -p blastp -e 1e-10 -b 5 -v 5 -m 8 -o xyz.blast

 

2、MCScanXの実行。テストランならdataディレクトリのos_sb.blastとos_sb.gffを対象に以下のようにランする。

 

./MCScanX data/os_sb

出力

f:id:kazumaxneo:20200620142515p:plain

> head -n 20 data/os_sb.collinearity

$ head -n 20 data/os_sb.collinearity 

############### Parameters ###############

# MATCH_SCORE: 50

# MATCH_SIZE: 5

# GAP_PENALTY: -1

# OVERLAP_WINDOW: 5

# E_VALUE: 1e-05

# MAX GAPS: 25

############### Statistics ###############

# Number of collinear genes: 33937, Percentage: 50.45

# Number of all genes: 67274

##########################################

## Alignment 0: score=559.0 e_value=1.8e-34 N=13 os1&os1 plus

  0-  0: Os01t0584200 Os01t0713850   7e-15

  0-  1: Os01t0584900 Os01t0714800   2e-30

  0-  2: Os01t0588200 Os01t0715500   9e-12

  0-  3: Os01t0589200 Os01t0716500   7e-62

  0-  4: Os01t0593200 Os01t0719000   7e-25

  0-  5: Os01t0593700 Os01t0719300       0

  0-  6: Os01t0595100 Os01t0719400   2e-09

  0-  7: Os01t0595201 Os01t0719600   1e-22

blast出力以外にペアワイズのアラインメントがあればそれを使用してランすることもできる(MCScanX_hを使用)。ファイル例はat_gm_pt_vv.homologyを参照。

 

MCScanXによるコリニアブロックのマルチプルアラインメントを表示したHTML出力。

f:id:kazumaxneo:20200620145939p:plain

sb9.html

f:id:kazumaxneo:20200620150030p:plain


 

 

もう1つのテストデータ

./MCScanX data/at_vv

 

 

他に重複遺伝子を見つけるためのduplicate_gene_classifierがある。

./duplicate_gene_classifier data/os_sb

 

追記

結果はSynVisioを使って視覚化できる。

SynVisio : Synteny Browser

f:id:kazumaxneo:20200620145835p:plain

MCScanXのテストデータを視覚化した。

 

引用
MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity
Yupeng Wang, Haibao Tang, Jeremy D. DeBarry, Xu Tan, Jingping Li, Xiyin Wang, Tae-ho Lee, Huizhe Jin, Barry Marler, Hui Guo, Jessica C. Kissinger, Andrew H. Paterson

Nucleic Acids Res. 2012 Apr; 40(7), Published online 2012 Jan 4

 

rRNA配列を探索し、分類、ツリー表示するSILVAのACTサービス

2020 6/19  説明追加

 

相同配列の解析では、マルチプルシーケンスアラインメント(MSA)の計算がボトルネックになっている。特にリボソームRNA(rRNA)のようなマーカー遺伝子の場合、数百万の配列がすでに公開されており、個々の研究で数十万の新しい配列を簡単に作成することができる。このような数に対応するための方法が開発されてきたが、精度の要求を満たすためにはさらなる改良が必要である。

  本研究では、SILVAリボソームRNAプロジェクトが提供するrRNA遺伝子データベースのアラインメントに使用したSILVA Incremental Aligner (SINA)を紹介する。SINAは、k-mer検索とpartial order alignment (POA) の組み合わせを用いて、非常に高いアライメント精度を維持しつつ、高いスループット性能要求を満足させる。SINAは、一般的に使用されている高スループットMSAプログラムPyNASTおよびmothurと比較して評価された。3つのBRAliBase IIIベンチマークMSAは99.3、97.6、96.1の精度で再現できた。38,772の配列からなるより大きなベンチマークMSAは、1,000および5,000の配列からなるリファレンスMSAを用いても98.9および99.3%の精度で再現できた。SINAは、実行されたすべてのベンチマークにおいて、PyNASTおよびmothurよりも高い精度を達成することができた。

 

 

webサービス

ここではSINAを使ったSILVAのrRNA探索と分類、およびツリーへの視覚化のサービス、ACTを簡単に紹介します。

https://www.arb-silva.de/aligner/ にアクセスする。

f:id:kazumaxneo:20200618234755p:plain

 

データのアップロード

1、配列をアップロードする。既にrRNAと分かっている配列だけ入力しても良いし、de novo transcriptomeのアセンブリ配列全体がそれほど大きくなければ、それを登録してSSU / LSU rRNAを探索することも可能(サイズ制限あり)。

f:id:kazumaxneo:20200618234858p:plain

 

2、SSUかLSUかを選択する。SSUの方がデータベースは充実している。

f:id:kazumaxneo:20200618235632p:plain

SILVAが全てのrRNA配列を管理しているわけではないです。真核生物のrRNA配列などを探索する場合は特に注意して下さい、 

 

3、Search and classifyにチェックを付け、Min identityを指定する。デフォルトでは0.95だが、データベースから遠い配列だと0.95では見つからない。

f:id:kazumaxneo:20200618235706p:plain

 

4、Compute treeにチェックをつける。Denovo including neighborを選ぶと、近いrRNA配列を含めて系統推定してくれる。

f:id:kazumaxneo:20200618235932p:plain

 

5、デフォルトでは近似のFastTree。非常に早く終わる。

f:id:kazumaxneo:20200619000645p:plain

 

6、SSUだと右下からDomainを選択できる。間違えないようにする。

f:id:kazumaxneo:20200619001526p:plain

LSUだと選べるデータベースが限定される。release138で選べるGTDBも選択できない。domainも選択できなくなる。

f:id:kazumaxneo:20200619001032p:plain

 

7、他のパラメータを決めてRn toolをクリック。サブミットに成功すると、下のジョブ管理ウィンドウに進捗が表示される。

f:id:kazumaxneo:20200619000934p:plain

 

 

出力内容

ジョブが終わると多重整列結果、系統推定結果のnewickファイル等をダウンロードしたり、その場で視覚化できるようになる。

f:id:kazumaxneo:20200619082704p:plain


View in wasabi => Aligned sequence with neighbors 

f:id:kazumaxneo:20200619083332p:plain

 

Denovo including neighborを選んでいれば、SILVAの近いrRNAも含めて多重整列、系統推定が行われる。

f:id:kazumaxneo:20200619094725p:plain

取得した配列がrRNA配列なのか、部分的な相同性しか示さない異なる配列なのか視覚化すれば明確にわかる。

 

SILBAの配列も含めた多重整列結果やnewickファイルはダウンロードできる。近似の最尤法などに頼りたくないならローカルマシンで計算を行えばよい。

引用

SINA: Accurate high-throughput multiple sequence alignment of ribosomal RNA genes
Elmar Pruesse, Jörg Peplies, Frank Oliver Glöckner
Bioinformatics, Volume 28, Issue 14, 15 July 2012, Pages 1823–1829

繰り返し配列を分析する RepeatProfiler

2020 7/7 リンクミス修正

2020 12/6  論文追記

 

 モデル生物におけるDNAリピートの研究は、ゲノムの進化や表現型の変化を促進する多くのプロセスにおけるリピートDNAの役割を浮き彫りにしている。反復配列はシングルコピーDNAよりもはるかにダイナミックであるため、進化の遅いゲノム領域の配列では明らかにならないような、短い時間スケールでの進化の歴史のシグナルを明らかにすることができる。リピートを研究するための多くのツールは、ゲノムアセンブリやリピートライブラリなど、既存のゲノムリソースを持つ生物を対象としている。しかし、ゲノム資源が限られている多様な非モデル生物群の進化の歴史を解明する上では、リピート変異のシグナルが特に有用であることが証明されるかもしれない。ここでは、カバー率の低いショートリードシークエンシングデータからリピートDNAプロファイルを生成、可視化、比較するためのツールであるRepeatProfilerを紹介する。RepeatProfilerは、リピートDNAカバレッジデプスプロファイルの生成と可視化を自動化し、サンプル間でプロファイル形状を統計的に比較することができる。さらに、RepeatProfilerは、プロファイル間の配列変異からシグナルを抽出し、系統解析を用いて分子形態学的特性として解析することで、プロファイルの比較を容易にする。著者らはカブトムシ(Bembidion)、ハエ(Drosophila)、トマト(Solanum)のデータセットでRepeatProfilerを検証した。種の区別、比較ゲノム、リピート生物学の研究のための高分解能データソースとしてのリピートDNAプロファイルの可能性を強調する。

 

https://twitter.com/search?q=RepeatProfiler&src=typed_query

 

インストール

brewでも導入できるが、checksumエラーが出たのでdocker imageを使ってテストした。

Github

#homebrew
brew install HounerX/homebrew-repo/repeatprof

#dockerhub
docker pull durberg7/repeatprof

 > docker run --rm -it durberg7/repeatprof bash repeatprof

$ sudo docker run --rm -it durberg7/repeatprof bash repeatprof

RepeatProfiler v 0.92 -prerelease- x Singlecopy x new x Correlation x

please use -h flag to view help menu for how to use this tool

____ ____ ___  ____ ____ ___    ___  ____ ____ ____ _ _    ____ ____

|__/ |___ |__] |___ |__|  |     |__] |__/ |  | |___ | |    |___ |__/ 

|  \ |___ |    |___ |  |  |     |    |  \ |__| |    | |___ |___ |  \ 

 

 

sudo docker run --rm -it durberg7/repeatprof bash repeatprof -h

$ sudo docker run --rm -it durberg7/repeatprof bash repeatprof -h

generating symlinks

 

Usage:

repeatprof profile <-p for paired reads or  -u for unpaired> <the refrence sequence path > <path of the folder containing reads> [opitonal flags]

 

optional flags:

-o <folder_path> use to specifiy an directory where the final output file  be directed to. Default: current directory

-corr type this flag to make the correlation analysis. user_groups.txt will be assumed to be in same directory

-corr <user_groups.txt path >   type this flag to make the correlation analysis providing User_provided.txt path

-k use this flag if you want to keep the sorted bam files of the alignments in the final output folder

 

other optional flags can be used:

--very-sensitive bowtie alignment setting. Default:--very-sensitive

--very-fast bowtie alignment setting. Default:--very-sensitive

--local bowtie alignment setting. Default:--very-sensitive

--fast bowtie alignment setting. Default:--very-sensitive

--very-fast-local bowtie alignment setting. Default:--very-sensitive

--fast-local bowtie alignment setting. Default:--very-sensitive

--sensitive bowtie alignment setting. Default:--very-sensitive

--very-sensitive bowtie alignment setting. Default:--very-sensitive

--very-sensitive-local         bowtie alignment setting. Default:--very-sensitive

--sensitive-local bowtie alignment setting. Default:--very-sensitive

 

Supported input formats:

reference: .fa .fasta .txt 

Paired reads: _R1.fastq _R1.fastq.gz _R1.fq  _R1.fq.gz  _1.fastq _1.fastq.gz _1.fq  _1.fq.gz  

              _R2.fastq _R2.fastq.gz _R2.fq  _R2.fq.gz  _2.fastq _2.fastq.gz _2.fq  _2.fq.gz  

Unpaired reads:  .fastq .fastq.gz .fq  .fq.gz  .fastq .fastq.gz .fq

 

references need to be in fasta format. They can be multi-sequence or single sequence fasta file.

 

--------------------------------other usages: 

 

repeatprof pre-corr <-u for unpaired reads  or -p paired reads> <path of the folder containing reads> use this command to help prepare user_groups.txt for correlation graph produced by -corr flag in profile 

 

repeatprof pre-corr -v  use this command to review your edited user_groups.txt with your desired groups. This ensures that it is still in the correct file format. Need to be in same directory as file

 

repeatprof clean use this command to cleanup a broken/terminated run

 

check our github page https://github.com/johnssproul/RepeatProfiler for more  detailed info  on how to run the tool/sample inputs/ output explanation/tutorial 

(repeatmodeler) kazu@kazu:~/Downloads/test/sample_data$ 

 

 

 

テストラン

ランにはシークエンシングデータのfastqとリファレンスのFASTAファイルが必要。テストデータを使って動作確認する。

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

wget https://github.com/johnssproul/RepeatProfiler/releases/download/0.96/sample_data.zip
unzip sample_data.zip && cd sample_data/

sample_dataの中身

f:id:kazumaxneo:20200618002112p:plain

(Valid extensions are '.fastq', '.fq'. Compressed reads ‘i.e., .fastq.gz’ also supported. For paired reads, the last string before the file extension should be ‘_1’ for Read1 and ‘_2’ for Read2 (alternatively ‘_R1’ and ‘_R2’ may be used).)

 

 2、dockerイメージを使ったRepeatProfilerのprofileコマンドの実行

リファレンスFASTAとfastqが置いてあるディレクトリを指定する。

sudo docker run --rm -itv $PWD:/data durberg7/repeatprof bash repeatprof profile -p /data/reference.fa /data/ -o /data

 出力

f:id:kazumaxneo:20200617235308p:plain

variant_profiles.pdf

f:id:kazumaxneo:20200618004239p:plain

拡大

f:id:kazumaxneo:20200618004427p:plain

このファイルには、このリファレンス内のすべてのサンプルの分散強調プロファイルが含まれている。これらのプロファイルは、サンプル間でリファレンス配列に対するバリアントのbp分解能を示しており、サンプル間で比較した場合にも興味深いパターンを明らかにすることができる。

 

scaled_profiles_allrefs/melR1.pdf

f:id:kazumaxneo:20200618004810p:plain

scaled_profiles_allrefsディレクトリには、同じカラースケールで表示されるラン全体(すなわち全リファレンス配列)のカラー強調プロファイルが出力される。これは,各リファレンスの出力フォルダにある 'scaled_profiles.pdf' ファイルと似ているが,カラースケールがリファレンスシーケンス内で観測された最大値ではなく,実行中のすべてのプロファイルの最大カバレッジに基づいて設定されている点が異なる。繰り返しにまたがる興味深いパターンを識別するのに役立つ。

 

実行方法

同じリファレンスの複数サンプルfastqが存在する場合(複数individualsのデータなど)、そのfastqディレクトリを指定してランすることでサンプル間のプロファイル形状を比較できる。その場合、user_groups.txtを指定してfastqのグループを指定する必要がある。このuser_groups.txtはrepeatprof pre-corrコマンドで自動作成することができる。

repeatprof pre-corr <'-u' for unpaired reads or '-p' for paired reads> <path reads folder>

 詳細はGithub参照。

 

引用

RepeatProfiler: a pipeline for visualization and comparative analysis of repetitive DNA profiles

S. Negm, A. Greenberg, A.M. Larracuente, J.S. Sproul

bioRxiv, Posted May 26, 2020

 

2020 12/16

RepeatProfiler: a pipeline for visualization and comparative analysis of repetitive DNA profiles
S. Negm A. Greenberg A.M. Larracuente J.S. Sproul
https://doi.org/10.1111/1755-0998.13305

 

TAMA

 

微生物は様々な環境の中で重要な役割を果たしている。微生物の組成を特定し、その存在量を推定することで、環境試料中の微生物の相互作用を理解することができる。微生物の環境をより深く理解するために、微生物ゲノムのメタゲノムアセンブリを用いて、環境試料中の微生物の組成を調べることが行われている。この分野では、様々なアルゴリズムに基づいた分類学解析ツールが開発されているが、同じメタゲノムデータセットからの解析結果のばらつきが、多くの研究者にとって大きな障害となっている。

 本研究では、3つの異なる分類学解析ツールの出力をインテリジェントに統合することで、メタゲノム分類学解析のための新しいメタ解析ツールTAMAを提案する。TAMAは、統合されたリファレンスデータベースを用いて、異なる分類ツールからの分類スコアを統合し、メタスコアに基づいて入力されたメタゲノムリードに対して分類を行うもので、既存の分類ツールよりも優れている。TAMAは、様々なベンチマークデータセットを用いて評価した結果、既存のツールよりも優れた性能を示した。また、チーズメタゲノムとヒト腸内メタゲノムの2種類のメタゲノムにおける微生物の相対的な種数分布と組成の違いを得るために適用することにも成功した。

 TAMA は簡単にインストールして利用することができる。TAMA は簡単に導入でき、複数の数や種類のメタゲノムリードサンプルからのメタゲノムリードの分類や相対的な種の豊富さの予測に利用することができる。特に、単一の分類学解析ツールでは信頼性が低い場合には、様々な環境から収集したメタゲノムサンプル中の微生物の組成をより正確に明らかにするためにTAMAを使用することができる。TAMAオープンソースのツールで、https://github.com/jkimlab/TAMA からダウンロードできる。

 

 

インストール

docker imageをビルドしてテストした。

依存

  • perl (v5.22.1)
  • python (2.7.12)
  • java (1.8.0)
  • git (2.7.4)
  • gcc, g++ (5.4.0)
  • make (GNU Make 4.1)
  • zip (3.0)
  • curl (7.47.0)
  • Sort::Key::Natural (perl library)

必要な計算リソース

Disk: (approximately) 300GB
* CLARK: 88GB, Kraken: 188GB, Centrifuge: 9.6GB
- Memory: (approximately) 185GB
* This large memory is required for running the taxonomy analysis tools,
especially CLARK and Kraken

本体 Github

git clone https://github.com/jkimlab/TAMA.git
cd TAMA
docker build -t tama .
docker run -it tama /bin/bash
cd TAMA/

perl TAMA.pl -h

# perl TAMA.pl -h

[Error] There is no proper parameter file.

Usage: ./TAMA.pl [option] --param param.txt 

Options:

-p The number of threads  (default: 1)

-o Path of output directory  (default: Current directory)

-t Save temporary files (default: False)

If you want to save, type 'True'

-h|help Print help page.

Input:

--param (Required) Path of paramter file

# Check the required perl libraries

./setup.pl --check

** Check the requirements..

ok       3.67 Cwd

ok       2.85 File::Basename

ok       1.51 FindBin

ok       2.49 Getopt::Long

ok       0.04 Sort::Key::Natural

 

# Install TAMA package

./setup.pl --install

 

テストラン

テストデータベースをダウンロードしてテストランを実行する。

./setup.pl --example
source ./src/env.sh
bash Example_run.cmd.sh

出力

f:id:kazumaxneo:20200617084734p:plain

f:id:kazumaxneo:20200617084740p:plain

 

データベースの準備

使用時に追記します。

  

引用

TAMA: Improved Metagenomic Sequence Classification Through Meta-Analysis

Mikang Sim, Jongin Lee, Daehwan Lee, Daehong Kwon , Jaebum Kim

BMC Bioinformatics. 2020 May 12;21(1):185

 

ゲノムスケッチを用いて迅速にコホートサンプルの関連性を推定する somalier

 

 複数の空間的または縦断的生検から得られたシーケンシングデータを解釈する際には、サンプルのmix upを検出することが不可欠であるが、生殖細胞変異の研究よりも困難である。腫瘍のほとんどのゲノム研究では、遺伝的変異は腫瘍とサンプル提供者の正常組織のペアワイズ比較によって検出されることが多く、多くの場合、体細胞変異のみが報告されている。その結果、遺伝子型情報がバラバラになるため、生殖細胞変異の遺伝子型のみに基づいてサンプル交換を検出する既存のツールの使用を妨げている。この問題を解決するために、著者らはアラインメント上で直接操作できるsomalierを開発した。Somalierは各サンプルについて情報量の多い遺伝的変異の小さなスケッチを抽出する。その後、何百もの生検や正常組織からのスケッチを1秒以内に比較することができる。この速度は、生殖細胞サンプルの大規模なコホートにおける関連性のチェックにも役立つ。Somalierは、テキスト出力とインタラクティブなビジュアルレポートの両方を生成し、複数の関連性メトリクスを使用してサンプルのスワップの検出と修正を容易にする。このツールを紹介し、正常、腫瘍、無細胞のDNAサンプルを含む5つの神経膠腫サンプルのコホートでその有用性を実証した。また、1000 Genomes Projectの高カバー率シーケンシングデータにSomalierを適用することで、いくつかの関連サンプルを同定することができた。また、そのデータは、ゲノムの構築や多様なシーケンシングデータに適用することができ、学術利用のためにgithub.com/brentp/somalierで自由に利用可能である。

 

 

somalierは既知の多型サイトのリストを取る。数百(または数十)のサイトでさえ、関連性の非常に良い指標になる。最良のサイトは、2つのサンプルが異なる確率を最大にするため、集団の対立遺伝子頻度が0.5に近いサイトである。そのようなサイトのリストは、GRCh37とhg38のリリースにある。これらの部位での遺伝子型を素早く計算するために、Somalierは正確な塩基をアッセイする。抽出ステップは、bam/cramファイルから1サンプルずつ直接行われる。

 relateステップは、extractコマンドの出力に基づいて実行される。新しいサンプルを追加して比較できるように、非常に高速に実行される。hom-ref, het, hom-altにはサンプルごとに3つのビットベクトルを使用する。各ビットベクタは64ビットの整数の配列で、サンプル中のインデックスのバリアントが例えばヘテロ接合体である場合、各ビットが設定される。この設定では、高速なビット演算とpopcountハードウェア命令を使用して、関連性を非常に迅速に計算することができる。

 

 

インストール

Github

リリースからstatic bainaryをダウンロードする。

./somalier -h

# ./somalier -h

somalier version: 0.2.10

Commands: 

  extract      :   extract genotype-like information for a single sample from VCF/BAM/CRAM.

  relate       :   aggregate `extract`ed information and calculate relatedness among samples.

  ancestry     :   perform ancestry prediction on a set of samples, given a set of labeled samples

  find-sites   :   create a new sites.vcf.gz file from a population VCF (this is rarely needed).

./somalier extract -h

# ./somalier extract -h

somalier version: 0.2.10

somalier extract

 

extract genotype-like information for a single-sample at selected sites

 

Usage:

  somalier extract [options] sample_file

 

Arguments:

  sample_file      single-sample CRAM/BAM/GVCF file or multi/single-sample VCF from which to extract

 

Options:

  -s, --sites=SITES          sites vcf file of variants to extract

  -f, --fasta=FASTA          path to reference fasta file

  -d, --out-dir=OUT_DIR      path to output directory (default: .)

  --sample-prefix=SAMPLE_PREFIX

                             prefix for the sample name stored inside the digest

  -h, --help                 Show this help

./somalier relate -h

# ./somalier relate -h

somalier version: 0.2.10

somalier relate

 

calculate relatedness among samples from extracted, genotype-like information

 

Usage:

  somalier relate [options] [extracted ...]

 

Arguments:

  [extracted ...]  $sample.somalier files for each sample. the first 10 are tested as a glob patterns

 

Options:

  -g, --groups=GROUPS        optional path  to expected groups of samples (e.g. tumor normal pairs).

specified as comma-separated groups per line e.g.:

    normal1,tumor1a,tumor1b

    normal2,tumor2a

  --sample-prefix=SAMPLE_PREFIX

                             optional sample prefixes that can be removed to find identical samples. e.g. batch1-sampleA batch2-sampleA

  -p, --ped=PED              optional path to a ped/fam file indicating the expected relationships among samples.

  -d, --min-depth=MIN_DEPTH  only genotype sites with at least this depth. (default: 7)

  --min-ab=MIN_AB            hets sites must be between min-ab and 1 - min_ab. set this to 0.2 for RNA-Seq data (default: 0.3)

  -u, --unknown              set unknown genotypes to hom-ref. it is often preferable to use this with VCF samples that were not jointly called

  -i, --infer                infer relationships (https://github.com/brentp/somalier/wiki/pedigree-inference)

  -o, --output-prefix=OUTPUT_PREFIX

                             output prefix for results. (default: somalier)

  -h, --help                 Show this help

./somalier ancestry -h

# ./somalier ancestry -h

somalier version: 0.2.10

somalier pca

 

dimensionality reduction

 

Usage:

  somalier pca [options] [extracted ...]

 

Arguments:

  [extracted ...]  $sample.somalier files for each sample. place labelled samples first followed by '++' then *.somalier for query samples

 

Options:

  --labels=LABELS            file with ancestry labels

  -o, --output-prefix=OUTPUT_PREFIX

                             prefix for output files (default: somalier-ancestry)

  --n-pcs=N_PCS              number of principal components to use in the reduced dataset (default: 5)

  --nn-hidden-size=NN_HIDDEN_SIZE

                             shape of hidden layer in neural network (default: 16)

  --nn-batch-size=NN_BATCH_SIZE

                             batch size fo training neural network (default: 32)

  --nn-test-samples=NN_TEST_SAMPLES

                             number of labeled samples to test for NN convergence (default: 101)

  -h, --help                 Show this help

./somalier find-sites -h

# ./somalier find-sites -h

somalier version: 0.2.10

somalier find-sites

 

Usage:

  somalier find-sites [options] vcf

 

Arguments:

  vcf              population VCF to use to find sites

 

Options:

  -x, --exclude=EXCLUDE      optional exclude files

  -i, --include=INCLUDE      optional include file. only consider variants that fall in ranges within this file

  --gnotate-exclude=GNOTATE_EXCLUDE

                             sites in slivar gnotation (zip) format to exclude

  --snp-dist=SNP_DIST        minimum distance between autosomal SNPs to avoid linkage (default: 10000)

  --min-AN=MIN_AN            minimum number of alleles (AN) at the site. (must be less than twice number of samples in the cohort) (default: 115_000)

  -h, --help                 Show this help

this will write output sites to: ./sites.vcf.gz

 

 

実行方法

1、extract  - extract genotype-like information for a single sample from VCF/BAM/CRAM.

cohort vcf or single samppleのvcfから抽出する。vcfはgzip圧縮してindexファイルが存在していること。

somalier extract -d extracted/ -s sites.vcf.gz -f g1k_v37_decoy.fa $cohort.vcf.gz
  • -s     sites vcf file of variants to extract
  • -f      path to reference fasta file
  • -d     path to output directory (default: .)

 

またはbam/cramファイル群から抽出する。forループで反復処理する。

for f in *.cram; do
somalier extract -d extracted/ -s sites.vcf.gz -f g1k_v37_decoy.fa $f
done

 

 

2、relate  - aggregate `extract`ed information and calculate relatedness among samples.

1の出力を指定する。任意でpedigreeファイルを指定する。

somalier relate -p $pedigree extracted/*.somalier
  • -p    optional path to a ped/fam file indicating the expected relationships among samples.

インタラクティブなHTML output(example)が出力される。

引用

Somalier: rapid relatedness estimation for cancer and germline studies using efficient genome sketches
Pedersen BS, Bhetariya PJ, Brown J, Marth G, Jensen RL, Bronner MP, Underhill HR, Quinlan AR
BioRxiv, 12 Nov 2019

 

ゲノム配列からウィルス配列を同定してアノテーションをつける VIBRANT

 

 細菌や古細菌に感染するウイルスは世界的に豊富であり、ほとんどの環境で宿主の数を上回っている [ref.1,2,3]。ウイルスは、感染時に宿主細胞の代謝状態を再プログラムすることができる義務的な細胞内病原性遺伝要素であり、多様な環境下で毎日20~40%の微生物を溶解させる可能性がある[ref.4, 5]。ウイルスは、その豊富さと広範な活動により、炭素、窒素、リン、硫黄などの必須栄養素の循環に寄与しているため、微生物群集の重要な要素となっている[ref.6,7,8,9,10]。ヒトのシステムでは、ウイルスは炎症性腸疾患などの様々な疾患を引き起こす可能性のある共生障害に関与していることが示唆されており、免疫システムとの共生的な役割を果たしていることさえある[ref.11,12,13]。

 ウイルスは、遺伝子の内容、配列、コードされた機能の多様性に大きな可能性を秘めている[ref.14,15,16,17]。その遺伝的多様性を認識し、新規抗微生物薬候補、バイオテクノロジー応用のための酵素バイオレメディエーションのためのこれらのウイルス配列を「マイニング」することに大きな関心が寄せられている[ref.18,19,20,21,22]。最近では、ウイルスが代謝プロセスを特異的に駆動することにより、栄養素の生物地球化学的循環を直接結びつける可能性があることがわかってきている[ref.23,24,25,26,27]。例えば、ウイルスは感染時に、宿主の代謝を引き継いで指示することで、必要な栄養素の40~90%を周囲の環境から獲得することができる[ref.28,29,30]。宿主の代謝フレームワークを操作するために、一部のウイルスは宿主から代謝遺伝子を選択的に「盗む」。これらの宿主由来の遺伝子は、補助代謝遺伝子(AMG)と総称され、感染時に積極的に発現し、ウイルスにフィットネス上の優位性を与えることができる [ref.31,32,33,34]。マイクロバイオームにおけるウイルスの役割を研究し、生態系機能のモデルにウイルスを組み込む必要性から、微生物群集全体の中でどのような配列がウイルスに由来するかを決定することが非常に重要になってきた。これらの配列には、遊離ウイルス、活動的な細胞内感染(これは、ある時点で全細菌の30%にも及ぶ可能性がある[ref.35])、粒子または宿主に付着したウイルス[ref.36]、および宿主に統合されたウイルスゲノムまたはエピソームウイルスゲノム(すなわち、プロウイルス)が含まれる。

(一部省略)

 より最近のツールはVirSorterの代替または補足ツールとして開発されてきた。VirFinder [ref.41]は、機械学習を実装し、ウイルス予測のアノテーションデータベースから完全に独立した最初のツールであり、後にPPR-Meta [ref.42]で実装されたプラットフォームであった。VirFinderは、ウイルスが8ヌクレオチド頻度の特徴的なパターン(8-merと呼ばれる)を示す傾向があることを考慮して構築されたが、これはウイルスが宿主と非常に類似したヌクレオチドパターンを共有できるという知識にもかかわらず提案されたものである[ref.43]。これらの8-merパターンは、500 bp以下の短い配列を迅速に分類し、モデルに基づいたスコアを生成するために使用されるが、スコアのカットオフを定義するのはユーザー次第である。VirFinderはVirSorterと比較してウイルスの回収能力を大幅に向上させることが示されたが、機械学習モデルのトレーニング中におそらくはリファレンスデータベースに関連したバイアスから、多様なウイルスを予測する際にホストとソース環境にかなりのバイアスがかかっていることも示された[ref.41]。さらに、特定の環境からのウイルスの回収率が低いことも確認されている[ref.44]。

 これまでのところ、VirSorterは、メタゲノムアセンブリ内の統合されたプロウイルスを同定するための最も効率的なツールである。他のツール、主にPHASTER [ref.47]とProphage Hunter [ref.48]は、メタゲノムアセンブリによって生成されたscaffoldsではなく、ゲノム全体からintegratiojnされたプロウイルスを同定することに特化したツールである。VirSorterと同様に、これら2つのプロウイルス予測器は、ウイルスに属する宿主ゲノムの領域を特定するために、スライドウィンドウを備えたリファレンス相同性とウイルス配列シグネチャに依存している。これらはゲノム全体には有用であるが、 溶菌性(すなわち非組み込み)ウイルスに属するscaffoldsを同定する能力はなく、大規模なデータセットでは動作が遅い。さらに、PHASTERとProphage Hunterはどちらもウェブベースのサーバーとしてのみ利用可能であり、スタンドアロンコマンドラインツールは提供されていない。

 ここで著者らはVIBRANT(Virus Identification By iteRative ANnoTation)を開発した。VIBRANTは、メタゲノムアセンブリとゲノム配列から遊離のウイルスとintegratedウイルスの両方を自動で回収、アノテーション、キュレーションするためのツールである。VIBRANTは、細菌と古細菌の両方に感染する多様なdsDNA、ssDNA、RNAウイルスを同定することができる。VIBRANTは、隠れマルコフモデル(HMM)を用いた非リファレンスベースの類似性検索から得られたタンパク質アノテーションシグネチャニューラルネットワークと独自の「v-score」指標を用いて、多様なウイルスや新規ウイルスの同定を最大化する。ウイルスを特定した後、VIBRANTは予測を検証するためのキュレーションステップを実行する。VIBRANTはさらに、AMGをハイライトすることでウイルス群集の機能を特徴づけ、ウイルス群集に存在する代謝パスウェイを評価する。すべてのウイルスゲノム、タンパク質、アノテーション代謝プロファイルは、ユーザーが使いやすい下流の解析と可視化のためのフォーマットにコンパイルされている。リファレンスウイルス、非リファレンスウイルスデータセット、様々なメタゲノムに適用した場合、VIBRANTはウイルスの回収率を最大化し、誤発見を最小限に抑えるという点で、VirFinder、VirSorter、MARVELを上回った。PHASTER、Prophage Hunter、VirSorterと比較した場合、VIBRANTは宿主scaffoldsからintegrationされたプロウイルス領域を抽出する能力において、同等の性能を示した。VIBRANTはまた、様々な環境に由来するウイルス間の代謝能力の違いを特定するためにも使用された。クローン病患者の3つのコホートにVIBRANTを適用したところ、健常者と比較してウイルスの種類が異なるだけでなく、病気の状態に影響を与えると考えられるウイルスがコードされた遺伝子を同定することができた。VIBRANTは https://github.com/AnantharamanLab/VIBRANT から無料でダウンロードできる。また、VIBRANTはCyVerse Discovery Environment(https://de.cyverse)を介して、ユーザーフレンドリーなウェブベースのアプリケーションとしても利用できる。

 

 

 

f:id:kazumaxneo:20200614172652p:plain

flowchart. Githubより

 

 

インストール

依存をcondaで導入し、VIBRANTのコードをpullしてテストした。テストしていないが、現在はbiocondaからインストールできる。

依存

VIBRANT has been tested and successfully run on Mac, Linux and Ubuntu systems.

Python Dependancies

本体 Github

#python3.7の仮想環境を作って入れる
conda create -n Vibrant -c bioconda -y hmmer Prodigal biopython python=3.7
conda activate Vibrant
conda install -c anaconda -y pandas matplotlib numpy seaborn scikit-learn

git clone https://github.com/AnantharamanLab/VIBRANT.git
chmod -R 777 VIBRANT
cd VIBRANT/

#or Bioconda(link)
conda install -c bioconda vibrant -y

python3 VIBRANT_run.py -h

$ python3 VIBRANT_run.py -h

usage: VIBRANT_run.py [-h] [--version] -i I [-f {prot,nucl}] [-t T] [-l L]

                      [-o O] [-virome] [-no_plot] [-k K] [-p P] [-v V] [-e E]

                      [-a A] [-c C] [-n N] [-s S] [-m M] [-g G]

 

Usage: python3 VIBRANT_run.py -i <input_file> [options]. VIBRANT identifies

bacterial and archaeal viruses (phages) from assembled metagenomic scaffolds

or whole genomes, including the excision of integrated proviruses. VIBRANT

also performs curation of identified viral scaffolds, estimation of viral

genome completeness and analysis of viral metabolic capabilities.

 

optional arguments:

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

  --version       show program's version number and exit

  -i I            input fasta file

  -f {prot,nucl}  format of input [default="nucl"]

  -t T            number of parallel VIBRANT runs, each occupies 1 CPU

                  [default=1, max of 1 CPU per scaffold]

  -l L            length in basepairs to limit input sequences [default=3000,

                  can increase but not decrease]

  -o O            number of ORFs per scaffold to limit input sequences

                  [default=4, can increase but not decrease]

  -virome         use this setting if dataset is known to be comprised mainly

                  of viruses. More sensitive to viruses, less sensitive to

                  false identifications [default=off]

  -no_plot        suppress the generation of summary plots [default=off]

  -k K            path to KEGG HMMs (if moved from default location)

  -p P            path to Pfam HMMs (if moved from default location)

  -v V            path to VOG HMMs (if moved from default location)

  -e E            path to plasmid HMMs (if moved from default location)

  -a A            path to viral-subset Pfam HMMs (if moved from default

                  location)

  -c C            path to VIBRANT categories file (if moved from default

                  location)

  -n N            path to VIBRANT annotation to name file (if moved from

                  default location)

  -s S            path to VIBRANT summary of KEGG metabolism file (if moved

                  from default location)

  -m M            path to VIBRANT neural network machine learning model (if

                  moved from default location)

  -g G            path to VIBRANT AMGs file (if moved from default location)

python3 scripts/VIBRANT_annotation.py -h

$ python3 scripts/VIBRANT_annotation.py -h

usage: VIBRANT_annotation.py [-h] [--version] -i I [-f {prot,nucl}] [-l L]

                             [-o O] [-virome] [-k K] [-p P] [-v V] [-e E]

                             [-a A] [-c C] [-n N] [-m M] [-g G]

 

See main wrapper script: VIBRANT_run.py. This script performs the bulk of the

work but is not callable on its own.

 

optional arguments:

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

  --version       show program's version number and exit

  -i I            input fasta file

  -f {prot,nucl}  format of input [default="nucl"]

  -l L            length in basepairs to limit input sequences [default=3000,

                  can increase but not decrease]

  -o O            number of ORFs per scaffold to limit input sequences

                  [default=4, can increase but not decrease]

  -virome         use this setting if dataset is known to be comprised mainly

                  of viruses. More sensitive to viruses, less sensitive to

                  false identifications [default=off]

  -k K            path to KEGG HMMs (if moved from default location)

  -p P            path to Pfam HMMs (if moved from default location)

  -v V            path to VOG HMMs (if moved from default location)

  -e E            path to plasmid HMMs (if moved from default location)

  -a A            path to viral-subset Pfam HMMs (if moved from default

                  location)

  -c C            path to VIBRANT categories file (if moved from default

                  location)

  -n N            path to VIBRANT annotation to name file (if moved from

                  default location)

  -m M            path to VIBRANT neural network machine learning model (if

                  moved from default location)

  -g G            path to VIBRANT AMGs file (if moved from default location)

 

 

データベースの準備

データベースをダウンロードする。11GB程度のディスクスペースが必要。

mkdir databases && cd databases
python3 VIBRANT_setup.py

#verify
python3 VIBRANT_test_setup.py

 

Done. Several new databases are now in this folder.

 

VIBRANT should be ready to go. You can verify this by running VIBRANT_test_setup.py within this folder (databases/)

# verify 

$ python3 VIBRANT_test_setup.py 

VIBRANT v1.0.0 is good to go!

 

 

テストラン

アセンブリFASTAファイルを指定する(ヘッダにfragmentという文字があってはならない)。アミノ酸FASTAも使用可能だが、ソートのされ方にルールがある。詳しくはGithub参照。

git clone https://github.com/AnantharamanLab/VIBRANT.git
cd VIBRANT/example_data/
python3 ../VIBRANT_run.py -i ../example_data/mixed_example.fasta -t 12
  • -t    number of parallel VIBRANT runs, each occupies 1 CPU [default=1, max of 1 CPU per scaffold] 
  • -i     input fasta file
  • -f    format of input [default="nucl"]
  • -l   length in basepairs to limit input sequences [default=3000, can increase but not decrease]

出力

f:id:kazumaxneo:20200614172434p:plain

 

VIBRANT_figure_phages_mixed_example.pdf

f:id:kazumaxneo:20200614172253p:plain

VIBRANT_figure_sizes_mixed_example.pdf

f:id:kazumaxneo:20200614172301p:plain

VIBRANT_figure_quality_mixed_example.pdf

f:id:kazumaxneo:20200614172332p:plain

VIBRANT_figure_pathways_mixed_example.pdf

f:id:kazumaxneo:20200614172349p:plain

(KEGG metabolic pathwayに基づく)

出力についてはGithubに置いてあるoutput_explanations.pdfを確認して下さい。詳しく説明されています。CyVerseでVIBRANTを使う手順についても簡単に書かれています。

引用

VIBRANT: automated recovery, annotation and curation of microbial viruses, and evaluation of viral community function from genomic sequences

Kristopher Kieft, Zhichao Zhou, Karthik Anantharaman
Microbiome volume 8, Article number: 90 (2020)

  

関連