33、致病性因子的基因组分析：从E值调整到数据可视化

致病性因子的基因组分析：从E值调整到数据可视化

1. 玩转E值

在基因组分析中，E值是一个关键的设定，它会影响我们的分析结果。E值本质上是用来判断两个蛋白质是否相等的依据。在BLAST搜索中，如果两个蛋白质被认为不相等，就不会出现匹配结果，而这是基于BLAST E值来判断的。在之前的例子中，我们将E值设定为0.00001。

E值的大小与序列的相似性密切相关。序列比对的E值越低，说明比对的序列越相似。在使用blastp时，我们设定的期望阈值（E值）越低，在大肠杆菌O157:H7中找到与给定大肠杆菌K12序列匹配的结果就越少。

为了验证这一假设，我们需要多次使用不同的参数执行blastp。这时，使用一个shell脚本就会让工作变得轻松许多。以下是脚本 autoblast.sh 的代码：

#!/bin/bash
# save as autoblast.sh
# changes E-value
for i in 1 0.1 0.001 0.0001 0.00001
do
./ncbi-blast-2.2.25+/bin/blastp -db ecolik12 -query ecoli-h7.faa -out h7vsk12-$i.txt -evalue $i
done

在这个脚本中，第7 - 11行使用了 for ...;do ...;done 结构，它会遍历第7行中给出的所有期望E值。实际的E值存储在变量 i 中，第10行的 $i 会被变量 i 的当前值所替换。

在运行脚本之前，不要忘记使用 chmod 命令使文件可执行：

$ chmod u+x autoblast.sh
$ time ./autoblast.sh

执行 autoblast.sh 脚本大约需要15分钟（具体时间取决于你的计算机系统），最终会得到五个包含不同E值结果的文件。我们可以观察到，E值越高，结果文件中的行数就越多。这是为什么呢？可以思考一下这个问题。

接下来，我们搜索那些没有匹配结果的序列。为了一次性处理所有结果文件，我们使用BASH shell的 for ...; do ...; done 结构：

$ for i in h7vsk12-*; do echo -n $i" : "; awk '/Query=/ || /No hits/{print}' $i | awk '{i++; line[i]=$0; if($0~/No hits/){print line[i-1]}}' | wc -l; done
$ for i in h7vsk12-*; do echo -n $i" : "; awk '/Query=/ || /No hits/{print}' $i | awk '{i++; line[i]=$0; if($0~/No hits/){print line[i-1]}}' | egrep -v "([Uu]nknown|[Pp]utative|[Hh]ypothetical)" | wc -l; done

通过这些操作，我们可以发现，随着E值的增加，被识别为大肠杆菌O157:H7特有的开放阅读框的数量会减少。你能解释这一现象吗？

2. 使用MySQL组织结果

当处理大型基因组并反复从BLAST+结果文件中提取信息时，将这些结果保存到MySQL数据库中是一个不错的选择。不过，这需要将数据转换为表格格式。幸运的是，BLAST+提供了 -outfmt 选项，可以将BLAST结果转换为不同的格式，具体如下表所示：
| 编号 | 输出格式 |
| ---- | ---- |
| 0 | 成对格式 |
| 1 | 显示匹配信息的查询锚定格式 |
| 2 | 不显示匹配信息的查询锚定格式 |
| 3 | 扁平查询锚定格式，显示匹配信息 |
| 4 | 扁平查询锚定格式，不显示匹配信息 |
| 5 | XML Blast输出 |
| 6 | 表格格式 |
| 7 | 带注释行的表格格式 |
| 8 | 文本ASN.1格式 |
| 9 | 二进制ASN.1格式 |
| 10 | 逗号分隔值格式 |
| 11 | BLAST存档格式（ASN.1） |

以下是修改后的 autoblast.sh 脚本，用于生成不同E值的表格输出：

