2024-06-16

让你的数据“开口说话”！从统计到可视化，把热图、PCA与箱线图一网打尽！

原创小果生信果 2024-05-25 19:00:54

收录于话题

大家好，我是小果~在生物科学的海洋中，基因表达数据如同复杂的密码，等待我们去解读。今天小果将带小伙伴们一起探索基因数据集处理过程，让数据可视化成为我们的“翻译官”！从直观的热图展示基因活跃度，到主成分分析（PCA）揭示数据的深层结构，再到箱线图揭示基因表达值的分布规律，每一步都让我们更深入地理解基因世界的奥秘。无需复杂的数学公式，只需掌握这些可视化工具，我们就能让数据“开口说话”，轻松解读基因表达背后的故事。快来跟随我们的步伐，一起探索基因表达的无限可能吧！同时呢，十年的生信之路，小果已经练就了一身扎实的本领，现在准备用这身本事为小伙伴们服务啦！如果你在生信分析上遇到了难题，那就来找小果吧！小果会用自己的专业知识和技能，为你解决困扰，助你一臂之力。期待你的联系哦~

一、用到的R包简要介绍

lBiobase 包：

它是Bioconductor 项目的核心组成部分，定义了如ExpressionSet这样的关键数据类（存储和操作基因表达数据的标准方式）。Biobase 包含多种功能，它支持数据查询、替换、聚合等操作，并且提供了丰富的方法和辅助函数。此外，Biobase 包还提供了数据操作和包管理的函数，如创建数据包和版本控制。

lMASS包：

它可以进行线性判别分析和二次判别分析（根据一组观测值的特征将其分为不同的类别）。可以使用lda()和qda()函数来进行这两种判别分析。此外，虽然MASS包本身并不直接提供神经网络建模的功能，但可以结合其他R包（如neuralnet包）和MASS包中的数据集来构建神经网络模型。MASS包还提供了一些用于数据可视化的函数，如parcoord()函数，它可以创建平行坐标图来展示多元数据。并且MASS包还提供了许多其他统计分析方法，如回归分析、方差分析等。

ltidyr包：

它是我们非常熟悉的数据处理包啦！主要用于数据的整理和清洗。tidyr包核心功能包括数据的长宽转换、处理数据框中的空值、根据一个表格衍生出其他表格，以及实现行或列的分隔和合并。

二、R语言代码实操

废话不多说，小果直接给小伙伴们上代码！本次介绍的R包操作会占用内存较大，小果建议使用服务器，欢迎小伙伴们联系小果租赁性价比居高的服务器！

公众号后台回复“111”

领取本篇代码、基因集或示例数据等文件

文件编号：240520-2

需要租赁服务器的小伙伴可以扫码添加小果，此外小果还提供生信分析，思路设计，文献复现等，有需要的小伙伴欢迎来撩~

Ø下载并加载Biobase包

if (!requireNamespace("Biobase", quietly = TRUE)) {   install.packages("Biobase")  } library(Biobase)

Ø生成基因表达数据：小果先创建了一个名为genes的向量，包含基因的名称，以及samples向量。然后，使用runif函数生成了一个随机数矩阵expression_data，模拟了基因在样本中的表达水平。矩阵的每一行代表一个基因，每一列代表一个样本。

genes <- c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5") samples <- c("Sample1", "Sample2", "Sample3") set.seed(123)  expression_data <- matrix(runif(15, min = 1, max = 10), nrow = 5, ncol = 3,                            dimnames = list(genes, samples))

Ø计算基因统计信息：小果使用了apply函数对expression_data矩阵的每一行（即每个基因）应用一个函数，计算该行的平均值和标准差，并将结果存储在gene_stats列表中。然后，使用DataFrame函数将基因名称和统计信息组合成一个数据框gene_info。

gene_stats <- apply(expression_data, 1, function(x) c(mean = mean(x), sd = sd(x))) gene_info <- DataFrame(Gene = genes, t(gene_stats)) colnames(gene_info)[-1:-2] <- c("mean", "sd")head(gene_info)

Ø安装和加载MASS包：检查是否已经安装了MASS包。

if (!requireNamespace("MASS", quietly = TRUE))  install.packages("MASS")library(MASS)

Ø绘制基因表达散点图矩阵：使用plot函数绘制基因表达数据的散点图矩阵。

plot(expression_data, main = "Gene Expression Scatterplot Matrix")

Ø数据标准化：对原始的表达数据expression_data进行转置和标准化处理，生成normalized_data。（标准化是通过减去每行的均值并除以标准差来实现的哦~）。然后，计算标准化后数据的均值和标准差，并将这些统计

信息添加到gene_info数据框中。normalized_data <- t(scale(t(expression_data)))gene_info <- cbind(gene_info, DataFrame(normalized_mean = rowMeans(normalized_data), normalized_sd = apply(normalized_data, 1, sd)))colnames(gene_info)[5:6] <- c("normalized_mean", "normalized_sd")head(gene_info)

