31、R语言统计分析与基因组比较：从数据探索到实践应用

R语言统计分析与基因组比较：从数据探索到实践应用

1. R语言中的数据探索

在数据分析领域，数据探索是至关重要的一步。以Anscombe四重奏为例，它提醒我们不能仅仅依赖简单的回归分析。虽然所有数据集的相关系数均为0.816，且线性回归线都是(y = 3.0 + 0.5x)，但通过绘图可以发现数据背后隐藏的不同模式。

以下是绘制Anscombe四重奏的R代码：

> anscombe<-read.table(
+ "http://www.hs-mittweida.de/wuenschi/data/media/
+ compbiolbook/anscombe.tab",header=TRUE,sep="\t")
> attach(anscombe)
> par(mfrow=c(2,2))
> plot(ax1,ay1,xlim=c(4,18),ylim=c(3,12),col="blue")
> abline(lm(ay1~ax1))
> plot(ax2,ay2,xlim=c(4,18),ylim=c(3,12),col="blue")
> abline(lm(ay2~ax2))
> plot(ax3,ay3,xlim=c(4,18),ylim=c(3,12),col="blue")
> abline(lm(ay3~ax3))
> plot(ax4,ay4,xlim=c(4,18),ylim=c(3,12),col="blue")
> abline(lm(ay4~ax4))
> par(mfrow=c(1,1))

在这段代码中，我们首先从互联网下载数据，然后使用 par() 函数将图形设置为2x2的布局，以便在一个图中显示四个子图。 xlim() 和 ylim() 参数用于设置坐标轴的范围。

2. 统计检验方法

2.1 t检验

在生命科学中，t检验是一种非常重要的统计检验方法。它是一种参数检验，假设数据呈正态分布且样本方差相似，并且数据集之间相互独立。t检验用于比较两个数据集，更具体地说，是寻找两个样本之间的差异。

t检验的原假设（(H_0)）是两个样本的均值没有差异。t检验会提供一个p值，它表示错误拒绝(H_0)的概率。通常，当p值小于或等于0.05时，我们会拒绝(H_0)，此时我们有95%的信心认为正确地拒绝了原假设。

以下是使用R进行t检验的示例，我们使用两个不同大肠杆菌菌株中所有基因的GC含量数据集：

> allk12 <- read.table("all-k12-id-gc.tab", header=FALSE, sep="\t")
> allo157 <- read.table("all-o157-id-gc.tab", header=FALSE, sep="\t")
> hist(ecolik12[[2]],main="E. coli K12",xlab="GC-Content [%]")
> hist(ecolio157[[2]],main="E. coli O157:H7",xlab="GC-Content [%]")
> sd(ecolik12[[2]])
[1] 5.125811
> sd(ecolio157[[2]])
[1] 6.01197
> t.test(allk12[[2]],allo157[[2]])

运行结果如下：

Welch Two Sample t-test

data:  allk12[[2]] and allo157[[2]]
t = 5.4471, df = 9983.783, p-value = 5.241e-08
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.3881420 0.8245378
sample estimates:
mean of x mean of y 
  50.67152   50.06518 

从结果可以看出，p值远低于0.05，因此我们可以拒绝原假设，即这两个大肠杆菌菌株中基因的GC含量均值存在显著差异。

2.2 卡方检验

卡方检验是生命科学中另一种重要的统计检验方法。它是一种非参数检验，用于评估观察频率和预期频率之间的差异是由于随机因素还是有其他有意义的原因。卡方检验只能用于频率数据，而不能用于百分比、尺度数据等。

原假设（(H_0)）通常表示两个数据总体是同质的，即频率没有差异。与t检验一样，卡方检验也会返回一个p值，用于判断错误拒绝原假设的概率。当p值小于或等于0.05时，我们拒绝原假设。

以下是使用卡方检验比较两个不同基因集的核苷酸频率的示例，数据集如下表所示：
| | 基因组主干 | 基因组岛 |
| — | — | — |
| 腺嘌呤 | 994360 | 702 |
| 胸腺嘧啶 | 988925 | 578 |
| 鸟嘌呤 | 1132739 | 936 |
| 胞嘧啶 | 1010405 | 322 |

以下是在R中进行卡方检验的代码：

> matrix(c(994360,988925,1132739,1010405,702,578,936,322), ncol=2)
     [,1]  [,2]
[1,] 994360  702
[2,] 988925  578
[3,] 1132739 936
[4,] 1010405 322
> chisq.test(matrix(c(994360,988925,1132739,1010405,702,578,936,322), ncol=2))

运行结果如下：

Pearson’s Chi-squared test

data:  matrix(c(994360, 988925, 1132739, 1010405,
 702, 578, 936, 322),ncol = 2)
X-squared = 241.2361, df = 3, p-value < 2.2e-16

由于p值远低于0.05，我们可以拒绝原假设，即基因组主干和基因组岛的核苷酸组成存在显著差异，这可能是水平基因转移的迹象。

3. R脚本

R是一种编程语言，我们可以将命令写入文件，加载到R中并执行。以下是一个简单的R脚本示例：

# save as import-export.r
data<-read.table("gene-exprA-exprB.tab",header=TRUE,sep="\t")
attach(data)
png("mioploto.png")
plot(log2(exprA),log2(exprB),col="blue")
dev.off()

在Bash命令行中执行以下命令可以运行该脚本：

R --slave --vanilla < script.r > /dev/null

这个命令会将存储在 script.r 中的R命令重定向到R，并将任何输出重定向到空设备 /dev/null 。 --slave 和 --vanilla 选项使R在静默模式下工作。

4. 扩展：安装和使用包

R的功能可以通过包来扩展。截至2012年9月，综合R存档网络（CRAN）包存储库中有4049个可用包。以RMySQL包为例，它可以让我们从R连接到MySQL数据库。

4.1 连接到MySQL

要连接到MySQL服务器，需要加载特定的库。如果在R命令 shell中执行 library("RMySQL") 没有错误消息，则表示该库已经安装。否则，需要在R命令行中执行以下命令来安装该包：

install.packages("RMySQL")

在安装过程中，需要指定一个CRAN镜像服务器。如果安装失败，可能是缺少一些MySQL库和头文件，需要以管理员权限在Bash命令行中执行以下命令来安装这些文件：

sudo apt-get install libmysqlclient-dev

安装成功后，再次执行 install.packages("RMySQL") 。

以下是连接到MySQL数据库的示例代码：

> library("RMySQL")
Loading required package: DBI
> con <- dbConnect(MySQL(), user="awkologist", password="awkology",
+ dbname="compbiol", host="localhost")
> dbListTables(con)
[1] "annotation" "expression"
> dbGetQuery(con, "SELECT * FROM expression LIMIT 3")
  gene expr_value
1 alr1207 1000
2 alr2938 10323
3 alr3395 1432
> expression<-dbGetQuery(con, "SELECT * FROM expression")
> names(expression)
[1] "gene"      "expr_value"
> attach(expression)
> mean(expr_value)
[1] 3984.571
> median(expr_value)
[1] 1432

4.2 生命科学相关包

除了RMySQL，还有许多其他可能感兴趣的R包，例如：
- nlstools ：用于细菌生长和酶动力学分析
- ape ：用于系统发育和进化分析
- seqRFLP ：用于模拟和可视化DNA序列的限制性酶切模式
- HardyWeinberg ：基于三元图进行Hardy-Weinberg平衡的图形测试
- Bioconductor ：用于DNA微阵列数据分析
- seqinR ：用于生物序列检索和分析

这些包可以从CRAN下载并使用 install.packages("...") 函数进行安装。

5. 基因组比较项目

5.1 项目背景

大肠杆菌有许多不同的菌株，其中一些是致病的，一些是非致病的。以大肠杆菌O157:H7为例，它是食源性疾病的新兴原因，感染后可能导致血性腹泻、腹痛等症状，甚至可能发展为肾衰竭。而在2011年，另一种致病菌株大肠杆菌O104:H4在德国和欧洲引起了关注，它比O157:H7更具毒性。

我们的目标是通过计算比较非致病菌株大肠杆菌K12和致病菌株大肠杆菌O157:H7以及大肠杆菌O104:H4的翻译注释基因组，找出可能导致致病性的蛋白质。

5.2 生物信息学工具和资源

• BLAST+ ：是一种标准程序，用于查找序列或序列集之间的相似性。可以从美国国家生物技术信息中心（NCBI）免费下载，通常人们通过网页界面运行该程序，这里我们将从命令行下载、安装和运行它。
• 基因组数据库 ：我们将使用NCBI基因组数据库和Pathosystems资源整合中心（PATRIC）。NCBI基因组数据库包含非致病实验室菌株大肠杆菌K12亚种MG1655和致病菌株大肠杆菌O157:H7亚种EDL933的基因组序列。PATRIC可用于下载大肠杆菌O104:H4亚种TY-2482的注释基因组序列。

5.3 详细步骤

5.3.1 下载和安装BLAST+

• 使用网页浏览器下载 ：打开网页浏览器，访问http://www.ncbi.nlm.nih.gov，点击“Downloads: Get NCBI data or software”，选择“BLAST (Stand-alone)”，找到最新版本（如果不是2.2.25，则进入2.2.25目录），下载适合操作系统的压缩存档文件（如 ncbi-blast-2.2.25+-ia32-linux.tar.gz 或 ncbi-blast-2.2.25+-universal-macosx.tar.gz ）到主目录。
• 安装 ：假设已经下载了 ncbi-blast-2.2.25+-ia32-linux.tar.gz 文件，在Bash命令行中执行以下命令进行解压：

$ gunzip ncbi-blast-2.2.25+-ia32-linux.tar.gz
$ tar -xf ncbi-blast-2.2.25+-ia32-linux.tar
$ cd ncbi-blast-2.2.25+
$ ls
bin ChangeLog doc LICENSE ncbi_package_info README
$ ls bin
blastdb_aliastool blastp makeblastdb segmasker
blastdbcheck blastx makembindex tblastn
blastdbcmd convert2blastmask psiblast tblastx
blast_formatter dustmasker rpsblast update_blastdb.pl
blastn legacy_blast.pl rpstblastn windowmasker
$ ls -l bin
...
$ rm ../ncbi-blast-2.2.25+-ia32-linux.tar

如果磁盘空间有限，可以只提取必要的程序：

tar -xf ncbi-blast-2.2.25+-ia32-linux.tar
+ ncbi-blast-2.2.25+/bin/makeblastdb +
ncbi-blast-2.2.25+/bin/blastp
+ ncbi-blast-2.2.25+/bin/blastdbcmd

通过以上步骤，我们可以完成R语言中的数据探索、统计检验、脚本编写、包的安装和使用，以及基因组比较项目的准备工作。这些知识和技能在生物信息学和数据分析领域具有重要的应用价值。

R语言统计分析与基因组比较：从数据探索到实践应用

6. 基因组比较的后续步骤

在完成BLAST+的下载和安装后，我们可以继续进行基因组比较的后续操作。

6.1 数据准备

我们已经从NCBI基因组数据库和PATRIC获取了大肠杆菌K12、O157:H7和O104:H4的基因组序列数据。接下来，需要将这些数据准备成适合BLAST+分析的格式。

通常，BLAST+需要一个数据库和一个查询序列。我们可以将其中一个基因组作为数据库，另一个作为查询序列。例如，将大肠杆菌K12的基因组作为数据库，大肠杆菌O157:H7的基因组作为查询序列。

6.2 创建BLAST数据库

使用 makeblastdb 程序可以将基因组序列文件转换为BLAST数据库。以下是一个示例命令：

$ makeblastdb -in ecoli_k12.fasta -dbtype nucl -out ecoli_k12_db

• -in ：指定输入的基因组序列文件，这里是 ecoli_k12.fasta 。
• -dbtype ：指定数据库类型， nucl 表示核苷酸数据库。
• -out ：指定输出的数据库名称，这里是 ecoli_k12_db 。

6.3 运行BLAST搜索

使用 blastn 程序可以在创建好的数据库中搜索查询序列。以下是一个示例命令：

$ blastn -query ecoli_o157.fasta -db ecoli_k12_db -out results.txt -outfmt 6

• -query ：指定查询序列文件，这里是 ecoli_o157.fasta 。
• -db ：指定要搜索的数据库，这里是 ecoli_k12_db 。
• -out ：指定输出结果文件，这里是 results.txt 。
• -outfmt ：指定输出格式， 6 表示以制表符分隔的表格格式输出。

6.4 结果分析

运行BLAST搜索后，会得到一个结果文件 results.txt 。可以使用R或其他工具对结果进行分析。以下是一个简单的R代码示例，用于读取和分析BLAST结果：

results <- read.table("results.txt", header = FALSE)
colnames(results) <- c("query_id", "subject_id", "percent_identity", "alignment_length", "mismatches", "gap_opens", "q_start", "q_end", "s_start", "s_end", "evalue", "bit_score")

# 筛选出相似度较高的匹配
high_similarity <- results[results$percent_identity > 90, ]

# 统计匹配的数量
match_count <- nrow(high_similarity)

# 输出结果
print(paste("相似度大于90%的匹配数量:", match_count))

7. 寻找独特的开放阅读框（ORFs）

在进行基因组比较后，我们的目标是找到在一个基因组中独特的开放阅读框（ORFs）。这些独特的ORFs可能与致病性相关。

7.1 预测ORFs

可以使用一些生物信息学工具来预测基因组中的ORFs，例如Prodigal。以下是一个使用Prodigal预测ORFs的示例命令：

$ prodigal -i ecoli_o157.fasta -o ecoli_o157_orfs.gff -a ecoli_o157_orfs.faa

• -i ：指定输入的基因组序列文件，这里是 ecoli_o157.fasta 。
• -o ：指定输出的GFF格式的ORF注释文件，这里是 ecoli_o157_orfs.gff 。
• -a ：指定输出的氨基酸序列文件，这里是 ecoli_o157_orfs.faa 。

7.2 筛选独特的ORFs

在得到ORF预测结果后，可以将这些ORFs与其他基因组进行比较，筛选出独特的ORFs。可以使用BLAST或其他序列比对工具进行比较。以下是一个简单的流程：
1. 将预测的ORFs作为查询序列，与其他基因组的数据库进行BLAST搜索。
2. 筛选出没有显著匹配的ORFs，这些ORFs就是独特的ORFs。

以下是一个使用R和BLAST结果筛选独特ORFs的示例代码：

# 读取BLAST结果
blast_results <- read.table("blast_results.txt", header = FALSE)
colnames(blast_results) <- c("query_id", "subject_id", "percent_identity", "alignment_length", "mismatches", "gap_opens", "q_start", "q_end", "s_start", "s_end", "evalue", "bit_score")

# 读取ORF列表
orfs <- read.table("ecoli_o157_orfs.faa", header = FALSE)
colnames(orfs) <- c("orf_id", "sequence")

# 筛选出没有显著匹配的ORFs
unique_orfs <- orfs[!(orfs$orf_id %in% blast_results$query_id), ]

# 输出独特的ORFs
write.table(unique_orfs, "unique_orfs.txt", sep = "\t", quote = FALSE, row.names = FALSE)

8. 总结与展望

通过以上步骤，我们完成了从数据探索、统计检验、R脚本编写、包的安装和使用，到基因组比较和寻找独特ORFs的整个过程。这些操作在生物信息学和基因组学研究中具有重要的应用价值。

在未来的研究中，我们可以进一步深入分析这些独特的ORFs，例如进行功能注释、结构预测等，以更好地理解它们在致病性中的作用。同时，还可以结合其他数据和方法，如转录组学、蛋白质组学等，进行更全面的研究。

以下是整个基因组比较流程的mermaid格式流程图：

graph LR
    A[数据准备] --> B[创建BLAST数据库]
    B --> C[运行BLAST搜索]
    C --> D[结果分析]
    D --> E[预测ORFs]
    E --> F[筛选独特ORFs]

通过这个流程图，我们可以清晰地看到整个基因组比较的流程，从数据准备到最终筛选出独特的ORFs。

总之，R语言和生物信息学工具为我们提供了强大的数据分析和处理能力，帮助我们在基因组学研究中取得更深入的成果。希望本文介绍的方法和步骤能够对读者在相关领域的研究有所帮助。