#!/bin/bash
# save as autoblast.sh
# changes E-value; output tabular
for i in 1 0.1 0.001 0.0001 0.00001
do
./ncbi-blast-2.2.25+/bin/blastp -db ecolik12 -query ecoli-h7.faa -out h7vsk12-$i.tab -evalue $i -outfmt 6
done

这个脚本与之前的脚本类似，只是做了两个修改：一是将文件扩展名改为 .tab ，二是添加了 -outfmt 6 参数，将输出格式设置为纯表格格式。

运行修改后的脚本后，查看其中一个结果文件，会看到长基因标识符和许多包含数字的字段。我们可以通过删除RefSeq ID来简化结果文件：

$ sed -e 's/|ref|NP_[0-9.]\+|//g' h7vsk12-0.00001.tab > h7vsk12-0.00001-gi.tab
$ head -5 h7vsk12-0.00001-gi.tab

BLAST+输出的表格格式各字段含义如下表所示：
| 字段 | 描述 | 缩写 |
| ---- | ---- | ---- |
| 01 | 查询ID | quid |
| 02 | 主题ID | suid |
| 03 | 比对一致性百分比 | iden |
| 04 | 比对长度 | alen |
| 05 | 比对中的错配数 | mism |
| 06 | 比对中的空位开口数 | gapo |
| 07 | 比对中查询序列的第一个核苷酸 | qsta |
| 08 | 比对中查询序列的最后一个核苷酸 | qend |
| 09 | 比对中主题序列的第一个核苷酸 | ssta |
| 10 | 比对中主题序列的最后一个核苷酸 | send |
| 11 | 比对的E值 | eval |
| 12 | 比对的比特分数 | bits |

2.1 将BLAST+结果加载到MySQL中

现在，我们将文件 h7vsk12-0.00001-gi.tab 导入到名为 h7vsk1200001 的MySQL表中。以下是在MySQL命令行中创建表、加载数据和提取特定条目的操作：

mysql> CREATE TABLE h7vsk1200001 (quid VARCHAR(15), suid VARCHAR(15), iden FLOAT, alen INT, mism INT, gapo INT, qsta INT, qend INT, ssta INT, send INT, eval FLOAT, bits INT);
mysql> DESCRIBE h7vsk1200001;
mysql> LOAD DATA LOCAL INFILE "h7vsk12-0.00001-gi.tab" INTO TABLE h7vsk1200001;
mysql> SELECT * FROM h7vsk1200001 WHERE iden = 100 AND bits > 1800;
mysql> SELECT quid, suid, iden, alen FROM h7vsk1200001 WHERE iden = 100 AND alen > 800;

在某些Linux安装中，你可能没有执行 LOAD DATA LOCAL INFILE ... INTO TABLE ... 命令的权限。在这种情况下，可以在终端中执行以下命令：

mysqlimport --local -D -u awkologist -p compbiol h7vsk1200001.tab

不过，需要将 h7vsk12-0.00001-gi.tab 的文件名改为 h7vsk1200001.tab 才能成功导入。

2.2 将基因注释加载到MySQL中

我们会发现表格形式的BLAST结果中并不包含任何基因注释信息。不过，BLAST+套件提供了一个很棒的命令 blastdbcmd ，它可以从BLAST数据库中检索序列和注释信息。在运行 makeblastdb 时，设置 parse_seqids 选项是很重要的。

以下是使用 blastdbcmd 从BLAST数据库中提取信息的示例：

$ ./ncbi-blast-2.2.25+/bin/blastdbcmd -db ecolik12 -entry "all" | head
$ ./ncbi-blast-2.2.25+/bin/blastdbcmd -db ecolik12 -entry "all" -outfmt "gi|%g#%t" | sed -e 's/\[.*/\tK12/' -e 's/#/\t/' | head
$ ./ncbi-blast-2.2.25+/bin/blastdbcmd -db ecolik12 -entry "all" -outfmt "gi|%g#%t" | sed -e 's/\[.*/\tK12/' -e 's/#/\t/' > annotation-k12.tab