Ø计算相关性矩阵：小果这里使用了cor函数计算expression_data的成对相关性，生成cor_matrix。（这里使用的是皮尔逊相关系数）

cor_matrix <- cor(expression_data, method = "pearson")head(cor_matrix)

Ø绘制相关性热图：小果使用了reshape2::melt函数将相关性矩阵cor_matrix转换为长格式，然后，使用了ggplot2包绘制相关性热图，其中颜色深浅表示相关性的强度。

cor_melt <- reshape2::melt(cor_matrix)ggplot(cor_melt, aes(Var1, Var2, fill = value)) +  geom_tile() +  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +  theme_minimal() +  theme(axis.text.x = element_text(angle = 45, hjust = 1),        axis.title = element_blank()) +  labs(x = "Gene", y = "Gene")

Ø安装和加载ggplot2包if (!requireNamespace("ggplot2", quietly = TRUE))  install.packages("ggplot2")library(ggplot2)

Ø转换数据框：将gene_info数据框转换为普通的数据框（as.data.frame），以确保数据类型正确。

gene_info <- as.data.frame(gene_info)

Ø数据重构：使用tidyr::pivot_longer函数将gene_info数据框从宽格式转换为长格式（为了更灵活的数据分析和可视化）。

long_data <- tidyr::pivot_longer(gene_info, cols = -Gene, names_to = "Sample", values_to = "Expression")

Ø绘制热图：使用ggplot2绘制长格式数据的热图。

ggplot(long_data, aes(x = Sample, y = Gene, fill = Expression)) +  geom_tile() +  theme_minimal() +  theme(axis.text.x = element_text(angle = 45, hjust = 1),        axis.title.y = element_blank(),        panel.border = element_blank()) +  scale_fill_gradient(low = "blue", high = "red")

Ø主成分分析（PCA）：小果使用了prcomp函数对转置的表达数据进行主成分分析，并对数据进行标准化处理。

pca_result <- prcomp(t(expression_data), scale. = TRUE)

ØPCA结果总结：使用summary函数查看PCA结果的摘要，包括每个主成分的方差比等。并且使用biplot函数绘制PCA的双标图，显示样本和变量之间的关系。

summary(pca_result)biplot(pca_result, scale = TRUE, col = c("red", "green", "blue"))

Ø安装和加载tidyr包

if (!requireNamespace("tidyr", quietly = TRUE))  install.packages("tidyr")library(tidyr)

Ø再次数据重构：使用pivot_longer函数将gene_info数据框的统计数据转换为长格式，以便于后续的可视化分析。

long_gene_info <- pivot_longer(gene_info, cols = -Gene, names_to = "Statistic", values_to = "Expression")

Ø绘制箱线图、散点图等进行可视化：使用ggplot2绘制散点图，展示不同基因的表达值。小果这里使用了facet_wrap函数来为每个基因创建一个单独的面板。

ggplot(long_gene_info, aes(x = Gene, y = Expression, fill = Gene)) +  geom_boxplot() +  theme_minimal() +  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  facet_wrap(~Gene, scales = "free_y") +  labs(x = "gene", y = "Expression Value")

ggplot(long_gene_info, aes(x = Gene, y = Expression, color = Gene)) +  geom_point(size = 2) +  theme_minimal() +  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  facet_wrap(~Gene ,scales = "free_y") +  labs(x = "Gene", y = "Expression Value")

long_gene_info$Gene <- as.numeric(as.factor(long_gene_info$Gene))ggplot(long_gene_info, aes(x = Gene, y = Expression, color = Gene)) +  geom_line() +  geom_point() +  theme_minimal() +  labs(x = "Time Point", y = "Expression Value")

三、文章小结

本期小果和小伙伴们一起学习了从数据生成、统计分析到数据可视化的完整流程。通过这些步骤，我们可以探索基因表达数据中的模式和关联性，为进一步的生物学研究提供相关依据。并且，掌握热图、PCA和箱线图，就是掌握了基因表达数据的钥匙！它们不仅能揭示基因间的微妙关系，还能让我们一眼看穿数据的内在规律。别再让数据沉默啦，将复杂的基因表达数据转化为直观、易懂的图表吧。不论是寻找研究突破口，还是展示研究成果，这些可视化工具都能让你事半功倍。现在就开始探索吧，让基因数据“开口说话”，为你的科研之路点亮明灯！

希望本期的内容对你有帮助，期待我们下期再见！最后如果各位小伙伴们觉得自己运行代码太麻烦，欢迎用我们的云生信小工具，只要输入合适的数据就可以直接出想要的图呢，附云生信链接

（http://www.biocloudservice.com/home.html）。

小果还提供思路设计、定制生信分析、文献思路复现；有需要的小伙伴欢迎直接扫码咨询小果，竭诚为您的科研助力！