软件教程 | StrainR2-宏基因组定量菌株丰度

2025-04-08

📅 官方地址:https://github.com/BisanzLab/StrainR2

🔬 教程版本:V2.2.1


📋 主要功能

量化微生物组中菌株丰度的传统方法(如16S rRNA测序)缺乏区分菌株的分辨率,只能粗略归类到物种层面。虽然宏基因组鸟枪法测序提供了替代方案,但FPKM等未标准化的丰度指标会因相似基因组获得较少独特比对而产生偏差。

为消除这种偏差开发的其他宏基因组方法(如NinjaMapStrainR1)存在精度不足的问题,且难以扩展至更大规模的菌群分析。尤其在菌株覆盖率较低时,丰度测算误差可达数量级差异,计算资源的扩展性也较差。

StrainR2作为解决方案,通过宏基因组鸟枪法测序能在传统方法所需时间的一小部分内,将菌株丰度定量误差控制在真实值的几个百分点以内。


🔧 依赖环境:

BBMap
fastp
GNU make (if compiling from source)
samtools
bedtools
R
R optparse
R tidyverse
zlib

🧷 安装方式:

StrainR2专为类Unix操作系统设计。依赖项列于本README文件底部,但强烈建议使用conda以确保可复现性。

# Option A: Bioconda (recommended) √
conda create -n strainr2 -c bioconda -c conda-forge strainr2
conda activate strainr2

# Option B: Docker
docker pull quay.io/biocontainers/strainr2:<tag>

# Option C: Installation from source
git clone https://github.com/BisanzLab/StrainR2.git
cd StrainR2/src
make release

🚩 主要命令介绍:

StrainR2分为两个步骤:PreProcessR和StrainR。

PreProcessR

PreProcessR为后续StrainR的运行创建数据库。它会将基因组contigs分割为subcontigs以确保相似的基因组构建质量,长度低于预设值(默认为10kbp)的contigs会被标记为”exluded”。独特k-mer的计数信息会存储在KmerContent.report中。此外,预处理过程会生成BBMap索引供后续使用。所有这些文件都将存放在指定的输出目录中。

对于特定的基因组群落,PreProcessR只需运行一次,其输出结果可被StrainR无限次重复使用。处理大规模数据时,可能需要增加交换空间大小以满足内存需求。这种情况仅在PreProcessR因超出内存限制崩溃时才需要调整。

可通过以下命令行方式调用PreProcessR:

PreProcessR -i <PATH_TO_GENOME_FILES> [OPTIONS]

参数说明:

-i 或 –indir [必需]:

包含待量化基因组文件的目录路径。仅接受.fna和.fasta格式的文件。

-o 或 –outdir:

输出文件的目标存储目录路径。运行PreProcessR前无需预先创建该目录。若输出目录在运行前已存在,请确保其为空目录。默认路径为./StrainR2DB。

-e 或 –excludesize:

子片段的最小排除尺寸。仅长度小于此值的contig会被排除。默认值=10,000。

-s 或 –subcontigsize:

此选项将覆盖默认的N50构建质量统计值(由StrainR2自动计算所有输入基因组得出)。StrainR2将以给定尺寸作为最大值创建最优子片段。建议保持默认最大尺寸。较小的子片段会降低StrainR2对低丰度菌株的敏感度,而过大的子片段会导致过度敏感——可能因索引跳读、测序错误和其他假读源而误报不存在菌株的丰度。

-r 或 –readsize:

双端测序中单端读长的尺寸。默认值150。用于计算k-mer大小。

StrainR

StrainR接收双端测序读段,并利用PreProcessR步骤准备的输出数据对菌株丰度进行标准化处理。该工具会生成.sam、.rpkm和.abundance文件。其中.sam和.rpkm文件通过BBMap工具生成,而.abundance文件则通过先前生成的KmerContent.report文件中的数据对.rpkm输出结果进行标准化。此外,StrainR还会生成一张丰度分布图。

可通过以下命令行方式调用StrainR命令:

StrainR -1 <PATH_TO_FORWARD_READS> -2 <PATH_TO_REVERSE_READS> -r <PATH_TO_OUTPUT_OF_PREPROCESSR>[OPTIONS]

参数说明:

-1 或 –forward [必需]:

正向读段的路径(读段无需为.gz格式)

-2 或 –reverse [必需]:

反向读段的路径(读段无需为.gz格式)