接下来，我们使用AWK评估 annotation-k12.tab 文件各字段的长度：

$ awk -F"\t" '{for(i=1;i<=NF;i++){if(length($i)>field[i]){field[i]=length($i)}}} END{for(x=1; x<=i; x++){print field[x]}}' annotation-k12.tab

根据评估结果，我们可以创建一个合适的MySQL表：

mysql> CREATE TABLE ecolianno (gene VARCHAR(12), anno VARCHAR(170), coli VARCHAR(3));
mysql> DESCRIBE ecolianno;
mysql> LOAD DATA LOCAL INFILE "annotation-k12.tab" INTO TABLE ecolianno;
mysql> SELECT * FROM ecolianno LIMIT 5;

在创建和填充表时，要确保在保存 annotation-k12.tab 文件的目录中登录MySQL，否则需要在 LOAD DATA 命令中提供该文件的完整路径。现在，你应该能够对大肠杆菌O157:H7的注释信息进行同样的操作了，你会怎么做呢？

2.3 提取带有注释的高度保守基因

到目前为止，我们已经创建了包含BLAST查询结果和基因注释的两种表。接下来，我们将这两个表（ h7vsk1200001 和 ecolianno ）结合起来，提取BLAST结果中高度保守的基因（至少800个碱基对的序列一致性达到100%）及其注释信息。

mysql> SELECT b.quid, b.suid, b.iden, b.alen, e.anno FROM h7vsk1200001 AS b, ecolianno AS e WHERE b.iden = 100 AND b.alen > 800 AND b.suid = e.gene;

从查询结果中可以看到，只有五个基因符合查询条件，其中两个与转录过程有关。

2.4 提取独特的开放阅读框（ORFs）

你可能会想，是否可以使用MySQL来完成之前用AWK做的事情，即提取大肠杆菌O157:H7或O104:H4致病菌株特有的ORFs。答案是肯定的。

为了找出致病菌株特有的ORFs，我们可以检查注释表 ecolianno.gene 中的哪些ORF ID不在BLAST+结果表 h7vsk1200001.suid 中。这是因为与默认的文本形式的BLAST+结果文件 h7vsk1200001.txt 不同，在表格输出文件 h7vsk1200001.tab 中，在BLAST+数据库（大肠杆菌K12 ORFs）中没有匹配结果的查询序列ID（大肠杆菌O157:H7 ORFs）不会被列出。

以下是尝试提取独特ID的示例：

mysql> SELECT count(e.gene) FROM ecolianno AS e, h7vsk1200001 AS h WHERE e.gene != h.suid;
mysql> SELECT count(e.gene) FROM ecolianno AS e LEFT JOIN h7vsk1200001 AS h ON h.suid = e.gene WHERE h.suid IS NULL LIMIT 5;
mysql> SELECT e.gene, h.suid, e.anno FROM ecolianno AS e LEFT JOIN h7vsk1200001 AS h ON h.suid = e.gene WHERE h.suid IS NULL limit 5;

简单地统计匹配ID的数量（如第一个查询结果为88264608），我们会发现结果存在问题。显然，我们需要优化查询，实际上，必须使用 JOIN 语句。

2.5 使用R可视化和分析结果

现在，让我们使用统计软件包R快速直观地查看BLAST+结果文件。在运行下面的示例之前，需要从CRAN存储库安装R包 DBI 和 RMySQL 。

以下是在R中进行数据可视化和分析的示例：

$ R
> library(RMySQL)
> drv = dbDriver("MySQL")
> con = dbConnect(drv,host="localhost",dbname="compbiol", user="awkologist",pass="awkology")
> data = dbGetQuery(con,statement="SELECT iden, alen, eval FROM h7vsk1200001")
> attach(data)
> plot(iden,alen,xlab="Identity [%]",ylab="Alignment Length [bp]") # Plot1
> abline(lm(alen~iden),col="red")
> summary(lm(alen~iden))
> hist(iden,xlab="Identity [%]",breaks=70) # Plot 2

