macでインフォマティクス

macでインフォマティクス

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

fastaを分割するUCSCの fasplitコマンド

 

タイトルの通りのコマンド。

 

インストール

Home

#conda (link)
mamba install -c bioconda ucsc-fasplit -y

faSplit

$ faSplit

faSplit - Split an fa file into several files.

usage:

   faSplit how input.fa count outRoot

where how is either 'about' 'byname' 'base' 'gap' 'sequence' or 'size'.

Files split by sequence will be broken at the nearest fa record boundary.

Files split by base will be broken at any base.

Files broken by size will be broken every count bases.

 

Examples:

   faSplit sequence estAll.fa 100 est

This will break up estAll.fa into 100 files

(numbered est001.fa est002.fa, ... est100.fa

Files will only be broken at fa record boundaries

 

   faSplit base chr1.fa 10 1_

This will break up chr1.fa into 10 files

 

   faSplit size input.fa 2000 outRoot

This breaks up input.fa into 2000 base chunks

 

   faSplit about est.fa 20000 outRoot

This will break up est.fa into files of about 20000 bytes each by record.

 

   faSplit byname scaffolds.fa outRoot/

This breaks up scaffolds.fa using sequence names as file names.

       Use the terminating / on the outRoot to get it to work correctly.

 

   faSplit gap chrN.fa 20000 outRoot

This breaks up chrN.fa into files of at most 20000 bases each,

at gap boundaries if possible.  If the sequence ends in N's, the last

piece, if larger than 20000, will be all one piece.

 

Options:

    -verbose=2 - Write names of each file created (=3 more details)

    -maxN=N - Suppress pieces with more than maxN n's.  Only used with size.

              default is size-1 (only suppresses pieces that are all N).

    -oneFile - Put output in one file. Only used with size

    -extra=N - Add N extra bytes at the end to form overlapping

pieces.  Only used with size.

    -out=outFile Get masking from outfile.  Only used with size.

    -lift=file.lft Put info on how to reconstruct sequence from

                   pieces in file.lft.  Only used with size and gap.

    -minGapSize=X Consider a block of Ns to be a gap if block size >= X.

                  Default value 1000.  Only used with gap.

    -noGapDrops - include all N's when splitting by gap.

    -outDirDepth=N Create N levels of output directory under current dir.

                   This helps prevent NFS problems with a large number of

                   file in a directory.  Using -outDirDepth=3 would

                   produce ./1/2/3/outRoot123.fa.

    -prefixLength=N - used with byname option. create a separate output

                   file for each group of sequences names with same prefix

                   of length N.

 

 

実行方法

faSplit base

指定したFASTAを10個のファイルに分割する。出力prefixはoutputとする。

faSplit base input.fa 10 output

 

faSplit size

指定した(Multi-)FASTAを100bpの長さのチャンクに分割する。

faSplit size input.fa 100 output

 

指定した(Multi-)FASTAを10000bpの長さ(チャンク)に分割する。Nが1000以上ある配列は捨てる。

faSplit size input.fa 10000 output -maxN=1000
  • -maxN=N - Suppress pieces with more than maxN n's.  Only used with size.
                  default is size-1 (only suppresses pieces that are all N).

 

指定した(Multi-)FASTAを10000bpの長さで分割する。オーバーラップする領域を500bpずつ加える。出力は1ファイルにする。

faSplit size input.fa 10000 output.fasta -extra=500 -oneFile
  • -extra=N - Add N extra bytes at the end to form overlapping
    pieces.  Only used with size.
  • -oneFile - Put output in one file. Only used with size

 

split about

指定した(Multi-)FASTAをだいたい20000byteのサイズになるよう分割する。

faSplit about input.fa 20000 output
  • -maxN=N - Suppress pieces with more than maxN n's.  Only used with size.

 

faSplit byname

指定した(Multi-)FASTAを配列名で分割する。

faSplit byname input.fa output

 

faSplit gap

指定した(Multi-)FASTAをギャップ部分で分割する。

faSplit gap input.fa output
  • -minGapSize=X  Consider a block of Ns to be a gap if block size >= X.
                      Default value 1000.  Only used with gap.
  •  -noGapDrops - include all N's when splitting by gap.

 

参考

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

メタゲノムのハイブリッドアセンブリとビニングのためのベスト・プラクティス・パイプライン nf-core/mag

2023/03/02 論文引用

 

 ショットガンメタゲノムデータを解析することで、微生物群集に関する貴重な知見が得られると同時に、個々のゲノムレベルでの解決が可能となる。しかし、完全なリファレンスゲノムが存在しない場合、シークエンスリードからメタゲノムアセンブルゲノム(MAG)を再構築する必要がある。本研究では、メタゲノムアセンブリ、ビニング、分類学的分類を行うnf-core/magパイプラインを紹介する。nf-core/magは、ショートリードとロングリードを組み合わせることでアセンブリの連続性を高め、サンプルごとのグループ情報を共アセンブリやゲノムビニングに利用することができる。パイプラインは、インストールが容易で、すべての依存関係がコンテナ内に用意されており、移植性と再現性に優れている。Nextflowで書かれており、パイプライン開発のベストプラクティスであるnf-coreイニシアチブの一環として開発されている。すべてのコードは、GitHubのnf-core organization(https://github.com/nf-core/mag)でホストされており、MITライセンスで公開されている。

 

usage

https://nf-co.re/mag/usage

 

Githubより

デフォルトでは、パイプラインは次の解析を実行する。ショートリードとロングリードの両方をサポートしている。

1、fastpとPorechopでリードとアダプターをクオリティートリムし、FastQCで基本的なQCを実行する。

2、Centrifugeおよび/またはKraken2を用いてリードにtaxonomyを割り当てる。
3、MEGAHITとSPAdesを用いてアセンブリを行い、Quastを用いて品質をチェックする。

4、MetaBAT2を用いてビニングを行い、Buscoを用いてゲノムビンの品質を確認する。

5、GTDB-TkやCATを用いてビンに分類を付与する。

6、指定されたresultsディレクトリに、結果の一部やソフトウェアのバージョンをまとめたMultiQCのレポートなどを作成する。

 

2023/03/02

 

インストール

依存

  • Nextflow (>=21.04.0)

Github

 

テストラン

conda、docker、Singularity、Shifter、Podman(Docker互換のコンテナエンジン)、Charliecloudなどに対応している。

#docker
nextflow run nf-core/mag -profile test,docker

#conda
nextflow run nf-core/mag -profile test,conda

 

出力

f:id:kazumaxneo:20210905223511p:plain

Taxonomy

f:id:kazumaxneo:20210905223720p:plain

Assembly

f:id:kazumaxneo:20210905223739p:plain

Genome Binning

f:id:kazumaxneo:20210905223806p:plain

MEGAHIT-test_minigut-binDepths.heatmap.png

f:id:kazumaxneo:20210905223842p:plain

SPAdes-test_minigut-binDepths.heatmap.png

f:id:kazumaxneo:20210905223917p:plain

Genome Binning/QC

f:id:kazumaxneo:20210905224044p:plain

multiqc

f:id:kazumaxneo:20210905224139p:plain

 

 

実際のランではprofileとfastqのパス、もしくはfastqのパスとサンプル名を記載したCSVファイルを指定する。

#docker
nextflow run nf-core/mag -profile docker --input '*_R{1,2}.fastq.gz'

#samplesheet.csv
nextflow run nf-core/mag -profile docker --input samplesheet.csv

カンマ区切りで最大5列の情報を記載する。ヘッダーはsample,group,short_reads_1,short_reads_2,long_readsとする。

sample,group,short_reads_1,short_reads_2,long_reads
sample1,0,data/sample1_R1.fastq.gz,data/sample1_R2.fastq.gz,data/sample1.fastq.gz
sample2,0,data/sample2_R1.fastq.gz,data/sample2_R2.fastq.gz,data/sample2.fastq.gz
sample3,1,data/sample3_R1.fastq.gz,data/sample3_R2.fastq.gz,



サンプルIDは一意でなければならない。2列目のグループ情報は、ビニングステップの共分散の計算にのみ使用され、共アセンブリには使用されない。共アセンブリには--coassemble_groupオプションを使う。3列目以降で指定するFastQファイルは圧縮されている必要がある(.fastq.gz, .fq.gz)。ロングリードもある場合、ペアエンドのshort readデータとの組み合わせでのみ提供可能。1つのサンプルシート内でシングルエンドとペアエンドの混在は不可。シングルエンドリードを指定する場合は、コマンドラインパラメータ -single_end も指定する。

 

引用

nf-core/mag: a best-practice pipeline for metagenome hybrid assembly and binning

Sabrina Krakau,  Daniel Straub,  Hadrien Gourlé,  Gisela Gabernet,  Sven Nahnsen

bioRxiv, Posted August 31, 2021

 

2023/01

nf-core/mag: a best-practice pipeline for metagenome hybrid assembly and binning
Sabrina Krakau, Daniel Straub, Hadrien Gourlé, Gisela Gabernet, and Sven Nahnsen

NAR Genom Bioinform. 2022 Mar; 4(1)

 

 

参考

file:///Users/kazu/Downloads/IPSJ-BIO18054047.pdf

 

DockerユーザーのためのPodmanとBuildahの紹介 - 赤帽エンジニアブログ

(メタゲノム)BAMのカバレッジ、polymorphic サイト率、リファレンスフリーのコンセンサス配列を計算する CMSeq

 

CMSeqは、SegataLabで公開されている、リファレンスのカバレッジ、polymorphic サイト率、BAMからのコンセンサス配列計算のための.bamファイルへのインターフェースを提供するコマンド群。

 

インストール

 

依存

  • Requires:
  • samtools (> 1.x)
  • numpy
  • pysam
  • pandas
  • Biopython with bcbio-gff module (warning: Biopython <= 1.76 is required)

Github

#conda (link)
mamba install -c bioconda cmseq -y

#pip(pypi)
pip install CMSeq

> breadth_depth.py -h

$ breadth_depth.py -h

usage: breadth_depth.py [-h] [-c REFERENCE ID] [-f] [--sortindex]

                        [--minlen MINLEN] [--minqual MINQUAL]

                        [--mincov MINCOV] [--truncate TRUNCATE] [--combine]

                        BAMFILE

 

calculate the Breadth and Depth of coverage of BAMFILE.

 

positional arguments:

  BAMFILE               The file on which to operate

 

optional arguments:

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

  -c REFERENCE ID, --contig REFERENCE ID

                        Focus on a subset of references in the BAM file. Can

                        be a list of references separated by commas or a FASTA

                        file (the IDs are used to subset)

  -f                    If set unmapped (FUNMAP), secondary (FSECONDARY), qc-

                        fail (FQCFAIL) and duplicate (FDUP) are excluded. If

                        unset ALL reads are considered (bedtools genomecov

                        style). Default: unset

  --sortindex           Sort and index the file

  --minlen MINLEN       Minimum Reference Length for a reference to be

                        considered

  --minqual MINQUAL     Minimum base quality. Bases with quality score lower

                        than this will be discarded. This is performed BEFORE

                        --mincov. Default: 30

  --mincov MINCOV       Minimum position coverage to perform the polymorphism

                        calculation. Position with a lower depth of coverage

                        will be discarded (i.e. considered as zero-coverage

                        positions). This is calculated AFTER --minqual.

                        Default: 1

  --truncate TRUNCATE   Number of nucleotides that are truncated at either

                        contigs end before calculating coverage values.

  --combine             Combine all contigs into one giant contig and report

                        it at the end

> polymut.py -h

$ polymut.py -h

usage: polymut.py [-h] [-c REFERENCE ID] [-f] [--sortindex] [--minlen MINLEN]

                  [--minqual MINQUAL] [--mincov MINCOV]

                  [--dominant_frq_thrsh DOMINANT_FRQ_THRSH]

                  [--gff_file GFF_FILE]

                  BAMFILE

 

Reports the polymorpgic rate of each reference (polymorphic bases / total

bases). Focuses only on covered regions (i.e. depth >= 1)

 

positional arguments:

  BAMFILE               The file on which to operate

 

optional arguments:

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

  -c REFERENCE ID, --contig REFERENCE ID

                        Focus on a subset of references in the BAM file. Can

                        be a list of references separated by commas or a FASTA

                        file (the IDs are used to subset)

  -f                    If set unmapped (FUNMAP), secondary (FSECONDARY), qc-

                        fail (FQCFAIL) and duplicate (FDUP) are excluded. If

                        unset ALL reads are considered (bedtools genomecov

                        style). Default: unset

  --sortindex           Sort and index the file

  --minlen MINLEN       Minimum Reference Length for a reference to be

                        considered. Default: 0

  --minqual MINQUAL     Minimum base quality. Bases with quality score lower

                        than this will be discarded. This is performed BEFORE

                        --mincov. Default: 30

  --mincov MINCOV       Minimum position coverage to perform the polymorphism

                        calculation. Position with a lower depth of coverage

                        will be discarded (i.e. considered as zero-coverage

                        positions). This is calculated AFTER --minqual.

                        Default:1

  --dominant_frq_thrsh DOMINANT_FRQ_THRSH

                        Cutoff for degree of `allele dominance` for a position

                        to be considered polymorphic. Default: 0.8

  --gff_file GFF_FILE   GFF file used to extract protein-coding genes

> poly.py -h

$ poly.py -h

usage: poly.py [-h] [-c REFERENCE ID] [-f] [--sortindex] [--minlen MINLEN]

               [--minqual MINQUAL] [--mincov MINCOV] [--pvalue PVALUE]

               [--seq_err SEQ_ERR] [--dominant_frq_thrsh DOMINANT_FRQ_THRSH]

               BAMFILE

 

Reports the polymorpgic rate of each reference (polymorphic bases / total

bases). Focuses only on covered regions (i.e. depth >= 1)

 

positional arguments:

  BAMFILE               The file on which to operate

 

optional arguments:

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

  -c REFERENCE ID, --contig REFERENCE ID

                        Focus on a subset of references in the BAM file. Can

                        be a list of references separated by commas or a FASTA

                        file (the IDs are used to subset)

  -f                    If set unmapped (FUNMAP), secondary (FSECONDARY), qc-

                        fail (FQCFAIL) and duplicate (FDUP) are excluded. If

                        unset ALL reads are considered (bedtools genomecov

                        style). Default: unset

  --sortindex           Sort and index the file

  --minlen MINLEN       Minimum Reference Length for a reference to be

                        considered. Default: 0

  --minqual MINQUAL     Minimum base quality. Bases with quality score lower

                        than this will be discarded. This is performed BEFORE

                        --mincov. Default: 30

  --mincov MINCOV       Minimum position coverage to perform the polymorphism

                        calculation. Position with a lower depth of coverage

                        will be discarded (i.e. considered as zero-coverage

                        positions). This is calculated AFTER --minqual.

                        Default:1

  --pvalue PVALUE       Binomial p-value threshold for the binomal-polymorphic

                        test. Default: 0.01

  --seq_err SEQ_ERR     Sequencing error rate. Default: 0.001

  --dominant_frq_thrsh DOMINANT_FRQ_THRSH

                        Cutoff for degree of `allele dominance` for a position

                        to be considered polymorphic. Default: 0.8

 

 

実行方法

breadth_depth.py  BAMアライメントファイルから各リファレンスのカバレッジ範囲の広さ、平均値、中央値を表にして出力

breadth_depth.py sorted.bam

#chrを指定
breadth_depth.py ―c chr1 sorted.bam

#unsoorted.bam, minimum depth 10
breadth_depth.py --sortindex --mincov 10 sorted.bam

 

 

polymut.py  -  タンパク質をコードする遺伝子上のpolymorphic サイト率を計算

polymut.py sorted.bam

#unsoorted.bam, minimum coverage 10, minimum quality 30
poly.py --sortindex --mincov 10 --minqual 30 unsort.bam

Githubより

この関数は、タンパク質をコードする遺伝子について、ヌクレオチドレベルでdominant およびsecond-dominant alleles を考慮し、ORFをタンパク質に変換し、タンパク質のdominant およびsecond-dominant間の(タンパク質レベルでの)同義的および非同義的な突然変異の数を計算して出力する。 Dominantとsecond-dominant のカバレッジの比率がdominant_frq_thrshよりも小さい位置は非変異とみなされる。この関数は、Pasolli et al., 2019において、メタゲノムにおける株の異質性を計算するためのアドホックな尺度として使用された。

 

 

consensus.py  -  BAMからReference Freeのコンセンサス配列を出力

consensus.py sort.bam > consensus.fasta

結果はFASTA形式で出力される。再構成された配列の長さは、リファレンスのオリジナルの長さに拘束される。その長さでは、すべての位置がカバーされていない可能性があり、次のような場合に起こる。

  • その位置に対応するリードが存在しない
  • その位置にマッピングされるリードの数が少なすぎる(つまりmincov未満)。
  • その位置にマッピングされているリードの品質が低い(すなわち、< minqual)。
  • その位置のヌクレオチドの分布に問題がある可能性がある(例:dominant_allele_frequency < dominant_frq_thrsh):この場合、ノイズを減らすために、その位置は除外される。

 

consensus_aDNA.py 

古代メタゲノム研究を想定して、BAMファイルからゲノムのコンセンサスを抽出する。カバレッジが5より小さく、ダメージ確率(MapDamage2のStats_out_MCMC_correct_prob.csv)が0.95より大きい位置は無視する。

consensus_aDNA.py --mincov 5 -r reference.fna --pos_specific_prob_tab Stats_out_MCMC_correct_prob.csv --pos_damage_prob_thrsh 0.95 sorted.bam

 

引用

https://github.com/SegataLab/cmseq#breadth-and-depth-of-coverage-with-breadth_depthpy

 

関連


 

タンパク質配列を使って ロングリードのフレームシフトエラー修正を行う Proovframe

 

 精度は向上しているものの、ロングリードデータの基本的な遺伝子予測は、small indelsから生じるフレームシフトによって損なわれることが多い。相補的なショートリードやロングリードを用いたコンセンサスポリッシュは、この影響を軽減することができるが、普遍的に深いシーケンスデプスが必要であり、コミュニティメンバーの大半が希少であるような複雑なサンプルでは達成が困難である。ここでは、ロングリードアセンブリや生のロングリードにおけるフレームシフトエラーを克服するための代替アプローチを実装したソフトウェアであるproovframeを紹介する。proovframeは、リファレンスデータベースに対するタンパク質とヌクレオチドのアラインメントを利用して、コンティグやリードのindelをピンポイントで特定し、1-2塩基を削除または挿入して修正することで、アラインメントされた領域のリーディングフレームの忠実度を保守的に回復させる。シミュレーションと実世界のベンチマークデータを用いて、proovframeは、アセンブルされたデータに対してショートリードベースのポリッシングと同等の性能を発揮すること、遠いタンパク質のホモログに対してもうまく機能すること、さらには生のリードにも直接適用できることを示した。これらの結果から、タンパク質によるフレームシフト補正は、一般的なポリッシング戦略と組み合わせても、またそれに代えても、ロングリードデータの解析性を大幅に向上させることが実証された。Proovframe は https://github.com/thackl/proovframe から入手できる。

 

Githubより

遠い関係にあるガイドタンパク質でも良好な結果が得られ、アミノ酸の同一性が60%未満のセットでのテストに成功しています。アセンブリーの従来のコンセンサスポリッシングアプローチに加えて、追加のポリッシングステップとして使用することができます。コンセンサスポリッシングのためのシーケンス深度が不足しているデータ(例えば、環境メタゲノムサンプルから得られる多くの希少生物に共通する問題)に使用することができます。

 

インストール

提供されているバイナリをダウンロードした。

依存

  • Requires DIAMOND v2.0.3 or newer for mapping.

Github

git clone https://github.com/thackl/proovframe.git
cd proovframe/bin/ #パスを通す

>  ./proovframe

$ ./proovframe

proovframe-v0.9.7 - frame-shift correction for long-read genomics data

 

Please cite:

- Hackl T, Trigodet F, Murat Eren A, Biller SJ, Eppley JM, Luo E, Burger A, DeLong EF,

  Fischer MG, "proovframe: frameshift-correction for long-read (meta)genomics",

  bioRxiv. 2021. p. 2021.08.23.457338. doi:10.1101/2021.08.23.457338

- Buchfink B, Reuter K, Drost HG, "Sensitive protein alignments at tree-of-life scale

  using DIAMOND", Nature Methods 18, 366–368 (2021). doi:10.1038/s41592-021-01101-x

 

Usage: proovframe <cmd> [options]

 

  map     produce read/contig-to-protein alignments with diamond (convenience wrapper)

  fix     correct reads/contigs guided by diamond protein mappings

 

auxiliary commands:

  prf     profile error rates from reference alignments at nucleotide-level

 

See `proovframe cmd -h` for more options

> ./proovframe-map  

$ ./proovframe-map 

Usage: proovframe map [-a|-d proteins] -o alignments.o6 seqs.fa -- extra-diamond-params

 -a/--aa              protein file, not required if db provided.

 -d/--db              created if not existing and --aa given [basename(aa).dmnd]

 -o/--out             write alignments to this file [basename(seqs).o6]

 -t/--threads         number of CPU threads

 -m/--diamond-mode    one of fast,sensitive,{mid,more,very,ultra}-sensitive' [more-sensitive]

 -y/--dry-run         print the diamond command, but don't run it

 -h/--help            show this help

 -D/--debug           show debug messages

 

For consensus sequences with rather low expected error rates

and if your reference database has a good represention of similar

sequences, you might want to switch to '-m fast' or '-m sensitive'

to speed things up.

Also note, I've experienced inefficient parallelization if

correcting a small number of Mb sized genomes (as opposed to thousands

of long-reads) - presumably because diamond threads on a per-sequence

basis

 

 

実行方法

1、エラー率推定(任意)

proovframe prf draft.fasta long-read.fq

 

2、タンパク質のリードへのマッピング

タンパク質配列(*1)、出力名、最後にONTリード、もしくはドラフトアセンブリ配列を指定する。

proovframe map -a proteins.faa -o raw-seqs.tsv raw-seqs.fa

raw-seqs.tsvが出力される。

 

3、 リードのフレームシフトの修正

proovframe fix -o corrected-seqs.fa raw-seqs.fa raw-seqs.tsv

f:id:kazumaxneo:20210904010244p:plain

corrected-seqs.faが出力される。

 

引用

proovframe: frameshift-correction for long-read (meta)genomics
Thomas Hackl,  Florian Trigodet,  A. Murat Eren,  Steven J. Biller, John M. Eppley, Elaine Luo, Andrew Burger, Edward F. DeLong,  Matthias G. Fischer

bioRxiv, Posted August 24, 2021

 

*1

ベンチマークをみるとuniref90などが使われているようです。

 

 

複数の実験で得られた機能的な遺伝子アノテーションを階層的に整理し、視覚的にナビゲートする FunMappOne

 

 オミックスデータの解析において、遺伝子の機能アノテーションは不可欠なステップである。現在、遺伝子群の機能をオントロジーや分子パスウェイなどの高次の表現にまとめるためのデータベースや手法が複数存在する。オミックス実験の結果を機能カテゴリにアノテーションすることは、根本的な制御ダイナミクスを理解するだけでなく、複数の実験条件をより高い抽象度で比較するためにも不可欠である。オミックス実験の機能プロファイルを表現したり、比較したりするためのツールは、すでにいくつか公開されている。しかし、実験の数や濃縮された機能語彙の数が多い場合、グラフで表現しても結果を解釈するのが難しくなる。そのため、高次元であっても結果の解釈を容易にするために、アノテーションをグラフィカルにナビゲートし、さらに要約するインタラクティブでユーザーフレンドリーなツールが必要とされている。
 著者らは、複数の機能アノテーションが持つ固有の階層構造を利用して、エンリッチメント分析で得られた結果をより高いレベルの解釈に要約し、要約された各レベルで遺伝子関連情報をマッピングする手法を開発した。1つまたは複数の実験の機能アノテーションを一度に視覚化できる、ユーザーフレンドリーなグラフィカルインターフェースを構築した。このツールはFunMappOneというR-Shinyアプリケーションとして実装されており、https://github.com/grecolab/FunMappOne で公開されている。
 FunMappOneはR-shinyのグラフィカルなツールで、ヒトやマウスの遺伝子の複数のリストを入力し、オプションでそれらの関連する修飾の大きさを加え、Gene Ontology、Kyoto Encyclopedia of Genes and Genomes、またはReactomeデータベースからエンリッチされたアノテーションを計算し、合理的なグループに整理された機能語彙とパスウェイのインタラクティブなマップを報告する。FunMappOneは、複数の実験を迅速かつ簡便に比較することができ、結果を簡単に解釈することができる。

 

論文表1に、DAVID,Enrichr,ToppGene,g:profiler,clusterProfiler,Goplot,BACA、を含む遺伝子機能解析ツールとFunMappOneの比較を示している。他のツールのほとんどは、KEGGパスウェイ、Reactomeパスウェイ、Gene Ontologyを分析する機能を提供しており、エンリッチメント結果をグラフィックで表示することもできる。Goplotのみが遺伝子関連の値を語彙にマッピングする機能を提供しており、Enrichrとg:profilerはウェブベースのグラフィカルユーザーインターフェースを提供する唯一のツールである。他のツールでは、結果を要約したり、複数の実験から得られた機能プロファイルをクラスター化する機能は提供されていない。著者の知る限り、FunMappOneはこれらの機能を全てユーザーフレンドリーなグラフィックインターフェースで提供する唯一のツールである。

 

Github

 

使い方

Docker image

docker pull grecolab/funmappone:3.5.3
docker run --rm -p 8787:3838 grecolab/funmappone:3.5.3

http://localhost:8787にアクセスする。

 

f:id:kazumaxneo:20210816092357p:plain

https://funmappone.cloud.nanosolveit.eu/でも利用できる。

 

関心のある各実験条件に関連する遺伝子のリストと、それらの修飾情報(例えば、DEGsのfold change値)を含むスプレッドシートファイルを用意する。

example (issues#1)

FunMappOne/Mouse_Test_Data.xlsx at master · Greco-Lab/FunMappOne · GitHub

スプレッドシートには各実験条件のシートが含まれており、条件IDで名付けられている。

f:id:kazumaxneo:20210902003501p:plain

各シートには2つの列があり、それぞれ遺伝子の識別子(Entrez Gene、Gene Symbol、またはEnsembl gene ID)、および任意でそれらの情報を含む。

 

追加のシート"Grouping"も必要で、Groupingにはそれぞれ条件IDと条件グループ情報(番号でグループ化)の2つの列がある。

f:id:kazumaxneo:20210902003615p:plain

 

左上のOrganisumから生物を選択して、準備したスプレッドシートファイルをアップロードする。写真は読み込んだexample データ(mouseを選択している)。

f:id:kazumaxneo:20210902002112p:plain

 

続いて、機能アノテーションGene Ontology - BP、Gene Ontology - CC、Gene Ontology - MF、KEGG、Reactome)の選択、p値補正法とカットオフ値を選択する。発現変動の振幅(fold change, p-value)が提供されている場合、濃縮された語彙の要約値をその値に関連したカラースケールでプロットするか、3色(マイナス、ゼロ、プラス)のみでプロットするかを選択する(後者の機能は、変動の支配的なサインを強調する場合に有用)。

f:id:kazumaxneo:20210902004306p:plain

発現変動の振幅(fold change, p-value)値を含まない遺伝子リストのみをアップロードした場合、エンリッチされた語彙のエンリッチメントP値が表示される。最後にめにゅー右下のGenerate Mapをクリックする。

 

エンリッチされているかどうか計算される。しばらく時間がかかる。

 

出力

f:id:kazumaxneo:20210902011155p:plain

出力する項目は編集できる。

マップを表示するレベルを選択する。

f:id:kazumaxneo:20210902011340p:plain

 

サンプルを選択する。

f:id:kazumaxneo:20210902011355p:plain

 

編集すると、選択したレベルのカテゴリーに対応する行が,レベルの上位のクラスによってグループ化された新しいマップが表示される。

f:id:kazumaxneo:20210902011604p:plain

マップのセルの色は、カテゴリーの行に属する実験条件の列にある、エンリッチされたすべての語彙の値に関連付けられている。

 

似たようなプロファイルを持つ実験条件をクラスタリングして並び替えることもできる。クラスタリング法と距離を選択する。

f:id:kazumaxneo:20210902011753p:plain

クラスタリングタブに切り替えるとデンドログラム表示される。

f:id:kazumaxneo:20210902012250p:plain

 

 

Heatmap genesタブでは遺伝子を視覚化して調べることができる。

f:id:kazumaxneo:20210902012159p:plain

 

引用

FunMappOne: a tool to hierarchically organize and visually navigate functional gene annotations in multiple experiments
Giovanni Scala, Angela Serra, Veer Singh Marwah, Laura Aliisa Saarimäki, Dario Greco
BMC Bioinformatics volume 20, Article number: 79 (2019)

メタゲノムアセンブリの品質評価を行う DeepMAsED

 

 アセンブリーの品質を評価する手法の多くは、リファレンスゲノム(アセンブリを比較するためにキュレートされたゲノムのセット)を必要とする。そのような手法として、コンティグを1つ以上のリファレンスゲノムにマッピングして、逆位、リアレンジメント、種間転位などの広範なミスアセンブリを推定するmetaQUAST(Mikheenko et al.、2016)がある。それにもかかわらず、新たに発見されたいくつかのMAGは、ゲノムデータベースに代表される生物とは遠縁の新規生物を表しており、このことが有効なリファレンスの選択をしばしば不可能にしている。例えば、最も近い既知のrelativeが異なる細菌門に属していることがあり、これは非常に深い進化的分岐を示唆している(Pasolli et al.、2019)。逆に、バクテリアの種にはゲノムの違いが大きい株が含まれることがあるため、株間の真のゲノム変異がミスアセンブリと誤解されることがある。

 そのため、メタゲノムからのフルゲノムまたは部分ゲノムのアセンブリには、リファレンスフリーの手法が有効だが、そのような手法はほとんど存在しない。最も広く使われている手法はCheckMであり、品質評価によく使われる2つのサブスコア(completenessとcontamination)を提供している。これらの指標は、特定の遺伝子座のカウントに基づいているため、個々のコンティグのレベルでゲノムアセンブリを評価するものではない。それでも、Almeidaら(2019)やPasolliら(2019)などの大規模なMAGデータセットでは、CompletenessとConaminationの閾値を用いて、アセンブリの品質をスクリーニングしている。

 コンティグのミスアセンブリを評価する既存のリファレンスフリー手法は、ALE、SuRankCo、REAPR、VALETである。ALEは、ゲノムに関する確率的な仮定を行い、アセンブリされたコンティグの各ヌクレオチドに対して尤度スコアを提供するが、コンティグレベルでの品質推定値は提供しない。さらに、これらはゲノムアセンブリ用に開発されており、メタゲノムアセンブリ用ではない。SuRankCoは、入力コンティグの品質推定値を予測するために、手作業で作成された特徴のセットにランダムフォレスト分類法を利用することで、ゲノムやアセンブラに関する仮定を少なくしている。REAPRは、リードをコンティグにマッピングし、カバレッジやリードペアの適切なマッピングなど、様々なメトリクスからアンサンブルスコアを計算する。VALETはREAPRに類似しており、コンティグにマッピングされたリードに基づく複数のメトリクスの組み合わせを使用する。重要なのは、ALE、SuRankCo、REAPRはすでに保守されておらず、互換性の問題で使用できないということである。また、これらの手法の精度は、複雑なメタゲノムアセンブリでは評価されていない。

 ここでは、コンティグレベルでのメタゲノムアセンブリの品質評価のための機械学習システムであるDeepMAsED (mased: a Middle English term for ‘misled’; Deep Metagenome Assembly Error Detection) を紹介する。DeepMAsEDの目的は、リファレンスゲノムがなくてもmetaQUASTから広範なミスアセンブリのラベルを予測することである(注;metaQUASTからの広範なミサアセンブリラベルをDeepMAsEDのターゲットとして使用)。本論文の貢献度は以下のようにまとめられる。(以下省略)

 

インストール

Github

#conda (link)
mamba create -n deepmased bioconda::deepmased -y
conda activate deepmased

#pip(pypi)
pip install DeepMAsED

> DeepMAsED -h

$ DeepMAsED -h

usage: DeepMAsED [-h] [--version] {train,predict,evaluate,features} ...

 

DeepMAsED: Deep learning for Metagenome Assembly Error Detection

 

positional arguments:

{train,predict,evaluate,features}

 

optional arguments:

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

--version show program's version number and exit

 

DESCRIPTION:

Usage: DeepMAsED <subcommand> <subcommand_params>

Example: DeepMAsED train -h

 

For general info, see https://github.com/leylabmpi/DeepMAsED/1

 

 

実行方法

ランにはアセンブリfastaファイル(1-kb以上の配列だけ使用することが推奨されている)と、その配列にリードをマッピングして得たbamファイルが必要(論文中ではBowtie2 (v2.3.5)を使用)。bamファイルを得たら、fastaとbamの関係を示したタブ区切りテキストを作成する。ここでは1組だけなので、次のようなbam_fasta.tableファイルを準備した。

bam fasta

map.bam final.contigs.fa

 

1、feature tableの作成。準備したタブ区切りテキストを指定する。

DeepMAsED features bam_fasta.table

map_feats.tsvができる。

 

2、ミスアセンブリの検出。

DeepMAsED predict map_feats.tsv

 deepmased_predictions.tsvができる。

 

出力例(Githubより転載)

f:id:kazumaxneo:20210831234630p:plain
DeepMAsEDのスコアのどの値でミスアセンブリ配列をカットオフするかですが、論文に掲載されている図が区切りの参考になります。

 

引用

DeepMAsED: evaluating the quality of metagenomic assemblies
Olga Mineeva, Mateo Rojas-Carulla, Ruth E Ley, Bernhard Schölkopf, Nicholas D Youngblut
Bioinformatics, Volume 36, Issue 10, 15 May 2020, Pages 3011–3017

 

関連


 

 

大量のタンデムリピート構造を含むゲノムをインタラクティブに可視化する StainedGlass

2022/01/13 論文引用

 

 ドットプロット解析は、配列の同一性や方向性の違いのような複雑なリピートの基礎構造を明らかにするためによく用いられる。ロングリードシーケンス技術の進歩により、最近ではますます連続したリファレンスゲノムのアセンブリやヒトの完全な染色体、完全なセントロメア、タンデムデュプリケーション、その他のヘテロクロマティックアレイが利用でき、より複雑なリピート構造や 遺伝的変異を体系的に解析することができるようになった。これらの構造の大きさと複雑さは、しばしば、3つの理由から従来のドットプロット解析では対応ない;1) 現在の可視化手法は、主に完全一致またはk-mer一致に基づいており、セントロメアに見られる複雑な高次のリピートや、これらの大きなリピート間に予想されるミスマッチなどの複雑な高次の繰り返し構造には適していない。2)ドットプロットは、メガベースの配列データからなるタンデムアレイの解像度は限られており、しばしば黒い四角になってしまい ドットプロットは、メガベースの配列データからなるタンデムアレイでは解像度が低く、しばしば黒い四角になり、リピートのサイズと存在以外の情報はほとんど伝わらない。タンデムアレイで完全に一致するものを特定する場合、可能なペアワイズマッチの数が急激に増加する。3)タンデムアレイ(MUMmerなど)で完全に一致するものを特定する場合、可能なペアワイズマッチの数が急速に増加し、より分岐したものを比較するために小さな最小マッチ長を使用すると、この問題はさらに悪化する。
 本研究では、StainedGlassを紹介する。このツールは、小さなk-merではなく配列アラインメントに基づいて同一性ヒートマップを生成するもので、簡単で拡張性があり、カスタマイズ可能なワークフローを採用しており、インタラクティブに使用できるだけでなく、出版物用の図を作成することもできます。このツールは、ゲノム全体のリピート構造の研究にも、特定の領域に絞って複雑な高次のリピート構造の特徴を調べるのにも適用できます。最近の8番染色体の解析では、この手法のプロトタイプを開発し、2Mbpセントロメアの高次反復構造をアイデンティティ・ヒートマップとして表示しました。このプロトタイプにより、セントロメアにおける高次の対称性と層状構造の発見が促進され、セントロメアの進化に関するより洗練されたモデルの開発や、コピー数変動のホットスポットの発見に役立ちました(Logsdon et al., 2021)。

 

