nf-coreのDeepVariantパイプライン

nf-core/deepvariantより

2017年12月にGoogleブレインチームがDeepLearningをベースにしたVariant Caller, DeepVariantをリリースした。DeepVariantはまずBAMファイルに基づいて画像を構築し、次にDeepLearningの画像認識アプローチを使用してバリアントを取得し、最終的には予測の出力を標準的なVCF形式に変換する。

NextflowパイプラインとしてのDeepVariantは、ユーザーにいくつかの利点を提供する。これは、DeepVariantに必要な入力であり、通常はユーザーが手動で作成する必要がある、いくつかの余分な必要なインデックス付き圧縮ファイルの作成を、前処理ステップを介して自動的に処理する。Variant Callingは、複数のBAMファイルに対して同時に実行でき、Nextflowの内部並列化のおかげで、リソースを無駄にすることがない。NextflowはDockerをサポートしているため、Dockerコンテナ内ですべてのステップを実行することで、計算上の再現性とクリーンな方法で結果を生成することができる。

DeepVariant: Highly Accurate Genomes With Deep Neural Networks

manual

https://nf-co.re/deepvariant/1.0/docs/usage

インストール

Github

# Make sure that Java v8+ is installed:
java -version
# Install Nextflow (持ってない人だけ、condaでも導入可能)
curl -fsSL get.nextflow.io | bash
mv nextflow ~/bin/

#test run (ユーザーにdocker実行権限がなければsudo実行する)
nextflow run nf-core/deepvariant -profile test,docker

出力

実行方法

アセンブリのバージョンとbam、ターゲットのbedファイルを指定する。複数のbamがある場合は--bam_folderオプションを使い、単体のbamでは--bamオプションを使う。

#WGS
nextflow run nf-core/deepvariant --genome hg19 --bam yourBamFile --bed yourBedFile -profile standard,docker

#Exome
nextflow run nf-core/deepvariant --exome --genome hg19 --bam_folder myBamFolder --bed myBedFile -profile standard,docker