在这个过程中，我们首先加载 RMySQL 库，然后建立与MySQL数据库的连接，执行查询并将结果存储在变量 data 中。接着，我们绘制了比对一致性百分比与比对长度的散点图，并添加了线性回归直线，同时查看了回归直线的拟合优度信息。最后，我们绘制了序列比对一致性的直方图。

通过直方图可以看到，大多数比对的一致性要么是100%，要么大约是28%。需要注意的是，序列比对的总数与进行比对的蛋白质序列数量并不匹配，你能解释这是为什么吗？

我们还可以绘制序列比对一致性与E值的关系图并进行回归分析：

> plot(iden,eval,xlab="Identity [%]",ylab="E-Value", main="E-Value/Sequence-Identity Relation") # Plot 3
> abline(lm(eval~iden),col="red")
> summary(lm(eval~iden))

从分析结果可以看出，E值和核苷酸一致性百分比之间存在显著的相关性。

最后，我们分析BLAST+的最大比对长度，将其分为三组：

> library(RMySQL)
> drv = dbDriver("MySQL")
> con = dbConnect(drv,host="localhost",dbname="compbiol", user="awkologist",pass="awkology")
> dna<-dbGetQuery(con,statement="SELECT a.alen FROM h7vsk1200001 AS a, ecolianno AS b WHERE anno LIKE '%DNA polymerase%' AND a.suid = b.gene")
> length(dna[[1]])
> dehydrogenase<-dbGetQuery(con,statement="SELECT a.alen FROM h7vsk1200001 AS a, ecolianno AS b WHERE anno LIKE '%dehydrogenase%' AND a.suid = b.gene")
> length(dehydrogenase[[1]])
> cytochrome<-dbGetQuery(con,statement="SELECT a.alen FROM h7vsk1200001 AS a, ecolianno AS b WHERE anno LIKE '%cytochrome%' AND a.suid = b.gene")
> length(cytochrome[[1]])
> boxplot(dehydrogenase[[1]],dna[[1]],cytochrome[[1]],notch=T, names=c("DH","DNA-Pol","CYT"),main="Maximum Alignment Length [bp]")
> t.test(dehydrogenase[[1]],cytochrome[[1]])
> mean(dehydrogenase[[1]])
> mean(cytochrome[[1]])
> q()

在这个过程中，我们分别获取了与“DNA聚合酶”、“脱氢酶”和“细胞色素”相关的蛋白质比对长度数据，创建了箱线图，并进行了t检验。从箱线图的直观观察可以发现，DNA-Pol和CYT的中位数没有差异，而它们似乎都与DH的中位数有显著差异。t检验结果显示，脱氢酶和细胞色素的最大比对长度的均值存在显著差异。你可以检查一下脱氢酶和DNA聚合酶的数据是否也存在相关性。最后，使用 q() 命令退出R时，你可以选择是否保存工作区的图像。

3. 总结

本文详细介绍了在基因组分析中如何调整E值、使用MySQL组织和分析BLAST+结果，以及利用R进行数据可视化和统计分析。通过这些方法，我们可以更深入地了解大肠杆菌不同菌株之间的基因差异和序列相似性。在实际应用中，这些技术可以帮助我们发现致病因子、研究基因功能等。希望这些内容对你有所帮助，你可以根据自己的需求进一步探索和实践。

4. 练习

• 当与大肠杆菌K12菌株相比时，大肠杆菌O157:H7菌株有1099个独特的蛋白质序列。那么，大肠杆菌O157:H7菌株编码了多少个带注释的独特蛋白质序列呢？
• 对于大肠杆菌O104:H4菌株，（a）与大肠杆菌K12菌株相比，（b）与大肠杆菌O157:H7菌株相比，分别能找到多少个独特的带注释蛋白质呢？
• 为什么在将一个基因组与另一个基因组进行比对，或者将一个序列与一组其他序列进行比对时，E值越高，得到的匹配结果就越多呢？
• 为什么随着E值的增加，被识别为查询菌株特有的开放阅读框的数量会减少呢？
• 如何将大肠杆菌O157:H7的注释信息添加到MySQL表 ecolianno 中呢？
• 大肠杆菌O157:H7的基因组大约编码5400个蛋白质，但表格形式的BLAST+输出文件 h7vsk12-0.00001.tab 有超过21000行（即比对结果），这是为什么呢？
• 脱氢酶和DNA聚合酶的最大比对长度之间是否也存在显著的相关性呢？

