VIRGO 是一个非冗余的基因目录，用于描述人类阴道微生物群落。它是通过结合宏基因组和泌尿生殖道分离基因组序列构建的。VIRGO 可用于描述阴道宏基因组和宏转录组数据集的分类和/或功能组成。可用的功能注释包括：KEGG、COG、eggNOG、CDD、GO、Pfam、TIGRfam 和 EC 编号。   VOG 是一个全面的阴道微生物群正交基因集合，通过一种新颖的 Jaccard 聚类方法定义。这个氨基酸序列的补充数据库可用于改进功能注释、比较基因组学和阴道正交蛋白家族的进化研究。  

如有任何问题或建议，请联系作者。谢谢并祝您使用愉快！

---

所需依赖项

BLAST 2.2.3+ Bowtie 1.1.2+ GNU Awk 3.1+ seqtk 1.0+

---

数据库结构

• _test_run 包含测试数据集和结果
• 0_db 包含可搜索的数据库文件和序列
  • VIRGO Bowtie 1 数据库
  • VIRGO Bowtie 2 数据库
  • VIRGO BLASTN 数据库
  • VIRGO BLASTP 数据库
• 1_VIRGO 包含基因的注释文件
  • 0.geneLength 包含基因长度（bp）
  • 1.taxon.table 包含分类信息
  • 2.A.proteinFamily.txt 包含来自 17 个来源的注释，每行一个来源
  • 2.B.proteinFamily.txt 包含来自 17 个来源的注释，每行一个基因
  • 3.eggnog.txt 包含 eggNOG 注释
  • 4.GC.txt 包含高或低基因计数信息，指示基因是否偏好高或低基因计数群落（>95%）
  • 5.EC.txt 包含酶委员会编号（EC 编号）
  • 6.product.txt 包含基因产物注释
  • 7.rxn.txt 包含 KEGG 反应注释
  • 8.A.kegg.ortholog.txt 包含 KEGG 正交注释
  • 8.B.kegg.module.txt 包含 KEGG 模块注释
  • 8.C.kegg.pathway.txt 包含 KEGG 通路注释
• 2_fungal_phage 包含 10,908 个真菌和 15,965 个噬菌体基因的注释
• 3_run_VIRGO 包含用于分析宏基因组和宏转录组数据集的脚本
• 4_Jaccard_index_clustering 包含用于执行 Jaccard 聚类以识别 VOG 的脚本
• 5_VOG 包含 VOG 序列和链接到 VIRGO 基因的文件，内容包括：
  • 0.VOG.VIRGO.tbl.txt VIRGO 和 VOG ID 之间的链接
  • 1.VOG.size.txt VOG 蛋白家族的大小
  • 2.VOG.GC.txt 基因计数注释
  • 3.VOG.funcat.txt eggNOG 功能类别
  • 4.VOG.taxonomy.txt 共识分类
  • 5.VOG.alignmentScore.txt 蛋白家族的对齐分数

---

本模块介绍如何测试 VIRGO 以确保其正常工作。

目录 _test_run 包含测试数据集：

cd _test_run

运行 step1 和 step2 脚本，确保指定 VIRGO 的下载路径（应包含数据库结构，如 /path/to/VIRGO/0_db 和 /path/to/VIRGO/1_VIRGO）

./runTesting.step1.sh -1 sub1.fq -2 sub2.fq -p test -d /path/to/VIRGO/ ./runTesting.step2.sh -p temp_testsample -d /path/to/VIRGO/

您应该看到以下输出文件！

ls temp_testsample/

• summary.Abundance.txt
• summary.Percentage.txt
• summary.Count.txt
• summary.geneRichness.txt
• test.1.reads2ref
• test.annotation.txt
• test.1.fq
• test.2.fq
• test.out

---

本模块演示如何将宏基因组或宏转录组样本的读取映射到 VIRGO 并展示结果。

请注意，VIRGO 尚不支持双端映射，请将双端读取单独或与单端读取合并为一个 fastq 文件。

runMapping.Step1.sh
参数：

-r 包含单端读取的 fastq 文件 -p 样本前缀 -d VIRGO 的路径

示例：

./runMapping.step1.sh -r sample1.fastq -p samplePrefix -d /full/path/to/VIRGO/

输出应为 sampleName.out 文本文件。第一列为 VIRGO 基因 ID，第二列为映射到 VIRGO 数据库的读取数，第三列为基因长度。文件按映射读取数列排序。

head sampleName.out

V1593031  2425  3663
V1607456  566 390
V1420626  348 1017
V1684181  327 1442
V1316648  324 480
V1541216  274 1386
V1064785  269 303
V1872750  246 1179
V1885045  227 1254
V1573770  226 291

对所有样本集重复此操作。以下是伪代码示例：

for sample in *.fq; do ./runMapping.step1.sh -r sample1.fastq -p $sample -d /virgo/path/; done

