PRÓPOLIS-UFG
Dayane Cristina da Costa, Leonardo Gonsalves, Lucas de Sousa Rodrigues
2026-04-09
1 Bibliotecas Utilizadas
library(DT)
library(knitr)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(plotly)
library(ggdendro)
library(patchwork)
library(reshape2)
library(tidyverse)
library(readxl)2 Objetivos
Comparar o perfil químico de compostos orgânicos voláteis (VOCs) de amostras de própolis provenientes de diferentes espécies de abelhas, coletadas em uma mesma região, utilizando análises estatísticas univariadas e multivariadas, a fim de identificar compostos discriminantes entre as amostras.
3 Análise Química
3.1 Preparo de Amostra
As amostras foram submetidas à temperatura de -10 °C por aproximadamente 24 horas. Em seguida, realizou-se a maceração com auxílio de gral e pistilo, com o objetivo de facilitar a pesagem. Após esse preparo, foram pesadas em vials de 10 mL aproximadamente 1,0 g em triplicata. Posteriormente, os frascos foram selados com o auxílio de um crimpador automático, utilizando tampas magnéticas.
3.2 Parâmetros Cromatográficos
As amostras foram analisadas via HS/GC-MS em cromatógrafo a gás da Agilent 7890B acoplado ao espectrômetro de massas 7000D. Os VOCs foram extraídos por 15 minutos a 750 rpm em incubadora mantida a 100 °C. O volume de injeção foi de 2,5 mL com seringa mantida a 100 °C. O gás de arraste utilizado foi o hélio 5.0 (99,999%) com velocidade linear de 30,405 cm/s e o injetor operou no modo splitless a 280 °C. Utilizou-se a coluna HP-5ms de 30 m x 250 μm x 0,25 μm. O forno operou a temperatura inicial de 40 °C por 5 minutos, seguido por rampa de 5 °C/min até 160 °C, em seguida, 1 °C/min até 170 °C, por fim, 10 °C/min até 250 °C, totalizando 47 minutos de análise. O espectrômetro de massas operou no modo EI com temperatura da fonte de íons a 300 °C, no modo SCAN a 70 eV com faixa de massa de 17-400 m/z. Os picos foram identificados com base na correlação de similaridade com os espectros padrões da biblioteca NIST 14.
4 Tratamento de Dados
Foi necessário converter os arquivos contendo os cromatogramas de (.D) para (mzML), utilizando o software da Proteo Wizard (Msconvert). Após a conversão, os arquivos foram carregados no software MZmine versão 4.9.14, onde foi realizado todo o processamento de dados com o seguinte workflow:
| Etapa | Descricao |
|---|---|
| 1. Importação | Entrada dos dados brutos |
| 2. Detecção de massas | Centroidização e threshold |
| 3. EIC | Construção de cromatogramas de íons extraídos |
| 4. Deconvolução | Separação de sinais sobrepostos |
| 5. Alinhamento | Comparação entre amostras (tempo de retenção e m/z) |
| 6. Identificação | Busca na biblioteca espectral NIST |
| 7. Exportação | O processamento é salvo em CSV |
4.1 Parâmetros do tratamento de dados no MZmine
A tabela completa contendo os parâmetros do processamento de dados no
MZmine (versão 4.9.14) está disponível no Material Suplementar:
Parâmetros
4.2 Codificação das Amostras
As amostras foram codificadas para facilitar o manuseio dos dados, da seguinte maneira:
| Nome_Científico | Código | Região |
|---|---|---|
| Scaptotrigona póstica | PSP | Câmpus II |
| Tetragonisca angustula | PTA | Câmpus II |
| Melipona fasciculata | PMF | Câmpus II |
| Frieseomelitta varia | PFV | Câmpus II |
4.3 Número de picos cromatográficos identificados
| Amostras | Picos cromatográficos |
|---|---|
| PSP | 31 |
| PTA | 14 |
| PMF | 79 |
| PFV | 80 |
4.4 Identificação dos Compostos
Após o processamento dos dados no software MZmine, foi exportada uma tabela contendo as features detectadas, caracterizadas pela razão massa/carga (m/z), tempo de retenção (RT) e área do pico. A identificação dos compostos foi realizada no próprio software, com base na comparação dos espectros de massas obtidos com aqueles disponíveis na biblioteca espectral NIST, utilizando critérios de similaridade espectral. Posteriormente, foi aplicada uma etapa de filtragem das features, na qual foram removidos os sinais que apareceram apenas uma única vez nas triplicatas de cada amostra, visando aumentar a robustez dos dados e reduzir a influência de possíveis ruídos.
dados <- read.csv2("TABELA_MODULAR_FILTRO.csv")
datatable(dados)4.5 CLASSIFICAÇÃO DOS COMPOSTOS
Foi realizada a contagem de classes químicas por amostra:
dados_class <- read_excel("TABELA_MODULAR_FILTRO_ORGCLASS2.xlsx")
colnames(dados_class)[6] <- "Class_sec"
filtrar_presentes <- function(df, prefixo){
cols <- grep(paste0("^", prefixo), colnames(df), value = TRUE)
df %>%
rowwise() %>%
mutate(
presenca = sum(c_across(all_of(cols)) > 0, na.rm = TRUE)
) %>%
filter(presenca >= 1) %>%
ungroup()}
PTA <- filtrar_presentes(dados_class, "PTA")
PSP <- filtrar_presentes(dados_class, "PSP")
PMF <- filtrar_presentes(dados_class, "PMF")
PFV <- filtrar_presentes(dados_class, "PFV")
set_PTA <- unique(PTA$Class_sec)
set_PSP <- unique(PSP$Class_sec)
set_PMF <- unique(PMF$Class_sec)
set_PFV <- unique(PFV$Class_sec)
count_PTA <- PTA %>% count(Class_sec)
count_PSP <- PSP %>% count(Class_sec)
count_PMF <- PMF %>% count(Class_sec)
count_PFV <- PFV %>% count(Class_sec)4.5.1 PTA
datatable(count_PTA)4.5.2 PVF
datatable(count_PFV)4.5.3 PSP
datatable(count_PSP)4.5.4 PMF
datatable(count_PMF)5 Análise estatística
5.1 NORMALIZAÇÃO
Os dados foram normalizados por meio da padronização (z-score), na qual se subtrai a média da intensidade de cada composto e divide-se pelo desvio padrão. Essa normalização é necessária porque compostos de baixa intensidade podem ser relevantes; no entanto, na presença de compostos mais intensos, os valores elevados tendem a dominar as análises estatísticas, suprimindo a contribuição dos compostos menos intensos.
dados_matriz <- read.csv2("MATRIZ-PROPOLIZ-UFG.csv")
dados_numericos <- dados_matriz[, -c(1, 2)]
dados_normalizados <- as.data.frame(scale(dados_numericos))5.2 ANOVA + teste post hoc
Após a normalização dos dados por meio da padronização (z-score), foi realizada a análise de variância (ANOVA), com o objetivo de avaliar diferenças na abundância dos compostos entre as diferentes classes de amostras. A ANOVA compara as médias entre grupos a partir da razão entre a variabilidade entre grupos e a variabilidade dentro dos grupos, resultando em um valor de p. Quando p < 0,05, rejeita-se a hipótese nula, indicando diferença estatisticamente significativa entre pelo menos dois grupos.
Após a ANOVA, foi aplicada uma correção para múltiplas comparações utilizando o método FDR (False Discovery Rate), com o objetivo de controlar a taxa de falsos positivos. Foram considerados significativos os compostos com valores de FDR < 0,05.
Adicionalmente, para cada composto significativo, foi identificada a classe que apresentou maior abundância média, permitindo a associação entre compostos discriminantes e grupos específicos de amostras.
5.3 HEATMAP
O heatmap foi utilizado para visualizar a distribuição da abundância relativa dos compostos entre as diferentes classes de própolis, com base nos dados normalizados por z-score. De modo geral, observou-se variabilidade na composição química entre as amostras, evidenciada por diferenças nos padrões de intensidade dos compostos.
A classe PSP apresentou menor número de compostos com alta intensidade relativa, entretanto, os compostos mais abundantes pertencem majoritariamente a uma mesma classe química, que são os álcoois, sugerindo um perfil mais específico. Por outro lado, a classe PTA não apresentou compostos com elevada intensidade relativa, indicando um perfil químico mais homogêneo e menos complexo.
As classes PMF e PFV exibiram padrões de abundância semelhantes, com regiões do heatmap nas quais os compostos apresentam intensidades próximas entre essas classes, além de regiões com distinção evidente. Esse comportamento é consistente com a análise hierárquica, que indicou um agrupamento parcial entre essas amostras, sugerindo similaridade no perfil químico.
# Preparar matriz
mat <- as.matrix(dados_normalizados)
rownames(mat) <- dados_matriz$Sample
colnames(mat) <- colnames(dados_normalizados)
# Clusterização hierárquica (HCA)
hc_col <- hclust(dist(t(mat)), method = "ward.D2") # compostos
hc_row <- hclust(dist(mat), method = "ward.D2") # amostras
ordem_col <- hc_col$order
ordem_row <- hc_row$order
# Reorganizar dados
mat_ord <- mat[ordem_row, ordem_col]
# Converter para formato longo
df_melt <- reshape2::melt(mat_ord)
colnames(df_melt) <- c("Sample", "Compound", "value")
df_melt$Sample <- factor(df_melt$Sample, levels = rownames(mat_ord))
df_melt$Compound <- factor(df_melt$Compound, levels = colnames(mat_ord))
# Limite da escala
limite <- max(abs(df_melt$value), na.rm = TRUE)
# Dendrograma (colunas)
dend_col <- ggdendro::dendro_data(as.dendrogram(hc_col), type = "rectangle")
p_dend_top <- ggplot2::ggplot(ggdendro::segment(dend_col)) +
ggplot2::geom_segment(aes(x = x, y = y, xend = xend, yend = yend),
linewidth = 0.4, color = "grey20") +
ggplot2::theme_void()
# Dendrograma (linhas)
dend_row <- ggdendro::dendro_data(as.dendrogram(hc_row), type = "rectangle")
p_dend_left <- ggplot2::ggplot(ggdendro::segment(dend_row)) +
ggplot2::geom_segment(aes(x = x, y = y, xend = xend, yend = yend),
linewidth = 0.4, color = "grey20") +
ggplot2::coord_flip() +
ggplot2::scale_y_reverse() +
ggplot2::theme_void()
# Heatmap
p_heat <- ggplot2::ggplot(df_melt, aes(x = Compound, y = Sample, fill = value)) +
ggplot2::geom_tile(color = "white", linewidth = 0.1) +
ggplot2::scale_fill_gradient2(
low = "#2166AC",
mid = "white",
high = "#B2182B",
midpoint = 0,
limits = c(-limite, limite),
name = "z-score"
) +
ggplot2::theme_minimal(base_size = 9) +
ggplot2::theme(
axis.text.x = element_text(angle = 90, hjust = 1, size = 6),
axis.text.y = element_text(size = 8),
axis.title = element_blank(),
panel.grid = element_blank()
)
# Barra de classe
cores_classe <- c(PSP = "#E6821E", PTA = "#2E8B57", PMF = "#7B52AB", PFV = "#D63A6B")
df_class <- data.frame(
Sample = factor(rownames(mat_ord), levels = rownames(mat_ord)),
Class = dados_matriz$Class[ordem_row]
)
p_class <- ggplot2::ggplot(df_class, aes(x = 1, y = Sample, fill = Class)) +
ggplot2::geom_tile() +
ggplot2::scale_fill_manual(values = cores_classe) +
ggplot2::theme_void()
# Combinar tudo
layout <- "
AABBB
CDEEE
"
p_final <- patchwork::wrap_plots(
A = patchwork::plot_spacer(),
B = p_dend_top,
C = patchwork::plot_spacer(),
D = p_class,
E = p_heat,
design = layout,
widths = c(0.05, 0.03, 1),
heights = c(0.15, 1)
)
p_final
5.4 PCA
Foi realizada a análise de componentes principais (PCA), que é uma análise multivariada que reduz a dimensionalidade dos dados, permitindo a visualização da distribuição das amostras em função dos componentes principais (PCs), os quais explicam a variância total dos dados. Observou-se que as amostras das classes PTA e PSP apresentaram maior proximidade, sugerindo similaridade em seus perfis químicos. Por outro lado, as classes PMF e PFV apresentaram maior separação ao longo dos componentes principais, indicando diferenças no perfil químico entre essas amostras.
Adicionalmente, os compostos com maiores cargas (loadings) nos componentes principais são os principais responsáveis pela separação observada entre as amostras.
# PCA
pca_res <- prcomp(dados_normalizados, scale. = FALSE)
# Scores (amostras)
scores <- as.data.frame(pca_res$x)
scores$Class <- Class
scores$Sample <- Sample
# Variância explicada
var_exp <- (pca_res$sdev^2 / sum(pca_res$sdev^2)) * 100
# PCA 3D
plot_ly(
data = scores,
x = ~PC1,
y = ~PC2,
z = ~PC3,
color = ~Class,
colors = c("PSP" = "#E07B8A", "PMF" = "#6BAE75", "PTA" = "black", "PFV" = "red"),
text = ~Sample,
type = "scatter3d",
mode = "markers",
marker = list(size = 10)
) %>%
layout(
title = "PCA",
scene = list(
xaxis = list(title = paste0("PC1 (", round(var_exp[1],1), "%)")),
yaxis = list(title = paste0("PC2 (", round(var_exp[2],1), "%)")),
zaxis = list(title = paste0("PC3 (", round(var_exp[3],1), "%)"))))5.5 BOXPLOT
5.5.0.1 COMP. COMPOSTOS SIGNIFICATIVOS
A partir dos resultados obtidos pela ANOVA, seguida de correção para múltiplas comparações (FDR), foram selecionados os três compostos mais significativos para cada classe, com base nos menores valores de p ajustados. A distribuição das intensidades normalizadas desses compostos foi avaliada por meio de boxplots, permitindo comparar sua variação entre as diferentes classes de própolis. Os resultados indicaram uma associação entre os compostos discriminantes e classes químicas específicas. As amostras da classe PSP apresentaram maior abundância relativa de compostos classificados como aldeídos, enquanto a classe PMF foi caracterizada principalmente por ésteres, e a classe PFV por terpenos.
gerar_boxplot_classe <- function(classe_alvo, top_n = 3) {
comps <- resultado_anova %>%
filter(FDR < 0.05,
Classe_Maior_Abundancia == classe_alvo) %>%
arrange(p_value) %>%
slice(1:top_n) %>%
pull(Composto)
if(length(comps) == 0){
message(paste("Sem compostos para:", classe_alvo))
return(NULL)
}
dados_plot <- dados_normalizados[, comps, drop = FALSE] %>%
as.data.frame() %>%
mutate(Class = Class) %>%
pivot_longer(
cols = -Class,
names_to = "Composto",
values_to = "Intensidade"
)
cores <- c(
"PSP" = "grey70",
"PMF" = "grey70",
"PFV" = "grey70",
"PTA" = "grey70"
)
cores[classe_alvo] <- "#E64B35"
ggplot(dados_plot, aes(x = Composto, y = Intensidade, fill = Class)) +
geom_boxplot(alpha = 0.85, outlier.shape = 21) +
scale_fill_manual(values = cores) +
labs(
title = paste("Compostos com maior abundância em", classe_alvo),
x = "Composto",
y = "Intensidade normalizada"
) +
theme_classic(base_size = 13) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "top"
)
}
p_PSP <- gerar_boxplot_classe("PSP", top_n = 3)
p_PMF <- gerar_boxplot_classe("PMF", top_n = 3)
p_PFV <- gerar_boxplot_classe("PFV", top_n = 3)
p_PSP
p_PMF
p_PFV
5.6 HISTOGRAMA
Foi construído um histograma sobreposto com o objetivo de visualizar a distribuição da intensidade dos compostos ao longo do tempo de retenção (RT) para as diferentes classes de amostras. Os resultados indicaram diferenças na distribuição das intensidades ao longo do cromatograma entre as classes, evidenciando variações no perfil químico global das amostras. Observou-se que determinadas regiões de RT apresentaram maior contribuição de intensidade para classes específicas, sugerindo a presença de compostos característicos em faixas distintas do cromatograma. Essa abordagem permite uma visualização global do perfil cromatográfico.
dados_long <- dados %>%
pivot_longer(
cols = matches("^(PFV|PMF|PSP|PTA)"),
names_to = "amostra",
values_to = "intensidade"
) %>%
filter(intensidade > 0) %>%
mutate(grupo = stringr::str_extract(amostra, "^[A-Z]+"))
bin_width <- 0.5
dados_bin <- dados_long %>%
mutate(bin = floor(rt / bin_width) * bin_width) %>%
group_by(grupo, bin) %>%
summarise(intensidade = sum(intensidade), .groups = "drop")
cores <- c(
"PFV" = "#534AB7",
"PMF" = "#1D9E75",
"PSP" = "#D85A30",
"PTA" = "#D4537E")
ggplot() +
geom_col(
data = dados_bin,
aes(x = bin, y = intensidade, fill = grupo),
position = "identity",
alpha = 0.6) +
geom_line(
data = dados_bin,
aes(x = bin, y = intensidade, color = grupo),
linewidth = 1.2,
alpha = 0.0) +
scale_fill_manual(values = cores) +
scale_color_manual(values = cores) +
scale_y_continuous(expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
labs(
title = "Distribuição de RT por grupo",
x = "RT (min)",
y = "Intensidade") +
theme_classic(base_size = 13) +
theme(
legend.position = "top",
plot.title = element_text(size = 14, face = "bold", hjust = 0.5))