同種の微生物細胞が共存すると、栄養嗜好から病原性までの表現型に影響を与える遺伝的変異を示すことが多い。本発表では、メタゲノムのペアエンドリードを用いて、全ゲノムにわたる集団内の遺伝的多様性（マイクロダイバーシティ）をプロファイリングし、マイクロダイバーシティを考慮した方法で微生物集団を比較するプログラムinStrainを紹介する。inStrainは既存の方法とベンチマークした場合のゲノム比較の精度を大幅に向上させる。inStrainを用いて、新生児未熟児の糞便メタゲノム1,000個以上をプロファイリングし、一卵性双生児は二卵性双生児よりも株を共有していないが、兄弟姉妹は無関係な幼児よりも有意に多くの株を共有していることを見出した。帝王切開で生まれた乳児は、経膣分娩で生まれた乳児よりも有意に高い塩基多様性を持つKlebsiella を保有しており、これは母親のマイクロバイオームではなく、病院からの取得を反映している可能性がある。個々の乳児で多様性を示すゲノム遺伝子座は、他の乳児の間で見られるバリアントを含んでおり、おそらく病院に関連した多様なソースからの接種を反映していると考えられる。
Githubより
inStrainは、メタゲノムから共起する（co-occurring）ゲノム集団を解析するためのpythonプログラムで、高精度なゲノム比較、カバレッジ、マイクロダイバーシティ、連鎖の解析、遺伝子局在や同義非同義識別による高感度SNP検出が可能です。
Documentation
https://instrain.readthedocs.io/en/latest/
The pub describing our new program, inStrain, is out today in Nature Biotechnology https://t.co/ScyoMH9maz . We benchmark inStrain with several real datasets (see docs for overview https://t.co/DcLroII3Pj ) and use it to identify patterns of strain sharing between infants (1/5)
— Matt Olm (@MattagenOlmics) January 18, 2021
We introduce the metric "popANI" (population ANI) to handle genomic comparisons between populations that contain microdiversity. Using a stringent 99.999% popANI "strain" definition, we find that infant strain acquisition can be linked to (1) family-specific sources... (2/5)
— Matt Olm (@MattagenOlmics) January 18, 2021
インストール
#conda(link)
mamba create -n instrain -y
conda activate instrain
mamba install -c conda-forge -c bioconda -c defaults instrain -y
#pip
pip install instrain
> inStrain -h
...::: inStrain v1.5.7 :::...
Matt Olm and Alex Crits-Christoph. MIT License. Banfield Lab, UC Berkeley. 2021
Choose one of the operations below for more detailed help. See https://instrain.readthedocs.io for documentation.
Example: inStrain profile -h
Main operations:
profile -> Create an inStrain profile (microdiversity analysis) from a mapping file
compare -> Compare multiple inStrain profiles (popANI, coverage_overlap, etc.)
Auxiliary operations:
check_deps -> Print a list of dependencies, versions, and whether they're working
quick_profile -> Quickly calculate coverage and breadth of a mapping using coverM
filter_reads -> Commands related to filtering reads from .bam files
plot -> Make figures from the results of "profile" or "compare"
other -> Other miscellaneous operations
profile_genes -> [DEPRECATED; functionality within PROFILE] Calculate gene-level metrics on an inStrain profile
genome_wide -> [DEPRECATED; functionality within PROFILE / COMPARE] Calculate genome-level metrics
> inStrain profile -h
usage: inStrain profile [-o OUTPUT] [--use_full_fasta_header] [--force_compress] [-p PROCESSES] [-d] [-h] [--version] [-l MIN_READ_ANI] [--min_mapq MIN_MAPQ] [--max_insert_relative MAX_INSERT_RELATIVE] [--min_insert MIN_INSERT]
[--pairing_filter {all_reads,paired_only,non_discordant}] [--priority_reads PRIORITY_READS] [--detailed_mapping_info] [-c MIN_COV] [-f MIN_FREQ] [-fdr FDR] [-g GENE_FILE] [-s [STB [STB ...]]] [--mm_level] [--skip_mm_profiling]
[--database_mode] [--min_scaffold_reads MIN_SCAFFOLD_READS] [--min_genome_coverage MIN_GENOME_COVERAGE] [--min_snp MIN_SNP] [--store_everything] [--scaffolds_to_profile SCAFFOLDS_TO_PROFILE] [--rarefied_coverage RAREFIED_COVERAGE]
[--window_length WINDOW_LENGTH] [--skip_genome_wide] [--skip_plot_generation]
bam fasta
REQUIRED:
bam Sorted .bam file
fasta Fasta file the bam is mapped to
I/O PARAMETERS:
-o OUTPUT, --output OUTPUT
Output prefix (default: inStrain)
--use_full_fasta_header
Instead of using the fasta ID (space in header before space), use the full header. Needed for some mapping tools (including bbMap) (default: False)
--force_compress Force compression of all output files (default: False)
SYSTEM PARAMETERS:
-p PROCESSES, --processes PROCESSES
Number of processes to use (default: 6)
-d, --debug Make extra debugging output (default: False)
-h, --help show this help message and exit
--version show program's version number and exit
READ FILTERING OPTIONS:
-l MIN_READ_ANI, --min_read_ani MIN_READ_ANI
Minimum percent identity of read pairs to consensus to use the reads. Must be >, not >= (default: 0.95)
--min_mapq MIN_MAPQ Minimum mapq score of EITHER read in a pair to use that pair. Must be >, not >= (default: -1)
--max_insert_relative MAX_INSERT_RELATIVE
Multiplier to determine maximum insert size between two reads - default is to use 3x median insert size. Must be >, not >= (default: 3)
--min_insert MIN_INSERT
Minimum insert size between two reads - default is 50 bp. If two reads are 50bp each and overlap completely, their insert will be 50. Must be >, not >= (default: 50)
--pairing_filter {all_reads,paired_only,non_discordant}
How should paired reads be handled?
paired_only = Only paired reads are retained
non_discordant = Keep all paired reads and singleton reads that map to a single scaffold
all_reads = Keep all reads regardless of pairing status (NOT RECOMMENDED; See documentation for deatils)
(default: paired_only)
--priority_reads PRIORITY_READS
The location of a list of reads that should be retained regardless of pairing status (for example long reads or merged reads). This can be a .fastq file or text file with list of read names (will assume file is compressed if ends in .gz
(default: None)
READ OUTPUT OPTIONS:
--detailed_mapping_info
Make a detailed read report indicating deatils about each individual mapped read (default: False)
VARIANT CALLING OPTIONS:
Minimum coverage to call an variant (default: 5)
-f MIN_FREQ, --min_freq MIN_FREQ
Minimum SNP frequency to confirm a SNV (both this AND the FDR snp count cutoff must be true to call a SNP). (default: 0.05)
-fdr FDR, --fdr FDR SNP false discovery rate- based on simulation data with a 0.1 percent error rate (Q30) (default: 1e-06)
GENE PROFILING OPTIONS:
-g GENE_FILE, --gene_file GENE_FILE
Path to prodigal .fna genes file. If file ends in .gb or .gbk, will treat as a genbank file (EXPERIMENTAL; the name of the gene must be in the gene qualifier) (default: None)
GENOME WIDE OPTIONS:
-s [STB [STB ...]], --stb [STB [STB ...]]
Scaffold to bin. This can be a file with each line listing a scaffold and a bin name, tab-seperated. This can also be a space-seperated list of .fasta files, with one genome per .fasta file. If nothing is provided, all scaffolds will be
treated as belonging to the same genome (default: [])
READ ANI OPTIONS:
--mm_level Create output files on the mm level (see documentation for info) (default: False)
--skip_mm_profiling Dont perform analysis on an mm level; saves RAM and time; impacts plots and raw_data (default: False)
PROFILE OPTIONS:
--database_mode Set a number of parameters to values appropriate for mapping to a large fasta file. Will set: --min_read_ani 0.92 --skip_mm_profiling --min_genome_coverage 1 (default: False)
--min_scaffold_reads MIN_SCAFFOLD_READS
Minimum number of reads mapping to a scaffold to proceed with profiling it (default: 1)
--min_genome_coverage MIN_GENOME_COVERAGE
Minimum number of reads mapping to a genome to proceed with profiling it. MUST profile .stb if this is set (default: 0)
--min_snp MIN_SNP Absolute minimum number of reads connecting two SNPs to calculate LD between them. (default: 20)
--store_everything Store intermediate dictionaries in the pickle file; will result in significantly more RAM and disk usage (default: False)
--scaffolds_to_profile SCAFFOLDS_TO_PROFILE
Path to a file containing a list of scaffolds to profile- if provided will ONLY profile those scaffolds (default: None)
--rarefied_coverage RAREFIED_COVERAGE
When calculating nucleotide diversity, also calculate a rarefied version with this much coverage (default: 50)
--window_length WINDOW_LENGTH
Break scaffolds into windows of this length when profiling (default: 10000)
OTHER OPTIONS:
--skip_genome_wide Do not generate tables that consider groups of scaffolds belonging to genomes (default: False)
--skip_plot_generation
Do not make plots (default: False)
実行方法
主な入力は、リファレンスゲノム配列のfastaファイル（MAG由来のbin配列など）と、これらの配列にマッピングされたリードを含むbamファイルの2つ。オプションとして、遺伝子.fnaファイルを提供することで、inStrainは遺伝子レベルのメトリックスを計算することができる。
inStrain profile mapping.bam genome_file.fasta -o inStrain_profile1
マニュアルより
de novo アセンブルゲノムの使用（推奨）
- 環境から収集した各サンプルについて、リードをコンティグにアセンブルする。推奨ソフトウェア：IDBA_UD, MEGAHIT, metaSPADES.
- 各アセンブリからゲノムをビニング（Differential coverage binning）する。推奨ソフトウェア：Bowtie2（マッピング用）、MetaBAT、CONCOCT、DasTOOL（ビニング用）。
- プロファイリングしたい全ゲノムセット（全環境の全ゲノム）を 97-99% アイデンティティでdereplicate し、低品質ゲノムをフィルタリングする。推奨ソフトウェア：dRep, checkM.
- ゲノムセットから scaffold-to-bin ファイル（タブで区切られた2列の.textファイルで、1列目はscaffolds名、2列目はそのscaffoldsが属するbin/genome名）を作成する。推奨ソフトウエア：parse_stb.py
- このdereplicatedセットから代表的なゲノムのBowtie2インデックスを作成し、各サンプルからこのセットへのリードをマッピングする。推奨ソフトウェア：Bowtie2
- 得られたマッピング.bamファイルをinStrainでプロファイリングし、binninbされた各ゲノムについて、ゲノムレベルのmicrodiveristyメトリックスを計算する。
このワークフローで重要なことは、同時に多くのゲノムにマッピングすることである。一度に1つのゲノムだけにマッピングするのは、他のゲノムからのミスマッピングがを引き起こすため、全く推奨されない。bowtie2インデックスに多くの（dereplicatedされていない）ゲノムを含めることで、ミスマップリードをより正確にフィルタリングし、偽陽性SNPを減らすことができる。
１、drepのparse_stb.pyを実行してステップ４のscaffold-to-bin ファイルを得る。binned.fastaのディレクトリを指定する。ここではdrepの出力を指定した。
parse_stb.py --reverse -f dereplicated_genomes/*fa -o representitve_genomes.stb
出力
２、リファレンスにマッピングする。
bowtie2-build -f assembly.fasta bowtie2_index
bowtie2 -p 36 -x bowtie2_index sample_R1.fq.gz sample_R2.fq.gz> sample.sam
３、（任意）inStrainに遺伝子レベルのプロファイリングをさせたい場合は、プロファイリングする遺伝子のリストを与える。prodigalを実行する。
prodigal -i binned_genome.fasta -d binned_genome_genes.fna
５、inStrainの実行。２のsamファイル、binned.fasta、出力ディレクトリ名、任意のアノテーションファイル、１からのscaffold-to-bin ファイルを指定する。
inStrain profile 32.sam MAG_bin5_selected.fa -o outdir -p 16 -g genes.fna -s representitve_genomes.stb
- -o Output prefix (default: inStrain)
- -p Number of processes to use (default: 6)
-
-g Path to prodigal .fna genes file. If file ends in .gb or .gbk, will treat as a genbank file (EXPERIMENTAL; the name of the gene must be in the gene qualifier) (default: None)
- -s Scaffold to bin. This can be a file with each line listing a scaffold and a bin name, tab-seperated. This can also be a space-seperated list of .fasta files, with one genome per .fasta file. If nothing is provided, all scaffolds will be treated as belonging to the same genome (default: [])
出力例（マニュアル）
outdir/
output/
scaffold_info.tsv
サンプルに含まれるスキャフォールドの基本的な情報；リードの同一性を確認するため
長さ、リードのデプスの平均と中央値と偏差、リードでカバーされている割合、全塩基の塩基多様性値など
mapping_info.tsv
各スキャfホールドにマッピングされたリードの数、およびその品質に関する基本的なメトリクス。
マッピングされたペアエンドリードの生の数、min_read_aniカットオフをペアエンドリード数、マッピングされたペアエンドリードのうち、mapQスコアの最小値カットオフをパスしたリード数、mapQスコアの平均値、ミスマッチの平均数、ペアのうち片方のリードだけがマップされた数、ペアエンドの平均インサートサイズ、
SNVs.tsv
gene_info.tsv
プロファイリング対象の遺伝子に関する基本的な情報
Genome_info.tsv
gene_info.tsvのメトリクスの多くを、スカフォールド単位ではなく、ゲノム単位で記述。
figures
CoverageAndBreadth_vs_readMismatch.pdf
Coverageとbread mismatchesの比較プロット
genomeWide_microdiveristy_metrics
塩基多様性とSNV密度のスパイクは、カバレッジの増加には対応していないため、リードのミスマッピングによるシグナルではないことが分かる。
GeneHistogram_plot.pdf
非同義のSNPのペアと全SNPのペアの連鎖プロット（上のが画像ではSNPs図が出力されていない）
eadANI_distribution.pdf
参照ゲノムにマッピングされたリードペアのANIレベルの分布。
ReadFiltering_plot.pdf
フィルタリングで除外されたリードの数を棒グラフで示したもの（パーセンテージはすべてペアのリードの数）。
MajorAllele_frequency_plot.pdf
bi-allelic SNVsの頻度スペクトラム。
LinkageDecay_plot.pdf
SNVの連鎖とSNV間の距離の指標。連鎖の減衰（一方のプロットで示され、他方では示されない）は、組換えの一般的なシグナルである。
ScaffoldInspection_plot
ゲノムワイドのマイクロダイバーシティメトリクスを細長くしたもの
inStrain compareコマンドを使うと、全サンプルを popANI に基づいて、また shared_bases に基づいて比較したデンドログラムも出力される。
popANI（population-level ANI）について（論文より）
- inStrainはゲノム比較の際にメジャーアレルとマイナーアレルの両方を考慮する。この新しい微小多様性を考慮したANI指標は「PopANI」（集団レベルのANI）と呼ばれ、コンセンサスベースのANI（「ConANI」）と共に報告される。両指標の計算では、まず、両サンプルで少なくとも5倍のカバレッジがあるゲノムの位置がすべて特定される。PopANIとConANIの計算では、これらの位置のみが考慮される。これらの位置のうち、サンプル間で対立遺伝子組成が異なる位置の数を列挙する。ConANIでは、コンセンサス塩基が2つのサンプル間で異なる場合、置換がコールされる。PopANIでは、両サンプルに対立遺伝子（メジャーまたはマイナー）がない場合のみ、その部位で置換がコールされる。
- 同一株と分類するために、論文中では 99.99996% popANIや99.999% popANIが選択されている。
実践的なチュートリアルが用意されています。マニュアルを確認して下さい。
引用
inStrain profiles population microdiversity from metagenomic data and sensitively detects shared microbial strains
Matthew R. Olm, Alexander Crits-Christoph, Keith Bouma-Gregson, Brian A. Firek, Michael J. Morowitz & Jillian F. Banfield
Nature Biotechnology volume 39, pages727–736 (2021)
関連