汇总多个样本的统计信息

./runMapping.step2.sh -p /path/to/output/of/step1/ -d /path/to/VIRGO/

输出包括：

• summary.Abundance.txt: 每个物种的读取数
• summary.Count.txt: 每个物种的基因数
• summary.Percentage.txt: 每个物种的标准化丰度百分比（总和为 100）
• summary.geneRichness.txt: 每个样本的基因数
• summary.NR.abundance.txt: 每个非冗余基因的读取数
• gene.lst.txt: 非冗余基因列表及其长度
• EggNOG.annotation.txt: 每个样本的 EggNOG 注释文件
• EC.annotation.txt: 非冗余基因列表及其 EC 编号
• GC.txt: 非冗余基因列表及其基因计数类别（HGC：高基因计数，LGC：低基因计数）
• geneProduct.txt: 非冗余基因列表及其基因产物注释
• Kegg.module.annotation.txt: 非冗余基因列表及其 KEGG 模块注释，包括模块 ID 和注释
• Kegg.ortholog.annotation.txt: 非冗余基因列表及其 KEGG 正交（KO）注释，包括正交 ID 和注释
• Kegg.pathway.annotation.txt: 非冗余基因列表及其 KEGG 通路注释，包括通路 ID 和注释
• proteinFamily.annotation.txt: 非冗余基因列表及其蛋白家族注释，来自 CDD、GO、Gene3D、Hamap、Interpro、MobiDBLite、PIRSF、PRINTS、Pfam、ProDom、ProSitePatterns、ProSiteProfiles、SFLD、SMART、SUPERFAMILY、TIGRFAM 数据库
• rxn.annotation.txt: 非冗余基因列表及其 KEGG 反应

---

需要对 bash 脚本和简化映射结果有基本了解

可调整的参数：

• "reads = " 默认值为 2：在 step1 脚本中，调用基因命中的读取数阈值。
• "bowtie -l" 默认值为 25：在 step1 脚本中，-l 是读取映射的种子长度。可以降低以获得更多命中，或提高以提高映射质量。
• "--VIRGO-path" VIRGO 目录的完整路径：硬编码路径，以便每次运行时无需指定。
• 输出目录和输出名称：脚本硬编码了输出名称和功能注释。可以自由命名更好的名称。
• 注释文件：可以为 VIRGO 数据库的原始 fasta 文件生成自己的注释。

---

使用单个 VIRGO geneID 可以快速提取您感兴趣的基因的注释！通过使用 "grep" 命令与 VIRGO 的注释文件，您可以获取大量信息。例如，对于基因 V1000057，您可以提取其长度、分类、GO、TIGRFAM、Interpro、SUPERFAMILY、CDD、Pfam、eggNOG、EC 编号、基因产物、KEGG 正交、KEGG 模块、KEGG 通路等。