通过完成这些练习，你可以更好地掌握本文所介绍的知识和技能。

5. 关键技术操作流程总结

为了更清晰地展示整个分析过程，下面用 mermaid 流程图呈现从 E 值调整到数据可视化的主要步骤：

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A([开始]):::startend --> B(调整 E 值):::process
    B --> C(生成 BLAST 结果文件):::process
    C --> D{是否需要 MySQL 存储}:::decision
    D -- 是 --> E(转换为表格格式):::process
    E --> F(创建 MySQL 表):::process
    F --> G(加载 BLAST 结果到 MySQL):::process
    G --> H(加载基因注释到 MySQL):::process
    H --> I(提取高度保守基因及注释):::process
    I --> J(提取独特 ORFs):::process
    D -- 否 --> K(直接分析结果文件):::process
    J --> L{是否需要可视化分析}:::decision
    K --> L
    L -- 是 --> M(使用 R 进行可视化和分析):::process
    L -- 否 --> N([结束]):::startend
    M --> N

6. 各步骤详细操作要点回顾

6.1 E 值调整

• E 值意义 ：E 值用于判断蛋白质是否相等，E 值越低，序列越相似。
• 操作步骤 ：使用 shell 脚本 autoblast.sh 多次执行 blastp，设置不同 E 值，观察结果文件行数变化。

#!/bin/bash
# save as autoblast.sh
# changes E-value
for i in 1 0.1 0.001 0.0001 0.00001
do
./ncbi-blast-2.2.25+/bin/blastp -db ecolik12 -query ecoli-h7.faa -out h7vsk12-$i.txt -evalue $i
done

6.2 MySQL 数据处理

• 表格格式转换 ：使用 BLAST+的 -outfmt 6 选项将结果转换为表格格式。

#!/bin/bash
# save as autoblast.sh
# changes E-value; output tabular
for i in 1 0.1 0.001 0.0001 0.00001
do
./ncbi-blast-2.2.25+/bin/blastp -db ecolik12 -query ecoli-h7.faa -out h7vsk12-$i.tab -evalue $i -outfmt 6
done

• 数据加载 ：创建 MySQL 表，将 BLAST 结果和基因注释文件加载到表中。

-- 创建 h7vsk1200001 表
mysql> CREATE TABLE h7vsk1200001 (quid VARCHAR(15), suid VARCHAR(15), iden FLOAT, alen INT, mism INT, gapo INT, qsta INT, qend INT, ssta INT, send INT, eval FLOAT, bits INT);
-- 加载 BLAST 结果
mysql> LOAD DATA LOCAL INFILE "h7vsk12-0.00001-gi.tab" INTO TABLE h7vsk1200001;
-- 创建 ecolianno 表
mysql> CREATE TABLE ecolianno (gene VARCHAR(12), anno VARCHAR(170), coli VARCHAR(3));
-- 加载基因注释
mysql> LOAD DATA LOCAL INFILE "annotation-k12.tab" INTO TABLE ecolianno;

• 数据查询 ：提取高度保守基因和独特 ORFs。

-- 提取高度保守基因
mysql> SELECT b.quid, b.suid, b.iden, b.alen, e.anno FROM h7vsk1200001 AS b, ecolianno AS e WHERE b.iden = 100 AND b.alen > 800 AND b.suid = e.gene;
-- 提取独特 ORFs
mysql> SELECT count(e.gene) FROM ecolianno AS e LEFT JOIN h7vsk1200001 AS h ON h.suid = e.gene WHERE h.suid IS NULL LIMIT 5;

