Zenn
📘

Strainy : ロングリードメタゲノムリードから株レベルのハプロタイプアセンブリとフェージング

2024/09/28に公開
bioinformatics
ゲノムアセンブリ
metagenome
はじ環
フェージング
tech

NanoporeとHiFiの両方のデータから株レベルのメタゲノムアセンブリとフェージングを行うStrainyを紹介します。

Strainy: phasing and assembly of strain haplotypes from long-read metagenome sequencing

https://www.nature.com/articles/s41592-024-02424-1

微生物群集内の細菌種はごく僅かな変異によって区別される株の混合物と言える。株内の異質性は機能的にも重要であるとされるが、メタゲノム配列からのデータからその特性を明らかにすることは困難とされる。

ショートリードだと、株感の小規模な変異を検出するために使用できるが、これらの変異を連続したハプロタイプのフェーズに分けることは難しい。

ロングリードのメタゲノムアセンブラーは連続したゲノムを生成できる一方、種レベルのコンセンサスを優先して動作することが多い。

ここでは、NanoporeとHiFiリードから株レベルのメタゲノムアセンブリとフェージングを行うStrainyを紹介する。Strainyはde novoメタゲノムアセンブリをInputとして、株の変異を識別してそれらをフェーズ分けして連続したハプロタイプにアセンブリする。

インストール

Githubの手順に従ってconda/mambaでインストールします。

https://github.com/katerinakazantseva/stRainy

mkdir strainy_test && cd strainy_test

# Githubリポジトリをクローン
git clone https://github.com/katerinakazantseva/strainy
cd strainy

git submodule update --init
make -C submodules/Flye

# condaを使ってインストール
mamba env create -f environment.yml -n strainy

cd ..

仮想環境の構築とstrainyのインストールができたら、activateしてhelp表示

mamba activate strainy
./strainy.py -h
strainy USAGE
usage: strainy.py [-h] -o OUTPUT -g GFA -m {hifi,nano} -q FASTQ [-s {phase,transform,e2e}] [--snp SNP] [-t THREADS]
                  [-b BAM] [--link-simplify] [--debug] [--unitig-split-length UNITIG_SPLIT_LENGTH]
                  [--only-split ONLY_SPLIT] [-d CLUSTER_DIVERGENCE] [-a ALLELE_FREQUENCY]
                  [--min-unitig-length MIN_UNITIG_LENGTH] [--min-unitig-coverage MIN_UNITIG_COVERAGE]
                  [--max-unitig-coverage MAX_UNITIG_COVERAGE] [-v]

options:
  -h, --help            show this help message and exit
  -s {phase,transform,e2e}, --stage {phase,transform,e2e}
                        stage to run: either phase, transform or e2e (phase + transform) (default: e2e)
  --snp SNP             path to vcf file with SNP calls to use (default: None)
  -t THREADS, --threads THREADS
                        number of threads to use (default: 4)
  -b BAM, --bam BAM     path to indexed alignment in bam format (default: None)
  --link-simplify       Enable agressive graph simplification (default: False)
  --debug               Generate extra output for debugging (default: False)
  --unitig-split-length UNITIG_SPLIT_LENGTH
                        The length (in kb) which the unitigs that are longer will be split, set 0 to disable
                        (default: 50)
  --only-split ONLY_SPLIT
                        Do not run stRainy, only split long gfa unitigs (default: False)
  -d CLUSTER_DIVERGENCE, --cluster-divergence CLUSTER_DIVERGENCE
                        cluster divergence (default: 0)
  -a ALLELE_FREQUENCY, --allele-frequency ALLELE_FREQUENCY
                        Set allele frequency for internal caller only (pileup) (default: 0.2)
  --min-unitig-length MIN_UNITIG_LENGTH
                        The length (in kb) which the unitigs that are shorter will not be phased (default: 1)
  --min-unitig-coverage MIN_UNITIG_COVERAGE
                        The minimum coverage threshold for phasing unitigs, unitigs with lower coverage will not be
                        phased (default: 20)
  --max-unitig-coverage MAX_UNITIG_COVERAGE
                        The maximum coverage threshold for phasing unitigs, unitigs with higher coverage will not be
                        phased (default: 500)
  -v, --version         show program's version number and exit