grep -w "V1000057" 1_VIRGO/*

• 0.geneLength.txt :V1000057 990
• 1.taxon.tbl.txt :Cluster_219398 V1000057 Lactobacillus_iners 990
• 2.A.proteinFamily.txt :V1000057 GO ; 磷酸核糖转移酶结构域；磷酸核糖转移酶类似；核糖-磷酸二磷酸激酶；核糖-磷酸焦磷酸激酶，N 端结构域
• 2.A.proteinFamily.txt :V1000057 TIGRFAM ; 核糖-磷酸二磷酸激酶
• 2.A.proteinFamily.txt :V1000057 Interpro ; 磷酸核糖转移酶结构域；磷酸核糖转移酶类似；核糖-磷酸二磷酸激酶；核糖-磷酸焦磷酸激酶，N 端结构域
• 2.A.proteinFamily.txt :V1000057 SUPERFAMILY ; PRTase 类似；PRTase 类似
• 2.A.proteinFamily.txt :V1000057 CDD ; PRTases_typeI
• 2.A.proteinFamily.txt :V1000057 Pfam ; 核糖磷酸焦磷酸激酶的 N 端结构域；磷酸核糖合成酶相关结构域
• 2.B.proteinFamily.txt :V1000057 ; PRTases_typeI ; 磷酸核糖转移酶结构域；磷酸核糖转移酶类似；核糖-磷酸二磷酸激酶；核糖-磷酸焦磷酸激酶，N 端结构域 ; 磷酸核糖转移酶结构域；磷酸核糖转移酶类似；核糖-磷酸二磷酸激酶；核糖-磷酸焦磷酸激酶，N 端结构域 ; 核糖磷酸焦磷酸激酶的 N 端结构域；磷酸核糖合成酶相关结构域 ; PRTase 类似；PRTase 类似 ; 核糖-磷酸二磷酸激酶
• 3.eggnog.NOG.txt :Cluster_219398 V1000057 PRS map00030,map00230,map01100,map01110,map01120,map01230 F 磷酸核糖焦磷酸合成酶 COG0462
• 5.EC.txt :Cluster_219398 V1000057 ec:2.7.6.1 核糖-磷酸二磷酸激酶；核糖-磷酸焦磷酸激酶；PRPP 合成酶；磷酸核糖焦磷酸合成酶；PPRibP 合成酶；PP-核糖 P 合成酶；5-磷酸核糖-1-焦磷酸合成酶；5-磷酸核糖焦磷酸化酶；5-磷酸核糖-α-1-焦磷酸合成酶；磷酸核糖-二磷酸合成酶；磷酸核糖焦磷酸合成酶；焦磷酸核糖磷酸合成酶；核糖磷酸焦磷酸激酶；核糖-5-磷酸焦磷酸激酶
• 6.product.txt :Cluster_219398 V1000057 核糖-磷酸焦磷酸激酶
• 8.A.kegg.ortholog.txt :V1000057 K00948 ljf:FI9785_163 dvvi:GSVIVP00018977001 GSVIVT00018977001; 组装 CDS；K00948 核糖-磷酸焦磷酸激酶 [EC:2.7.6.1]
• 8.B.kegg.module.txt :V1000057 M00005 PRPP 生物合成，核糖 5P => PRPP
• 8.C.kegg.pathway.txt :V1000057 map00230 嘌呤代谢

---

可以使用 VIRGO ID 列表提取感兴趣的基因信息。

步骤 1. 生成 VIRGO ID 的基因列表。例如：

cat lst

V1891288
V1891259
V1891167
V1891137
V1891104

步骤 2. 提取这些基因的注释，如分类、功能注释、eggNOG、KEGG 等。

awk -F"\t" 'NR==FNR{a[$2]=$3;next} ($1 in a) {print $0"\t"a[$1]}' 1_VIRGO/1.taxon.tbl.txt lst

V1891288  Lactobacillus_crispatus
V1891259  Lactobacillus_crispatus
V1891167  Lactobacillus_crispatus
V1891137  Lactobacillus_crispatus
V1891104  Lactobacillus_crispatus

---

VIRGO 也可以用于 BLASTN 和 BLASTP。只需运行常规的 blastn 或 blastp 命令。

兼容 blast 2.2.3+

---

步骤 1 在蛋白质集上运行全对全 BLASTP

步骤 2 使用修改后的 clusterBsmlPairwiseAlignments 脚本运行 Jaccard 聚类

./run-jaccard.sh

脚本/程序：

clusterBsmlPairwiseAlignments-mod.pl

输入：

• Wu-blastp.bsml.list（BLASTP BSML 文件列表）
• 从多肽对生成簇的可执行文件位置
• 各种参数（linkscore、percent_identity、percent_coverage、p_value）

输出：

• Jaccard-clusters.out(.gz)

clusterBsmlPairwiseAlignments-mod.pl:

• 读取并过滤所有 BLASTP 数据，生成多肽对列表
• 将多肽对写入文件（20808.pairs.gz）
• 调用 /usr/local/projects/ergatis/package-latest/bin/cluster 将多肽对转换为簇

注意，所有聚类/过滤参数在多肽对生成步骤中应用，而不是在多肽对聚类步骤中。

示例：

clusterBsmlPairwiseAlignments-mod.pl
--bsmlSearchList=wu-blastp.bsml.list
--log=test-jaccard.log
--linkscore=0.6
--percent_identity=80
--percent_coverage=60
--p_value=1e-5
--outfile=jaccard-clusters.out
--cluster_path=/path/to/cluster/files/

步骤 3. 在 Jaccard 簇上运行双向最佳命中分析。

./run-best-hit.sh

3a. 将 BLASTP BSML 转换为最佳命中的单个 btab 文件。

脚本/程序： CogBsmlLoader-mod.pl

输入：

• Wu-blastp.bsml.list（BLASTP BSML 文件列表）
• Jaccard-clusters.out（来自 Jaccard 聚类步骤）
• 参数（最小 p 值和最小覆盖率）

输出：

• 汇总 btab 文件（/path/to/btab） （最后一列包含 Jaccard 簇信息）

3b. 合并所有具有相互最佳命中的簇。

脚本/程序： best_hit.pl

输入：

• 来自步骤 3a 的 All-best-hits.btab
• 最小 Jaccard 系数（如果设置为 > 0，用于额外的簇修剪）

输出：

• 最终簇文件（final-clusters.cog.gz）

4. 生成 FASTA 文件。

脚本/程序： makeClusterFASTAFiles.pl

输入：

• 原始多肽文件（如 total.faa）
• 来自步骤 2 的 Kaccard-clusters.out
• 来自步骤 3b 的 Final-clusters.cog
• 输出目录的路径

输出：

• 每个大小 > 1 的簇的 FASTA 文件
• 一个 singletons.fa 文件，包含所有大小为 1 的簇
• 簇大小直方图（在 stdout 上 - 参见 make-fasta-files-2.log.gz）