• --bed   Path to bedfile, specifying region to be analysed must also be supplied
• --bam_folder   Use this to specify a folder containing BAM files. Allows multiple BAM files to be analyzed at once. All BAM files will be analyzed unless --bame_file_prefix is used (see below). For example:
• --bam   Use this to specify the BAM file
• --genome   Standard versions of the genome are prepared with all their compressed and indexed file in a lifebit s3 bucket (hg19| h38 | grch37primary | hs37d5 | hg19chr20).
• --exome    For exome bam files
• --max_memory    Use to set a top-limit for the default memory requirement for each process. Should be a string in the format integer-unit. eg. `--max_memory '8.GB'``

• --max_time    Use to set a top-limit for the default time requirement for each process. Should be a string in the format integer-unit. eg. --max_time '2.h'

• --max_cpus    Use to set a top-limit for the default CPU requirement for each process. Should be a string in the format integer-unit. eg. --max_cpus 1

引用

https://nf-co.re/deepvariant

A universal SNP and small-indel variant caller using deep neural networks
Ryan Poplin, Pi-Chuan Chang, David Alexander, Scott Schwartz, Thomas Colthurst, Alexander Ku, Dan Newburger, Jojo Dijamco, Nam Nguyen, Pegah T Afshar, Sam S Gross, Lizzie Dorfman, Cory Y McLean, Mark A DePristo

Nat Biotechnol. 2018 Nov;36(10):983-987

参考

生殖細胞バリアントや体細胞バリアントを検出する自動化されたパイプライン Sarek

2021 2/9 タイトル修正

2021 2/12, 2/15コマンド追記

2021 5/16 コメント追記

　全ゲノムシークエンシング(WGS)は、精密医療の発展のための研究の基盤技術であるが、WGS解析のためのポータブルで使いやすいワークフローが限られていることが、多くの研究グループにとって大きな課題となっており、科学の進歩を阻害する要因となっている。ここでは、WGS、全ゲノムシークエンシング（WES）、遺伝子パネルからのシークエンシングデータに基づいて、生殖細胞変異や体細胞突然変異を検出するためのオープンソースのワークフローであるSarekを紹介する。Sarekの特徴は、(i)簡単なインストール、(ii)異なるコンピュータ環境での堅牢な移植性、(iii)包括的なドキュメント、(iv)透明で読みやすいコード、(v)広範な品質指標の報告である。SarekはNextflowワークフロー言語で実装されており、DockerやSingularityコンテナ、Conda環境の両方をサポートしているため、POSIX互換のコンピュータやクラウドコンピュート環境に簡単に導入することができる。Sarekは、リードアライメントと前処理に関するGATKのベスト・プラクティスに準拠しており、生殖細胞および体細胞のシングルヌクレオチドバリアント、挿入・欠失バリアント、構造バリアント、腫瘍サンプルの純度、およびploidyとコピー数のバリエーションの同定とアノテーションのための幅広いソフトウェアが含まれている。Sarekは、簡単、効率的、再現性のあるWGS解析を提供し、シーケンシング施設での生産ワークフローとして、また、個々の研究グループのための強力なスタンドアロンツールとして、容易に使用することができる。Sarek のソースコード、ドキュメント、インストール手順書は https://github.com/nf-core/sarek および https://nf-co.re/sarek/ で自由に入手できる。

usage

https://nf-co.re/sarek/usage#introduction

インストール

ubuntu18.04LTSでテストした。

Github

# Make sure that Java v8+ is installed:
java -version
# Install Nextflow (持ってない人だけ、condaでも導入可能)
curl -fsSL get.nextflow.io | bash
mv nextflow ~/bin/

#test run (権限がないならsudo実行、数時間かかる。docker/singularity/podman、に対応している。condaも対応しているが推奨されていない）
#dockerを使用
nextflow run nf-core/sarek -profile test,docker

result

Preprocessing/

Reports/

Reports/9876T/

VariantCalling/

pipeline_info/

実行方法

実際にサンプルをランする際には、fastq（または処理済みのbam）のパス（--input）、プロファイル（-profile）、サンプルグループを記載したTSVファイルを指定する。

nextflow run nf-core/sarek --input sample.tsv --step mapping --use_gatk_spark --save_bam_mapped --generate_gvcf -profile docker

• --input　 Path to input file(s).
• --outdir　The output directory where the results will be saved.
• --step   Starting step. default:'mapping'
• --cpus　default:8

• --use_gatk_spark   Enable usage of GATK Spark implementation

• --save_bam_mapped   Save Mapped BAMs

• --generate_gvcf   Generate g.vcf output from GATK HaplotypeCaller

• --single_cpu_mem   Use to set memory for a single CPU.　default:'7 GB'
• --max_memory   Maximum amount of memory that can be requested for any single job.　default:'128.GB'
• -profile  Use this parameter to choose a configuration profile. Profiles can give configuration presets for different compute environments （docker|singularity|podman|conda|test|test_annotation|test_use_gatk_spark|

test_split_fastq|test_targeted|test_tool|test_trimming|test_umi_qiaseq|test_umi_tso）.

サンプルグループを記載するTSVファイルについてはmanualを確認して下さい。 

 メモ

メモリが足りない可能性があるので、多く利用できるなら増やしておく(option)。

Sentieon (commercial solution)

Germline SNV/INDEL Variant Calling - DNAseq

nextflow run nf-core/sarek --input sample.tsv -profile docker --sentieon and --tools DNAseq 

•  --sentieon   Enable Sentieon if available.

• --tools   Tools to use for variant calling and/or for annotation.

• --step  Starting step. default:'mapping'

メモリは多めに使います。最低128GB程度は必要のようです。

2021 2/16

バリアントコールまで行うには商用のオプションが必須のようです。勘違いしておりました。

引用

Sarek: A portable workflow for whole-genome sequencing analysis of germline and somatic variants

 Maxime Garcia, Szilveszter Juhos, Malin Larsson, Pall I. Olason,
Marcel Martin, Jesper Eisfeldt, Sebastian DiLorenzo, Johanna Sandgren,
Teresita Díaz De Ståhl, Philip Ewels, Valtteri Wirta, Monica Nistér,
Max Käller, Björn Nystedt

 F1000Res, Version 2. doi: 10.12688/f1000research.16665.2

配列をアセンブリグラフにマッピングしてグラフを拡張する minigraph

最近のシーケンシング技術の進歩により、個々のゲノムを参照ゲノムの質に合わせて組み立てることが可能になった。同一種からの複数のゲノムを統合し、統合された表現を生物学者が利用できるようにするにはどうすればよいのかは、依然として未解決の課題である。ここでは、直線的なリファレンスゲノムの座標を維持しながら複数のゲノムを表現するためのグラフベースのデータモデルと関連するフォーマットを提案する。このアイデアをminigraphツールキットに実装し、パンゲノムグラフを効率的に構築し、現在のリファレンスゲノムから欠落している何万もの構造変異をコンパクトにエンコードできることを実証した。

Github

Minigraphはシーケンスからグラフへのマッパーであり、グラフのコンストラクタである。Minigraphは，シーケンスグラフ内の問い合わせシーケンスの近似位置を見つけ，既存のグラフから分岐した長い問い合わせ部分シーケンスで既存のグラフを増分的に拡張する。Minigraphはminimap2から多くのアイデアとコードを借りている。かなり効率的で、24CPUコアを使って40のヒトアセンブリから半日でグラフを構築することができる。ベースアライメントがないこともあり、minigraphは最適ではないマッピングやローカルグラフを生成することがあり得る。

2022/05/11

Minigraph v0.18 released with improved performance for species of higher diversity. https://t.co/7fHB8CtXCf

— Heng Li (@lh3lh3) 2022年5月9日

インストール

Github

git clone https://github.com/lh3/minigraph
cd minigraph && make

#gfatools
git clone https://github.com/lh3/gfatools
cd gfatools && make

> ./minigraph

$ ./minigraph

Usage: minigraph [options] <target.gfa> <query.fa> [...]

Options:

Indexing:

-k INT k-mer size (no larger than 28) [15]

-w INT minizer window size [10]

Mapping:

-f FLOAT ignore top FLOAT fraction of repetitive minimizers [0.0002]

-U INT[,INT] choose the minimizer occurrence threshold within this interval [50,250]

-j FLOAT expected sequence divergence [0.1]

-g NUM stop chain enlongation if there are no minimizers in INT-bp [10000]

-F NUM max fragment length (effective with -xsr or in the fragment mode) [800]

-r NUM bandwidth used in chaining [2000]

-n INT[,INT] minimal number of minimizers on a graph/linear chain [3,2]

-m INT[,INT] minimal graph/linear chaining score [50,30]

-p FLOAT min secondary-to-primary score ratio [0.8]

-N INT retain at most INT secondary mappings [5]

-D skip self diagonal matches

Graph generation:

--ggen perform incremental graph generation

-q INT min mapping quality [5]

-l NUM min alignment length [100000]

-d NUM min alignment length for depth calculation [10000]

-L INT min variant length [50]

--call call the graph path in each bubble and output BED

Input/output:

-t INT number of threads [4]

-o FILE output mappings to FILE [stdout]

-K NUM minibatch size for mapping [500M]

-S output linear chains in * sName sLen nMz div sStart sEnd qStart qEnd

--vc output in the vertex coordinate

Preset:

-x STR preset []

- lr: noisy long read mapping (the default)

- asm: asm-to-ref mapping

- sr: short reads

- ggs: simple algorithm for graph generation

 > ./gfatools

$ ./gfatools 

Usage: gfatools <command> <arguments>

Commands:

  view        read a GFA file

  stat        statistics about a GFA file

  gfa2fa      convert GFA to FASTA

  gfa2bed     convert rGFA to BED (requiring rGFA)

  blacklist   blacklist regions

  bubble      print bubble-like regions (EXPERIMENTAL)

  asm         miniasm-like graph transformation

  sql         export rGFA to SQLite (requiring rGFA)

  version     print version number

実行方法

GFAフォーマットのリファレンスにサンプルのFASTA配列をマッピングして、グラフを拡張する。配列をグラフにマッピングするには、GFA フォーマット（ 推奨はrGFA フォーマット）のリファレンスグラフファイルを用意する必要がある。

minigraph -x lr graph.gfa query.fa > out.gaf

グラフがない場合は、複数のサンプルの配列からグラフを生成することもできる。

#リファレンスに１サンプルの配列をマッピング
minigraph -x ggs -t16 ref.fa sample1.fa > sample1.gfa

#リファレンスに2サンプルの配列をマッピング
minigraph -x ggs -t16 ref.fa sample1.fa sample2.fa > out.gfa

#GFAのグラフに配列をマッピング
minigraph -x ggs -t16 ref.gfa sample.fa > out.gfa

minigraph GFAパーサーはFASTAをシームレスに解析するため、FASTA中の配列をリファレンスとして提供することができる。この場合、minigraphはminimap2のように動作するが、ベースレベルのアラインメントではないと説明されている。

テスト

2つの配列を用意して挙動を確認する。リファレンスの配列、その中央に６kbくらいのリファレンスにない配列を挿入した配列（変異リファレンス）の２つを用意した。

minigraphを実行して、その２つの配列からグラフを出力。

minigraph -xggs -t 12 -l 1k ref.fa mt.fa > output.gfa

• -l   min alignment length [100000]

Bandageで視覚化した。想定通り、中央にバブルが発生している。

GFAファイル(link)

最初の3行が配列情報（S）、次の3行がリンク情報（L）。s1が正方向でそれぞれs2の正方向とs3の正方向とギャップなしで繋がる。s2の正方向がs3の正方向とギャップなしで繋がる。s1とs3についてもリンクの記述があれば、環状化して閉じた環になる。今回は線状配列を提供したので記載はない。

*リファレンスよりクエリの方が長い配列だとしても、完全にクエリにしかない領域はリファレンスンにアラインメントされず、出力されるグラフには表現されない。リファレンスに部分的にはアラインメントされ、途中で構造が違う領域だけがグラフに反映される。

FASTAに変換(lossy)

gfatools gfa2fa -s output.gfa > output.fa

SVコール

minigraph -xggs --call -l 1k ref.fa mt.fa > call.bed

出力はBEDライクなフォーマットになる。

引用

The design and construction of reference pangenome graphs with minigraph
Heng Li, Xiaowen Feng & Chong Chu
Genome Biology volume 21, Article number: 265 (2020)

アセンブリグラフを用いたメタゲノムコンティグのビニングを行う GraphBin2

　メタゲノムシークエンシングは、微生物群集の構造、多様性、生態を純粋な培養物を得ることなく研究することを可能にする。多くのメタゲノム研究では、メタゲノムシークエンシングから得られたリードは、最初に長いコンティグにアセンブリされ、これらのコンティグは、クラスタ内のコンティグが同じ種に由来すると予想されるコンティグのクラスタにビン分けされる。異なる種がゲノム中で共通の配列を共有している場合があるため、1つのコンティグが複数の種に属している可能性がある。しかし、既存のツールでは、コンティグのビン分けは非オーバーラップビン分けしかサポートしておらず、各コンティグは最大で1つのビン（種）に割り当てられている。本論文では、既存のツールで得られたビニング結果を改良し、コンティグを複数のビンに割り当てることができる GraphBin2 を紹介する。GraphBin2は、アセンブリグラフの接続性とカバレッジ情報を用いて、既存のコンティグのビニング結果を調整し、複数の種が共有するコンティグを推定する。シミュレーションデータと実データを用いた実験結果から、 GraphBin2は既存のツールのビニング結果を改善するだけでなく、コンティグを複数のビンに割り当てることも可能であることが示された。

GraphBin

インストール

付属のyamlファイルでcondaの仮想環境を作成してテストした（ubuntu18.04LTS）。

Github

git clone https://github.com/Vini2/GraphBin2.git
cd GraphBin2/
conda env create -f environment.yml
conda activate graphbin2

$ ./graphbin2 -h

¥usage: graphbin2 [-h] --assembler ASSEMBLER --graph GRAPH --contigs CONTIGS [--paths PATHS] [--abundance ABUNDANCE] --binned BINNED --output OUTPUT [--prefix PREFIX] [--depth DEPTH] [--threshold THRESHOLD] [--delimiter DELIMITER] [--nthreads NTHREADS]

GraphBin2 Help. GraphBin2 is a tool which refines the binning results obtained from existing tools and, more importantly, is able to assign contigs to multiple bins. GraphBin2 uses the connectivity and coverage information from assembly graphs to adjust existing binning

results on contigs and to infer contigs shared by multiple species.

optional arguments:

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

--assembler ASSEMBLER

name of the assembler used (SPAdes, SGA or Flye)

--graph GRAPH path to the assembly graph file

--contigs CONTIGS path to the contigs file

--paths PATHS path to the contigs.paths file

--abundance ABUNDANCE

path to the abundance file

--binned BINNED path to the .csv file with the initial binning output from an existing tool

--output OUTPUT path to the output folder

--prefix PREFIX prefix for the output file

--depth DEPTH maximum depth for the breadth-first-search. [default: 5]

--threshold THRESHOLD

threshold for determining inconsistent vertices. [default: 1.5]

--delimiter DELIMITER

delimiter for input/output results. Supports a comma (,), a semicolon (;), a tab ($'\t'), a space (" ") and a pipe (|) [default: , (comma)]

実行方法

１、メタゲノムアセンブリとビニングを実行する（Github参照）。

２、graphbin2のラン。アセンブラ（metaspades、SGA, metaFlye）とコンティグ、グラフファイル、パスファイル、ビニングレポートファイルを指定する。

python graphbin2.py --assembler spades --contigs contigs.fasta --graph graph.gfa --paths paths_file.paths --binned binning_result.csv --output output_dir

引用

GraphBin2: Refined and Overlapped Binning of Metagenomic Contigs Using Assembly Graphs
 

Mallawaarachchi, Vijini G. ; Wickramarachchi, Anuradha S. ; Lin, Yu

DOI: 10.4230/LIPIcs.WABI.2020.8

20th International Workshop on Algorithms in Bioinformatics (WABI 2020)

ロングリードからトランスポゾンを検出する TELR

　TELR（Tellerと発音）は、ロングリードシーケンシングデータ（PacBioまたはOxford Nanopore）からの高速な非リファレンストランスポーザブルエレメント（TE）検出器である。TELRは、リファレンスゲノムにマッピングされたロングリードを使用してSnifflesを使用して挿入を識別し、挿入を支持するリードとユーザーが提供したTEコンセンサス配列を照合することで挿入をフィルタリングする。各TE挿入候補遺伝子座について、TELRはTE挿入をサポートするすべてのリードのローカルアセンブリを行い、アセンブリーされたコンティグにTE配列をアノテーションした後、そのフランクをリファレンスゲノムにマップする。最後に、TELRは非参照TE挿入の座標とアセンブルされたTE配列を生成する。

　現在のバージョンのTELRは、ヘテロ接合性のTE挿入を含む実際のDrosophila melanogasterのデータセットで良好なパフォーマンスを示している。

インストール

付属のyamlファイルでcondaの仮想環境を作成してテストした（ubuntu18.04LTS）。

Github

git clone git@github.com:bergmanlab/TELR.git
cd TELR
conda env create -f envs/telr.yml
conda activate TELR_env

> python telr.py -h

$ python telr.py -h

usage: telr.py [-h] -i READS -r REFERENCE -l LIBRARY [-x PRESETS] [-p POLISH]

               [-o OUT] [-t THREAD] [-g GAP] [-v OVERLAP] [-k]

Script to detect TEs in long read data

required arguments:

  -i READS, --reads READS

                        reads in fasta/fastq format or read alignments in bam

                        format

  -r REFERENCE, --reference REFERENCE

                        reference genome in fasta format

  -l LIBRARY, --library LIBRARY

                        TE consensus sequences in fasta format

optional arguments:

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

  -x PRESETS, --presets PRESETS

                        parameter presets for different sequencing

                        technologies (default = 'pacbio')

  -p POLISH, --polish POLISH

                        rounds of contig polishing (default = 1)

  -o OUT, --out OUT     directory to output data (default = '.')

  -t THREAD, --thread THREAD

                        max cpu threads to use (default = '1')

  -g GAP, --gap GAP     max gap size for flanking sequence alignment (default

                        = '20')

  -v OVERLAP, --overlap OVERLAP

                        max overlap size for flanking sequence alignment

                        (default = '20')

  -k, --keep_files      If provided then all intermediate files will be kept

                        (default: remove intermediate files)

実行方法

ロングリードとリファレンスのfasta、TEのコンセンサス配列を指定する。

python3 telr.py -i long_read.fq -r ref.fasta -t 20 -x pacbio -l test/library.fasta

• -i    reads in fasta/fastq format or read alignments in bam format

• -r    reference genome in fasta format

• -l    TE consensus sequences in fasta format

• -x    parameter presets for different sequencing technologies (default = 'pacbio')

• -o    directory to output data (default = '.')

• -t     max cpu threads to use (default = '1')

Outputについて

https://github.com/bergmanlab/TELR/blob/master/docs/03_Output_Files.md

引用

https://github.com/bergmanlab/TELR

2020 Shunhua Han and Casey M. Bergman

RNAのロングリードを分析する IsoQuant

　IsoQuantは、PacBioやOxford Nanoporesのような長いRNAリードのリファレンスベース解析のためのツールである。IsoQuantは、リファレンスゲノムにリードをマッピングし、それらのイントロンとエクソンの構造に基づいて、アノテーションされたアイソフォームに割り当てる。IsoQuantはまた、イントロンの保持、代替スプライスサイト、スキップされたエクソンなどの様々な修飾を発見することができる。IsoQuant はさらに、遺伝子、アイソフォーム、エクソン、イントロンの定量を行う。リードがグループ化されている場合（例えば、細胞タイプに応じて）、カウントは提供されたグループ化に従って報告される。さらに、IsoQuant は新規のものも含めて、発見された転写モデルを生成する。IsoQuantバージョン1.1は、2020年12月11日にGPLv2の下でリリースされ、https://github.com/ablab/IsoQuant からダウンロードできる。

HP

http://cab.spbu.ru/software/isoquant/

Meet IsoQuant - a new tool for long-read RNA analysis! It can detect structural alternations in individual reads, quantify isoforms and construct novel transcript models. Works well on PacBio and ONT reads, and has features for single-cell data analysis. https://t.co/KzYpajdq74

— Center for Algorithmic Biotechnology (@bioinf_spbu) 2020年7月14日

インストール

condaを使ってpython3.7の仮想環境を作ってテストした（OSはubuntu18.04LTS）。

依存
IsoQuant requires a 64-bit Linux system or Mac OS and Python (3.7 and higher) to be pre-installed on it. You will also need

• gffutils
• pysam
• biopython
• pyfaidx
• pandas
• numpy
• minimap2
• samtools)
• STAR (optional)

Github

#bioconda (link)
conda create -n isoquant -y python=3.7
conda activate isoquant
conda install -c bioconda isoquant -y
#or 
conda install -c bioconda -c conda-forge -c isoquant isoquant -y

#gffutilsｙとpybedtoolsが導入されなかったので、追加で入れた
pip install gffutils pybedtools

#from source
git clone https://github.com/ablab/IsoQuant.git
cd IsoQuant
git checkout latest
pip install -r requirements.txt
#=> You also need minimap2 to be in the $PATH variable.

> isoquant.py -h

$ isoquant.py -h

usage: isoquant.py [-h] [--output OUTPUT] --genedb GENEDB [--complete_genedb]

                   [--reference REFERENCE] [--index INDEX]

                   (--bam BAM [BAM ...] | --fastq FASTQ [FASTQ ...] | --bam_list BAM_LIST | --fastq_list FASTQ_LIST)

                   --data_type {assembly,pacbio_raw,pacbio_ccs,nanopore}

                   [--stranded STRANDED] [--has_polya] [--fl_data]

                   [--full_help] [--test] [--threads THREADS]

                   [--labels LABELS [LABELS ...]] [--read_group READ_GROUP]

                   [--sqanti_output] [--count_exons]

                   [--matching_strategy {exact,precise,default,loose}]

                   [--model_construction_strategy {reliable,default,fl,all,assembly}]

optional arguments:

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

  --output OUTPUT, -o OUTPUT

                        output folder, will be created automatically

                        [default=isoquant_output]

  --genedb GENEDB, -g GENEDB

                        gene database in gffutils DB format or GTF/GFF format

  --complete_genedb     use this flag if gene annotation contains transcript

                        and gene metafeatures, e.g. with official annotations,

                        such as GENCODE; speeds up gene database conversion

  --reference REFERENCE, -r REFERENCE

                        reference genome in FASTA format, should be provided

                        to compute some additional stats and when raw reads

                        are used as an input

  --index INDEX         genome index for specified aligner, should be provided

                        only when raw reads are used as an input

  --bam BAM [BAM ...]   sorted and indexed BAM file(s), each file will be

                        treated as a separate sample

  --fastq FASTQ [FASTQ ...]

                        input FASTQ file(s), each file will be treated as a

                        separate sample; reference genome should be provided

                        when using raw reads

  --bam_list BAM_LIST   text file with list of BAM files, one file per line,

                        leave empty line between samples

  --fastq_list FASTQ_LIST

                        text file with list of FASTQ files, one file per line,

                        leave empty line between samples

  --data_type {assembly,pacbio_raw,pacbio_ccs,nanopore}, -d {assembly,pacbio_raw,pacbio_ccs,nanopore}

                        type of data to process, supported types are:

                        assembly, pacbio_raw, pacbio_ccs, nanopore

  --stranded STRANDED   reads strandness type, supported values are: forward,

                        reverse, none

  --has_polya           set if reads were not polyA trimmed; polyA tails will

                        be detected and further required for transcript model

                        construction

  --fl_data             reads represent FL transcripts; both ends of the read

                        are considered to be reliable

  --full_help           show full list of options

  --test                run IsoQuant on toy dataset

  --threads THREADS, -t THREADS

                        number of threads to use

  --labels LABELS [LABELS ...], -l LABELS [LABELS ...]

                        sample names to be used; input file names are used if

                        not set

  --read_group READ_GROUP

                        a way to group feature counts (no grouping by

                        default): by BAM file tag (tag:TAG), using additional

                        file (file:FILE:READ_COL:GROUP_COL:DELIM), using read

                        id (read_id:DELIM)

  --sqanti_output       produce SQANTI-like TSV output (requires more time)

  --count_exons         perform exon and intron counting

  --matching_strategy {exact,precise,default,loose}

                        matching strategy to use from most strict to least

  --model_construction_strategy {reliable,default,fl,all,assembly}

                        transcritp model construnction strategy to use

インストールチェック

> isoquant.py --test

$ isoquant.py --test

=== Running in test mode === 

Any other option is ignored 

2020-08-18 00:37:15,601 - INFO -  === IsoQuant pipeline started === 

2020-08-18 00:37:15,601 - INFO - Indexing reference

2020-08-18 00:37:15,633 - INFO - Aligning /home/kazu/anaconda3/share/isoquant-1.0.0-0/tests/toy_data/MAPT.Mouse.ONT.simulated.fastq to the reference

2020-08-18 00:37:17,727 - INFO - Converting gene annotation file to .db format (takes a while)...

2020-08-18 00:37:18,263 - INFO - Gene database written to /home/kazu/Documents/isoquant_test/MAPT.Mouse.genedb.db

2020-08-18 00:37:18,264 - INFO - Provide this database next time to avoid excessive conversion

2020-08-18 00:37:18,264 - INFO - Loading gene database from /home/kazu/Documents/isoquant_test/MAPT.Mouse.genedb.db

2020-08-18 00:37:18,266 - INFO - Loading reference genome from /home/kazu/anaconda3/share/isoquant-1.0.0-0/tests/toy_data/MAPT.Mouse.reference.fasta

2020-08-18 00:37:18,283 - INFO - Processing 1 sample

2020-08-18 00:37:18,283 - INFO - Processing sample 00_MAPT.Mouse.ONT.simulated

2020-08-18 00:37:18,283 - INFO - Sample has 1 BAM file: isoquant_test/00_MAPT.Mouse.ONT.simulated/00_MAPT.Mouse.ONT.simulated.bam

2020-08-18 00:37:18,283 - INFO - Processing chromosome chr11

2020-08-18 00:37:18,367 - INFO - Combining output

2020-08-18 00:37:18,373 - INFO - Finished processing sample 00_MAPT.Mouse.ONT.simulated

2020-08-18 00:37:18,373 - INFO - Gene counts are stored in isoquant_test/00_MAPT.Mouse.ONT.simulated/00_MAPT.Mouse.ONT.simulated.gene_counts.tsv

2020-08-18 00:37:18,373 - INFO - Transcript counts are stored in isoquant_test/00_MAPT.Mouse.ONT.simulated/00_MAPT.Mouse.ONT.simulated.transcript_counts.tsv

2020-08-18 00:37:18,373 - INFO - Transcript model file isoquant_test/00_MAPT.Mouse.ONT.simulated/transcript_models.gff

2020-08-18 00:37:18,373 - INFO - Processed sample 00_MAPT.Mouse.ONT.simulated

2020-08-18 00:37:18,373 - INFO - Read assignment statistics

2020-08-18 00:37:18,373 - INFO - empty: 15

2020-08-18 00:37:18,373 - INFO - unique: 117

2020-08-18 00:37:18,374 - INFO - Transcript model statistics

2020-08-18 00:37:18,374 - INFO - known: 10

2020-08-18 00:37:18,374 - INFO - Processed 1 sample

2020-08-18 00:37:18,374 - INFO -  === IsoQuant pipeline finished === 

2020-08-18 00:37:18,375 - INFO -  === TEST PASSED CORRECTLY === 

O.K 

実行方法

ナノポアdRNAリードを指定。ポリAはトリミングされていない(--has_polya)。ゲノムとGTF(--genedb)のオフィシャル(--complete_genedb)アノテーションを指定、ファイル名の代わりにサンプルラベルを使用。サンプル名 My_ONT。

isoquant.py --data_typ nanopore --has_polya --stranded forward --fastq ONT.raw.fastq.gz --reference reference.fasta --genedb annotation.gtf --complete_genedb --output output_dir --threads 12 --labels My_ONT

• --data_type {assembly,pacbio_raw,pacbio_ccs,nanopore} 　 type of data to process, supported types are:  assembly, pacbio_raw, pacbio_ccs, nanopore

• --has_polya　 set if reads were not polyA trimmed; polyA tails will be detected and further required for transcript model construction

• --stranded 　reads strandness type, supported values are: forward, reverse, none

• --fastq   input FASTQ file(s), each file will be treated as a separate sample; reference genome should be provided when using raw reads
•  --reference    reference genome in FASTA format, should be provided to compute some additional stats and when raw reads are used as an input
• --genedb 　 gene database in gffutils DB format or GTF/GFF format
• --complete_genedb     use this flag if gene annotation contains transcript and gene metafeatures, e.g. with official annotations, such as GENCODE; speeds up gene database conversion

• --output    output folder, will be created automatically [default=isoquant_output]

• --threads    number of threads to use

• --labels    sample names to be used; input file names are used if not set

出力ディレクトリにサンプルごとのサブフォルダができる。

引用

https://github.com/ablab/IsoQuant

関連

ハイブリッドアセンブリとビニング及び下流解析を行う自動化されたパイプライン MUFFIN

2021 2/5,2/6 出力例追記

2021 2/11 論文引用

メタゲノミクスは微生物学の多くの分野を再定義した。しかし、メタゲノムアセンブルゲノム（MAG）は、主にショートリードでシーケンスが行われた場合、断片化されていることが多い。最近のロングリードシーケンシング技術は、ゲノム再構成を改善することを約束している。しかし、2つの異なるシークエンシング技術を統合すると、下流の解析が複雑になる。そこで著者らは、ショートリードとロングリードを用いて高品質のビンとそのアノテーションを生成する完全なメタゲノムワークフローであるMUFFINを開発した。このワークフローは、ワークフローオーケストレーションソフトウェアであるNextflowを用いて記述されており、高い再現性と高速かつ簡単な使用を実現している。このワークフローはまた、ビンの分類学的分類とKEGGパスウェイを生成し、定量化とアノテーションのためのRNA-Seqデータ（オプション）を提供することができる。20 個のバイオガスリアクターサンプルを用いてワークフローをテストし、微生物群集とその機能を解析するために必要なファイルを処理して出力する MUFFIN の能力を評価した。MUFFINは機能パスウェイの予測を行い、提供された場合には、メタゲノムサンプルと各ビンのデノボ転写産物のアノテーションを生成する。

The Workflow. Githubより転載

公開当初は一部プロセスが開発中でランできなかったが、その後バージョンアップによって全プロセスが通しでランできるようになった（らしい）。

インストール

依存

• nextflow
• conda

本体　Github

nextflow run RVanDamme/MUFFIN --parameters

> nextflow run RVanDamme/MUFFIN --parameter -profile local,docker

$ nextflow run RVanDamme/MUFFIN --parameter -profile local,docker

N E X T F L O W  ~  version 20.07.1

Launching `RVanDamme/MUFFIN` [extravagant_cuvier] - revision:

347367e00d [master]

[-        ] process > sourmash_download_db   -

[-        ] process > sourmash_download_db   [  0%] 0 of 1

[-        ] process > discard_short          -

[-        ] process > merge                  -

[-        ] process > fastp                  -

[-        ] process > spades                 -

[-        ] process > minimap2               -

[-        ] process > bwa                    -

[-        ] process > metabat2               -

[-        ] process > maxbin2                -

[-        ] process > concoct                -

[-        ] process > refine3                -

[-        ] process > checkm                 -

[-        ] process > sourmash_bins          -

[-        ] process > sourmash_checkm_parser -

[-        ] process > eggnog_download_db     [  0%] 0 of 1

[-        ] process > eggnog_bin             -

[-        ] process > parser_bin             -

[-        ] process > readme_output          [  0%] 0 of 1

*********Start running MUFFIN*********

MUFFIN is a hybrid assembly and differential binning workflow for

metagenomics, transcriptomics and pathway analysis.

If you use MUFFIN for your research pleace cite:

https://www.biorxiv.org/content/10.1101/2020.02.08.939843v1

or

Van Damme R., Hölzer M., Viehweger A., Müller B., Bongcam-Rudloff E.,

Brandt C., 2020

"Metagenomics workflow for hybrid assembly, differential coverage

binning, transcriptomics and pathway analysis (MUFFIN)",

doi: https://doi.org/10.1101/2020.02.08.939843

**************************************

No files match pattern `*_R{1,2}.fastq{,.gz}` at path: ./illumina/

No files match pattern `*.fastq{,.gz}` at path: ./nanopore/

テストラン

#test locally with conda
nextflow run RVanDamme/MUFFIN --output results_dir --cpus 8 --memory 32g -profile local,conda,test

#test locally with docker (docker実行権限がユーザに無いならsudo実行)
nextflow run RVanDamme/MUFFIN --output results_dir --rna -profile local,docker,test

実行するとdockerイメージのダウンロードやcondaの環境が作られ、ジョブが進行してい（データベースのダウンロードで100GB以上のディスクスペースを使うようなのでコマンド実行時は注意）。

テストランには半日かかった。

出力

assemble

classify

classify_step_summary.csv

annotate

MUFFIN_sample_result.html

実際にfastqを指定する時は--ontと--illuminaフラグを立て、fastqのディレクトリを指定する。アセンブラはmetaspadesとmetaflyeから選ぶ。

nextflow run RVanDamme/MUFFIN --output result --ont nanopore/ --illumina illumina/ --assembler metaspades --rna rna/ -profile local,docker --modular full

fastqは解凍し、イルミナはxxx_R1.fastqとxxx_R2.fastq、ONTはxxx.fastqのファイル命名則に従う必要がある。

引用

Metagenomics workflow for hybrid assembly, differential coverage binning, transcriptomics and pathway analysis (MUFFIN)
Renaud Van Damme, Martin Hölzer,Adrian Viehweger,Bettina Müller, Erik Bongcam-Rudloff, Christian Brandt

BioRxiv, Posted February 08, 2020

2021 2/11

Metagenomics workflow for hybrid assembly, differential coverage binning, metatranscriptomics and pathway analysis (MUFFIN)
Renaud Van Damme , Martin Hölzer, Adrian Viehweger, Bettina Müller, Erik Bongcam-Rudloff

PLOS Computational Biology, Published: February 9, 2021

（2/11現在は原稿のページ校正が反映されていない早期版）

関連