-r 或 –reference [必需]:

PreProcessR生成的输出目录路径(即PreProcessR的输出目录)

-c 或 –weightedpercentile:

用于最终丰度计算的菌株子contig FUKMs的加权百分位数(按独特k-mer数量加权)。提高此值将增加灵敏度但可能增加假阳性。降低此值则可能增加假阴性。默认=60

-s 或 –subcontigfilter:

基于独特k-mer数量过少而需要过滤掉的子contig百分比。默认=0

-b1 或 –background1:

需从分析中排除的背景正向读段路径。通常不用于合成群落的标准场景,该选项用于过滤未定义群落的背景读段(读段无需为.gz格式)

-b2 或 –background2:

需从分析中排除的背景反向读段路径。通常不用于合成群落的标准场景,该选项用于过滤未定义群落的背景读段(读段无需为.gz格式)

-o 或 –outdir:

存放标准化丰度结果的输出目录路径。默认=当前目录

-p 或 –prefix:

群落名称(用于输出文件名)。默认=”sample”

-t 或 –threads number:

运行fastp、BBMap和samtools时使用的线程数。最大值为16,默认=8

-m 或 –mem number:

运行BBMap时使用的内存(GB)。默认=8

📰 软件测试:

这里将使用单元测试数据进行最小化演示,测试数据包括位于tests/genomes/multiple_complete目录中的8个参考基因组,以及tests/inputs/目录下的模拟输入文件。

注:如使用Bioconda方式安装,则需要另外下载测试数据。

git clone https://github.com/BisanzLab/StrainR2.git

📋 测试数据:

# 基因组文件:
StrainR2/tests/genomes/multiple_complete/
# 宏基因组测序文件:
StrainR2/tests/inputs/mock_reads_testing_R1.fastq.gz
StrainR2/tests/inputs/mock_reads_testing_R2.fastq.gz

📑 测试代码:

PreProcessR --indir StrainR2/tests/genomes/multiple_complete/ --outdir database --readsize 150

StrainR \
 --forward StrainR2/tests/inputs/mock_reads_testing_R1.fastq.gz \
 --reverse StrainR2/tests/inputs/mock_reads_testing_R2.fastq.gz \
 --reference database/ \
 --prefix testrun \
 --outdir strain2_output

📣 结果解释:

在生成的文件输出目录strainr2_output中,可以找到testrun_abundance_summary.tsv结果文件,其中有菌株的丰度信息:

StrainID	weighted_percentile_FUKM	median_FUKM	sd_FUKM	subcontigs_detected	subcontigs_total	percent_detected	percent_abundance
JEB00015	15.002476052670442	14.949939922104155	0.46622826476923757	38	38	100	13.524840656749035
JEB00022	14.885728923252971	14.865410790723015	0.6712969280564985	52	52	100	13.41959227528449
JEB00024	15.005580606237396	14.972944694458326	0.9742879632461824	42	42	100	13.527639434241228
JEB00030	15.034270854543474	14.829354067500985	0.5823678254771962	28	28	100	13.553503900571984
JEB00031	15.124317636874952	15.107007945963842	0.6698836185646367	30	30	100	13.634681726046297
JEB00036	15.088684933746867	14.775182972247737	1.7746910771776707	48	48	100	13.602558586486888
JEB00041	14.71064817309631	14.689071798986676	0.6518258370527644	26	26	100	13.261755712865087
JEB00052	6.073637031925493	5.977622518629116	0.27601498880837394	50	50	100	5.4754277077549895

StrainID: 将提供汇总统计数据的fasta文件名。

weighted_percentile_FUKM, median_FUKM, sd_FUKM: 菌株中所有子contig FUKM的加权百分位数(根据StrainR的-c选项)、中位数和标准差

subcontigs_detected: 检测到FUKM>0的子contig数量

subcontigs_total: 该菌株的子contig总数

percent_detected: subcontigs_detected / subcontigs_total

percent_abundance: weighted_percentile_FUKM / 群落中所有weighted_percentile_FUKM的总和

还可在testrun.pdf中查看结果的可视化图表,该图表对丰度进行了汇总。在此特定测试数据中,除JEB00052丰度偏低外,多数菌株的丰度水平相近。

image

🏷️ 标签:StrainR2、宏基因组、相对丰度

📅 创建:2025-04-08

🔗 原文:查看文献