HP

https://mrvollger.github.io/StainedGlass/

 

 

インストール

Github

mamba create -n snakemake -c conda-forge -c bioconda snakemake -y
conda activate snakemake
git clone https://github.com/mrvollger/StainedGlass.git
cd StainedGlass/

 

 

テストラン

StainedGlass/config/config.yamlを編集する(テストランでは編集は不要)。テストゲノムは.test/という隠しディレクトリに収納されている。実際の解析でStainedGlass/直下にゲノムを置くなら、.test/を消してinput.fastaとすれば認識する。

f:id:kazumaxneo:20210831013305p:plain

パラメータはconfig/README.mdで詳しく書かれています。

 

実行する。

cd StainedGlass/
snakemake --use-conda --cores 24

出力

f:id:kazumaxneo:20210831011447p:plain

results/{sample}.{d+}.{d+}.bedには、パイプラインで同定されたすべてのアラインメントが格納されている。未処理のアラインメントを含むbamファイルも含まれる。複数の染色体配列などを提供した場合は各配列ごとのフォルダができて、その中に出力される。

 

cofigの内容はオプションを立てることで上書きできる。

snakemake --use-conda --cores 24 --config sample=test2 fasta=/some/fasta/path.fa

 

 

 

特定の領域の画像を生成するには、make_figuresを追加する。最大で40Mbpの合計5つの領域を比較するのに適している。