Required named arguments:
  -o OUTPUT, --output OUTPUT
                        output directory (default: None)
  -g GFA, --gfa GFA     input gfa to uncollapse (default: None)
  -m {hifi,nano}, --mode {hifi,nano}
                        type of reads (default: None)
  -q FASTQ, --fastq FASTQ
                        fastq file with reads to phase / assemble (default: None)

Tutorialに沿って実行

５つの大腸菌株のモックメタゲノムデータをダウンロード

# wgetが遅いのでaria2cでダウンロード
aria2c -x 15 -ctrue https://zenodo.org/records/11187925/files/strainy_ecoli_example.tar.gz

tar -xvf strainy_ecoli_example.tar.gz

tree -L 1

以下のようなフォルダ構造のはず。

.
├── strainy
├── strainy_ecoli_example
└── strainy_ecoli_example.tar.gz

2 directories, 1 file

Strainy は、PacBio HiFi、Nanopore R9 (Guppy5+)、および R10 シーケンスをサポートしています。

Inputファイルは、metaFlyeから取得可能なGFAと、アセンブル/フェーズ化する対象のfasta/fastqです。提供するfasta/fastqは、もし1つのコンティグに複数のハプロタイプが混ざっているアセンブリ結果を改善したい場合は、同じリードをアセンブラに使用する必要があります。


Direct Phased Genome Assembly Using Nighthawk on HiFi Reads Fig.1

python3 strainy/strainy.py \
--gfa strainy_ecoli_example/ecoli_5strain_metaflye_hap.gfa \
--fastq strainy_ecoli_example/ecoli_5strain_sim_badread.fastq.gz \
--mode nano \
-t 30 \
--output strainy_out

出力ファイル

  • alignment_phased.bam

    • フェーズドアライメントファイル。フェーズはYCタグで定義され、IGVでの視覚化にも使用可能
    • 対応するReference配列は、preprocessing_data/gfa_converted.fasta
    • フェーズなしのアライメントが入力として提供されなかった場合、Strainyはminimap2を使用してgfaに対するアライメントを生成

  • strain_unitigs.gfa

    • 株レベルのハプロタイプアセンブリグラフ。対応する情報はphased_unitig_info_table.csvに記載

  • phased_unitig_info_table.csv

    • 各フェーズの株 unitig の統計 (strain_unitigs.gfaファイルに一致)。各株 unitig について、その長さ、カバレッジ、SNP 率、およびその他の統計が報告されます。

  • strain_contigs.gfa

    • 最終的なアセンブリグラフ。一部の株のユニティグはより長いコンティグに結合されているため、これらの配列はphased_unitig_info_table.csvにとは完全一致しない場合がある。

  • strain_variants.vcf

    • アセンブルされた株のハプロタイプから作成されたバリアント情報ファイル
    • INFOフィールドでは、ALT_HAPはALTバージョンのバリアントをサポートするストレインユニティグを指し、REF_HAPはバリアントを含まないユニティグ(参照状態)のリストに対応

  • multiplicity_stats.txt

    • 重複度やDepth, SNP密度などをまとめたファイル

結果の確認

igvを使って、以下のファイルの視覚化。

  • strainy_out/preprocessing_data/gfa_converted.fasta
  • strainy_out/alignment_phased.bam
  • strainy_out/intermediate/strain_utgs.fasta_ref_aln.bam

図はedge_422のもの。

それっぽく7つにフェージングされているが、全長がカバーできていないフェーズがあり、strain_utgs.fasta_ref_aln.bamでは4つにクラスタリングされてハプロタイプが区別されている。

また、turorialのようにbandageを使ってアセンブリグラフデータの視覚化する。

mamba install -c bioconda -c conda-forge bandage -y

起動

Bandage &

strain_unitigs.gfaをインポートして、IGVで確認していたedge_422を表示

strain_utgs.fasta_ref_aln.bamと同様に４つのハプロタイプを含むunitigが作成されていました。

edge422について、unitigの情報をphased_unitig_info_table.csvで確認することも可能。

Strain_unitig Reference_unitig Length Coverage Abundance_Ratio #SNP SNP_density Start_positioin End_position
edge_422_20109 edge_422 8541 15 13 108 0.0126 0 8548
edge_422_20111 edge_422 8548 29 25 110 0.0128 0 8548
edge_422_30111 edge_422 8548 23 20 110 0.0128 0 8548
edge_422_2010004 edge_422 8548 26 23 110 0.0128 0 8548

おわり。

Discussion