6.3 R 数据可视化

• 建立连接 ：使用 RMySQL 库连接 MySQL 数据库，获取数据。

> library(RMySQL)
> drv = dbDriver("MySQL")
> con = dbConnect(drv,host="localhost",dbname="compbiol", user="awkologist",pass="awkology")
> data = dbGetQuery(con,statement="SELECT iden, alen, eval FROM h7vsk1200001")

• 绘制图表 ：绘制散点图、直方图和箱线图，进行回归分析和 t 检验。

# 绘制散点图和回归直线
> plot(iden,alen,xlab="Identity [%]",ylab="Alignment Length [bp]")
> abline(lm(alen~iden),col="red")
# 绘制直方图
> hist(iden,xlab="Identity [%]",breaks=70)
# 绘制箱线图
> boxplot(dehydrogenase[[1]],dna[[1]],cytochrome[[1]],notch=T, names=c("DH","DNA-Pol","CYT"),main="Maximum Alignment Length [bp]")
# t 检验
> t.test(dehydrogenase[[1]],cytochrome[[1]])

7. 注意事项

• 文件权限 ：在使用 shell 脚本时，要确保文件具有可执行权限，使用 chmod u+x 命令进行设置。
• MySQL 权限 ：在某些 Linux 系统中，执行 LOAD DATA LOCAL INFILE 命令可能需要额外的权限，可使用 mysqlimport 命令替代。
• R 包安装 ：在使用 R 进行分析前，要确保 DBI 和 RMySQL 包已正确安装。

8. 拓展思考

• 多菌株分析 ：可以将分析扩展到更多的大肠杆菌菌株，比较它们之间的基因差异和相似性，进一步研究致病机制。
• 其他分析方法 ：除了本文介绍的方法，还可以尝试使用其他生物信息学工具和算法，如聚类分析、系统发育分析等，深入挖掘基因数据。
• 数据整合 ：将基因组数据与其他类型的数据（如转录组数据、蛋白质组数据）进行整合分析，以获得更全面的生物学信息。

9. 总结

通过本文的介绍，我们了解了在基因组分析中调整 E 值、使用 MySQL 组织数据以及利用 R 进行可视化和分析的详细方法。这些技术可以帮助我们更好地理解大肠杆菌不同菌株之间的基因特征，为研究致病因子和基因功能提供有力支持。在实际应用中，要根据具体需求灵活运用这些方法，并注意操作过程中的细节和注意事项。同时，不断探索和尝试新的分析思路和方法，以获得更深入的生物学见解。

10. 练习答案提示（仅供参考）

• 问题 1 ：需要从 BLAST 结果中筛选出大肠杆菌 O157:H7 特有的蛋白质序列，并结合注释信息进行统计。
• 问题 2 ：分别对大肠杆菌 O104:H4 与大肠杆菌 K12、O157:H7 进行 BLAST 比对，然后筛选和统计独特的带注释蛋白质。
• 问题 3 ：E 值越高，允许的序列差异越大，因此匹配的结果就越多。
• 问题 4 ：E 值增加，更多的序列会被认为是匹配的，所以被识别为查询菌株特有的开放阅读框数量会减少。
• 问题 5 ：可以使用类似加载大肠杆菌 K12 注释信息的方法，将大肠杆菌 O157:H7 的注释文件转换为合适的格式，然后加载到 ecolianno 表中。
• 问题 6 ：一个蛋白质序列可能与多个序列有比对结果，导致输出文件的行数多于蛋白质序列的数量。
• 问题 7 ：可以使用 R 中的相关函数（如 cor.test ）对脱氢酶和 DNA 聚合酶的最大比对长度数据进行相关性分析。

希望这些内容能帮助你更好地掌握基因组分析的相关知识和技能，通过不断练习和实践，提升自己的生物信息学分析能力。