snakemake --use-conda --cores 24 make_figures

 

f:id:kazumaxneo:20210831133540p:plain

T2Tアセンブリのchr8の3Mbの領域を比較。

 

 

ゲノム全体をインタラクティブに可視化するには、HiGlassプログラムとWebブラウザを使う。

pip install higlass-manage #*1
higlass-manage view results/small.5000.10000.strand.mcool

 テストデータ

f:id:kazumaxneo:20210831012228p:plain

シロイヌナズナ

f:id:kazumaxneo:20210831015343p:plain

 

 

高解像度のインタラクティブな視覚化のパラメータ。各ビンにマッピングされたリードの数に比例して色が変わる。計算時間はより多くかかる。

snakemake --use-conda --cores 24 cooler_density --config window=32 cooler_window=100

シロイヌナズナ

f:id:kazumaxneo:20210831020835p:plain

 

 

コマンドの使い方は論文中でも説明されています。確認して下さい。

引用

StainedGlass: Interactive visualization of massive tandem repeat structures with identity heatmaps
Mitchell R. Vollger, Peter Kerpedjiev,Adam M. Phillippy, Evan E. Eichler

bioRxiv, Posted August 21, 2021

 

StainedGlass: Interactive visualization of massive tandem repeat structures with identity heatmaps
Mitchell R Vollger, Peter Kerpedjiev, Adam M Phillippy, Evan E Eichler
Bioinformatics, btac018, https://doi.org/10.1093/bioinformatics/btac018
Published: 10 January 2022

 

HiGlass: web-based visual exploration and analysis of genome interaction maps
Peter Kerpedjiev, Nezar Abdennur, Fritz Lekschas, Chuck McCallum, Kasper Dinkla, Hendrik Strobelt, Jacob M. Luber, Scott B. Ouellette, Alaleh Azhir, Nikhil Kumar, Jeewon Hwang, Soohyun Lee, Burak H. Alver, Hanspeter Pfister, Leonid A. Mirny, Peter J. Park & Nils Gehlenborg
Genome Biology volume 19, Article number: 125 (2018)

 

*1

pip install higlass-pythonも追加で実行した。