Projeto de Transcriptômica: FASD e Dieta Hiperlipídica em Zebrafish

Definição da Pergunta Biológica

Objetivo Central: Avaliar como a Exposição Alcoólica Fetal (FAE) altera o transcriptoma basal de embriões de Danio rerio aos 13 dpf e investigar como essas perturbações se refletem na topologia de redes de interação gênica em condições de dieta normal (ND).

O fluxo de análise buscará responder: 1. Quais genes são diferencialmente expressos pela FAE sob dieta normal? 2. Quais processos biológicos estão enriquecidos nessa assinatura transcricional? 3. Como a conectividade das redes de coexpressão (LIONESS) e de interação proteína-proteína (PPI) é remodelada pela exposição ao álcool?

---

1. Setup e Preparação dos Dados

# Manipulação e Visualização
library(dplyr)
library(tidyr)
library(ggplot2)
library(tibble)
library(RColorBrewer)
library(cowplot)

# Transcriptômica e Estatística
library(DESeq2)
library(matrixStats)

# Heatmaps
library(pheatmap)
library(ComplexHeatmap)
library(circlize)

# Enriquecimento e Redes
library(clusterProfiler)
library(org.Dr.eg.db)
library(STRINGdb)
library(igraph)
library(tidygraph)
library(ggraph)

# Semente global para reprodutibilidade
set.seed(123)

2. Importação e alinhamento

# 1. Metadados
meta <- read.table("metadata_fasd.csv", sep = ",", fill = TRUE, header = TRUE, stringsAsFactors = FALSE) %>%
  transmute(
    sample = Run,
    age = Developmental_Stage,
    treatment = treatment,
    diet = Diet
  ) %>%
  dplyr::filter(age == "13 dpf")

# Limpeza dos fatores para o DESeq2
meta <- meta %>%
  mutate(
    treatment = factor(ifelse(treatment == "0% EtOH", "Control", "EtOH"), levels = c("Control", "EtOH")),
    diet = factor(ifelse(diet == "normal", "ND", "HFHC"), levels = c("ND", "HFHC"))
  )


# 2. Contagens
counts <- read.csv("GSE142311_20180702_OGWFAE_raw.data.csv", sep = ";", header = TRUE, check.names = FALSE)
if (!is.numeric(counts[[1]])) {
  rownames(counts) <- counts[[1]]
  counts <- counts[, -1]
}

if (nrow(meta) == ncol(counts)) {
  meta$sample <- colnames(counts)
  rownames(meta) <- meta$sample
} else {
  stop("Erro: O número de amostras no metadado (", nrow(meta), 
       ") não bate com as colunas do count matrix (", ncol(counts), ").")
}

# 3. Alinhamento Estrito
common_samples <- intersect(colnames(counts), meta$sample)
counts <- counts[, common_samples, drop = FALSE]

meta <- meta %>% 
  filter(sample %in% common_samples) %>% 
  arrange(match(sample, colnames(counts)))

rownames(meta) <- meta$sample
stopifnot(all(meta$sample == colnames(counts)))

3. Controle de Qualidade (QC)

lib_size <- colSums(counts)
qc_tbl <- meta %>% mutate(lib_size = lib_size, genes_detected = colSums(counts > 0))

p_lib <- ggplot(qc_tbl, aes(x = sample, y = lib_size, fill = treatment)) +
  geom_col() + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Tamanho de Biblioteca", y = "Reads", x = "", fill = "Tratamento")

p_genes <- ggplot(qc_tbl, aes(x = sample, y = genes_detected, fill = treatment)) +
  geom_col() + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Genes Detectados", y = "Genes", x = "", fill = "Tratamento")

plot_grid(p_lib, p_genes)

# Transformação VST
dds_tmp <- DESeqDataSetFromMatrix(countData = counts, colData = meta, design = ~1)
vsd <- vst(dds_tmp, blind = TRUE)
vst_mat <- assay(vsd)

# PCA
pca <- prcomp(t(vst_mat))
pca_tbl <- as.data.frame(pca$x) %>% mutate(sample = rownames(pca$x)) %>% left_join(meta, by = "sample")

p_pca <- ggplot(pca_tbl, aes(x = PC1, y = PC2, color = diet, shape = treatment)) +
  geom_point(size = 4, alpha = 0.8) +
  scale_color_manual(values = c("ND" = "#FC8D62", "HFHC" = "#66C2A5")) +
  theme_bw() +
  labs(title = "PCA (VST) - Foco no Agrupamento Global")

p_pca

4. Expressão Diferencial (DESeq2)

# Carregue este pacote junto com os outros no início do script
library(ggrepel)

# Isolando a análise para a Dieta Normal (ND)
idx_nd <- meta$diet == "ND"
dds_nd <- DESeqDataSetFromMatrix(
  countData = counts[, idx_nd, drop = FALSE],
  colData = meta[idx_nd, ],
  design = ~ treatment
)

# Filtro de baixa contagem
dds_nd <- dds_nd[rowSums(counts(dds_nd)) >= 10, ]
dds_nd <- DESeq(dds_nd)

# Resultados com Shrinkage
res_nd <- lfcShrink(dds_nd, coef = "treatment_EtOH_vs_Control", type = "apeglm")

res_nd_tbl <- as.data.frame(res_nd) %>%
  rownames_to_column("gene") %>%
  # Remove NAs para não distorcer os limites do gráfico
  filter(!is.na(padj)) %>%
  mutate(sig = case_when(
    padj < 0.05 & log2FoldChange > 0 ~ "Up-regulado",
    padj < 0.05 & log2FoldChange < 0 ~ "Down-regulado",
    TRUE ~ "NS"
  ))

# 1. Separando os 10 genes mais significativos
top_10 <- res_nd_tbl %>%
  arrange(padj) %>%
  slice_head(n = 10)

# 2. Encontrando os limites reais do eixo X para não sobrar espaço vazio
limites_x <- range(res_nd_tbl$log2FoldChange, na.rm = TRUE)

# Volcano Plot
ggplot(res_nd_tbl, aes(x = log2FoldChange, y = -log10(padj))) +
  
  # Camada base: todos os genes
  geom_point(aes(color = sig), alpha = 0.6, size = 2) +
  
  # Camada de destaque: apenas os top 10, usando shape 21 (preenchimento + borda preta)
  geom_point(data = top_10, aes(fill = sig), shape = 21, color = "black", size = 2.5, stroke = 0.8) +
  
  # Camada de texto: labels apenas para os top 10
  geom_text_repel(data = top_10, aes(label = gene), 
                  size = 3.5, 
                  fontface = "bold",
                  box.padding = 0.5, 
                  max.overlaps = Inf) +
  
  # Definindo cores
  scale_color_manual(values = c("Down-regulado" = "#2166AC", "NS" = "grey80", "Up-regulado" = "#B2182B")) +
  scale_fill_manual(values = c("Down-regulado" = "#2166AC", "NS" = "grey80", "Up-regulado" = "#B2182B"), guide = "none") +
  
  # Cortando o excesso de espaço nos eixos
  coord_cartesian(xlim = limites_x) +
  
  theme_bw() +
  labs(title = "Volcano plot: FAE vs Control (Dieta Normal)", color = "Significância")

5. Enriquecimento Funcional e Redes

# -------------------------------------------------------------------------
# 1. ORA (Over-Representation Analysis)
# -------------------------------------------------------------------------
sig_genes <- res_nd_tbl %>% filter(sig != "NS") %>% pull(gene)

entrez <- bitr(sig_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Dr.eg.db)

if(nrow(entrez) > 0) {
  ego <- enrichGO(
    gene = entrez$ENTREZID, OrgDb = org.Dr.eg.db, ont = "BP",
    pAdjustMethod = "BH", pvalueCutoff = 0.05, readable = TRUE
  )
  
  if(!is.null(ego) && nrow(as.data.frame(ego)) > 0) {
    p_ora <- dotplot(ego, showCategory = 15) + 
      labs(title = "ORA - Processos Biológicos (FAE vs Control)")
    print(p_ora)
  } else {
    cat("Nenhum termo significativo encontrado no ORA.\n")
  }
}

# -------------------------------------------------------------------------
# 2. GSEA (Gene Set Enrichment Analysis)
# -------------------------------------------------------------------------
# Ranking pelo log2FoldChange (usando todos os genes testados, sem filtro de p-valor)
rank_tbl <- res_nd_tbl %>%
  filter(!is.na(log2FoldChange)) %>%
  arrange(desc(log2FoldChange))

rank_vec <- rank_tbl$log2FoldChange
names(rank_vec) <- rank_tbl$gene

# Conversão de SYMBOL para ENTREZID para o GSEA
rank_entrez <- bitr(names(rank_vec), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Dr.eg.db)

# Filtragem e substituição de nomes no vetor
rank_vec <- rank_vec[names(rank_vec) %in% rank_entrez$SYMBOL]
names(rank_vec) <- rank_entrez$ENTREZID[match(names(rank_vec), rank_entrez$SYMBOL)]

# O vetor obrigatoriamente precisa estar ordenado do maior LFC para o menor
rank_vec <- sort(rank_vec, decreasing = TRUE)

# Execução do GSEA
gsea <- gseGO(
  geneList = rank_vec,
  OrgDb = org.Dr.eg.db,
  ont = "BP",
  pAdjustMethod = "BH",
  verbose = FALSE
)

# Separação dos top pathways (Up e Down) e plot customizado
if(!is.null(gsea) && nrow(as.data.frame(gsea)) > 0) {
  gsea_tbl <- as.data.frame(gsea)
  
  top_up <- gsea_tbl %>% filter(NES > 0) %>% arrange(desc(NES)) %>% slice_head(n = 10)
  top_down <- gsea_tbl %>% filter(NES < 0) %>% arrange(NES) %>% slice_head(n = 10)
  
  gsea_top <- bind_rows(top_up, top_down) %>%
    mutate(Description = factor(Description, levels = Description[order(NES)]))
  
  p_gsea <- ggplot(gsea_top, aes(x = NES, y = Description, color = p.adjust, size = setSize)) +
    geom_point() +
    scale_color_viridis_c(option = "C") +
    theme_bw() +
    labs(
      title = "GSEA - Top Pathways Ativados e Reprimidos", 
      x = "Normalized Enrichment Score (NES)", 
      y = NULL
    )
  print(p_gsea)
} else {
  cat("Nenhum termo significativo encontrado no GSEA.\n")
}

5.1. Rede de Interação Proteína-Proteína (PPI)

string_db <- STRINGdb$new(version = "11.5", species = 7955, score_threshold = 400, input_directory = "")

# 1. Mapeamento
mapped <- string_db$map(data.frame(gene = sig_genes), "gene", removeUnmappedRows = TRUE) %>%
  dplyr::distinct(STRING_id, .keep_all = TRUE)

Warning:  we couldn't map to STRING 4% of your identifiers

ppi <- string_db$get_interactions(mapped$STRING_id)

# 2. Preparar os nós trazendo o log2FoldChange do DESeq2 (res_nd_tbl)
nodes_tbl <- data.frame(STRING_id = unique(c(ppi$from, ppi$to))) %>%
  dplyr::left_join(mapped, by = "STRING_id") %>%
  # Trazendo os resultados de expressão (certifique-se de que res_nd_tbl está carregado)
  dplyr::left_join(res_nd_tbl %>% dplyr::select(gene, log2FoldChange), by = "gene") %>%
  dplyr::mutate(label = ifelse(is.na(gene), STRING_id, gene))

# 3. Construir o grafo base
g <- graph_from_data_frame(d = ppi, vertices = nodes_tbl, directed = FALSE)
g <- igraph::simplify(g)

# 4. Calcular métricas e filtrar rótulos com tidygraph
tg <- as_tbl_graph(g) %>%
  activate(nodes) %>%
  mutate(
    degree = centrality_degree(),
    # Define quem terá o nome plotado. Aqui, filtramos os top 15 genes com maior grau (hubs)
    is_hub = rank(-degree) <= 15,
    plot_label = ifelse(is_hub, label, NA)
  )

# 5. Plotagem Estética
set.seed(123) # Fixa o layout para não mudar a cada execução
ggraph(tg, layout = "fr") +
  # Arestas mais transparentes para não roubar a atenção
  geom_edge_link(alpha = 0.15, color = "gray60") +
  
  # Nós usando shape 21 (permite preenchimento interno com cor + borda branca)
  geom_node_point(aes(size = degree, fill = log2FoldChange), shape = 21, color = "white", stroke = 0.5) +
  
  # Escala divergente: Down = Azul, Neutro = Branco, Up = Vermelho
  scale_fill_gradient2(
    low = "#2166AC", mid = "white", high = "#B2182B", midpoint = 0,
    name = "Log2 Fold Change"
  ) +
  
  # Escala de tamanho ajustada para não engolir o plot
  scale_size_continuous(range = c(2, 9), name = "Grau (Degree)") +
  
  # Textos apenas nos Hubs, com sombra branca (bg.color) para destacar por cima das linhas
  geom_node_text(
    aes(label = plot_label), 
    repel = TRUE, 
    size = 4, 
    fontface = "bold", 
    color = "black", 
    bg.color = "white", 
    bg.r = 0.15,
    max.overlaps = Inf
  ) +
  
  theme_void() +
  labs(
    title = "Rede PPI - Genes Significativos na Dieta Normal",
    subtitle = "Top 15 Hubs mapeados pela direção da regulação (Up/Down)"
  ) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "right"
  )

5.2 Redes Específicas por Amostra (LIONESS)

# Função otimizada LIONESS (Pearson)
lioness_cor_optimized <- function(se) {
  X <- assay(se); n <- ncol(X)
  S <- rowSums(X); SS <- rowSums(X^2); SP <- X %*% t(X)
  num_full <- SP - outer(S, S)/n
  den_full <- sqrt((SS - S^2/n) %o% (SS - S^2/n))
  C_full <- num_full / den_full
  ut_idx <- which(upper.tri(C_full), arr.ind = TRUE)
  genes <- rownames(X)
  reg_names <- genes[ut_idx[, 1]]; tar_names <- genes[ut_idx[, 2]]
  v_full <- C_full[ut_idx]
  out <- matrix(NA_real_, nrow = length(reg_names), ncol = n)
  colnames(out) <- colnames(X)
  
  for (q in seq_len(n)) {
    xq <- X[, q]; S_q <- S - xq; SS_q <- SS - xq^2; SP_q <- SP - xq %*% t(xq)
    num_q <- SP_q - outer(S_q, S_q)/(n - 1)
    den_q <- sqrt((SS_q - S_q^2/(n - 1)) %o% (SS_q - S_q^2/(n - 1)))
    C_q <- num_q / den_q
    out[, q] <- n * (v_full - C_q[ut_idx]) + C_q[ut_idx]
  }
  list(W = out, reg = reg_names, tar = tar_names)
}

# Preparando os dados (Usaremos os top 1000 genes mais variáveis da Dieta ND para performance)
vst_nd <- vst_mat[, idx_nd]
gene_var <- rowVars(vst_nd)
top_idx <- order(gene_var, decreasing = TRUE)[1:1000]
expr_top <- vst_nd[top_idx, , drop = FALSE]

se_filtered <- SummarizedExperiment(
  assays = list(counts = expr_top),
  colData = DataFrame(sample = colnames(expr_top), row.names = colnames(expr_top)),
  rowData = DataFrame(gene = rownames(expr_top), row.names = rownames(expr_top))
)

# Rodando LIONESS
lion_res <- lioness_cor_optimized(se_filtered)

# Calculando o Grau (Degree) por amostra
W <- lion_res$W
all_genes <- unique(c(lion_res$reg, lion_res$tar))
deg_list <- lapply(seq_len(ncol(W)), function(i) {
  col_w <- abs(W[, i])
  thr <- mean(col_w, na.rm = TRUE) + (1.5 * sd(col_w, na.rm = TRUE)) # Threshold adaptativo
  keep <- !is.na(col_w) & col_w > thr
  active_genes <- c(lion_res$reg[keep], lion_res$tar[keep])
  gene_counts <- table(active_genes)
  sample_deg <- setNames(integer(length(all_genes)), all_genes)
  sample_deg[names(gene_counts)] <- as.integer(gene_counts)
  sample_deg
})

deg_mat <- as.data.frame(do.call(cbind, deg_list))
colnames(deg_mat) <- colnames(W)
rownames(deg_mat) <- all_genes

# Extraindo Differentially Interacting Genes (DIGs)
meta_nd <- meta[idx_nd, ]
sf_idx <- which(meta_nd$treatment == "EtOH")
gc_idx <- which(meta_nd$treatment == "Control")

dig_res <- lapply(seq_len(nrow(deg_mat)), function(i) {
  gene_name <- rownames(deg_mat)[i]
  group_sf <- deg_mat[i, sf_idx]; group_gc <- deg_mat[i, gc_idx]
  if (sd(group_sf) == 0 && sd(group_gc) == 0) return(NULL)
  
  tryCatch({
    t_test <- t.test(group_sf, group_gc)
    data.frame(gene = gene_name, p_val = t_test$p.value, diff_deg = mean(group_sf) - mean(group_gc))
  }, error = function(e) NULL)
}) %>% bind_rows() %>% mutate(padj = p.adjust(p_val, method = "BH")) %>% arrange(p_val)

# 1. Checa se algum gene passou no padj rigoroso. Se não, avisa no console.
qtd_sig <- dig_res %>% filter(padj < 0.05) %>% nrow()
if (qtd_sig == 0) {
  message("Aviso: Nenhum gene atingiu padj < 0.05. Plotando os top 6 pelo p-valor nominal (p_val).")
}

# 2. Pega os top 6 genes pelo ranking de significância bruta (mesmo que padj > 0.05)
top_digs <- dig_res %>% 
  arrange(p_val) %>% 
  slice_head(n = 6) %>% 
  pull(gene)

# Plot dos Top DIGs
if (length(top_digs) > 0) {
  plot_tbl <- as.data.frame(t(deg_mat[top_digs, , drop = FALSE])) %>%
    rownames_to_column("sample") %>%
    pivot_longer(cols = -sample, names_to = "gene", values_to = "degree") %>%
    left_join(meta_nd %>% dplyr::select(sample, treatment), by = "sample")

  p_lioness <- ggplot(plot_tbl, aes(x = treatment, y = degree, fill = treatment)) +
    geom_boxplot(outlier.size = 0.6) + 
    geom_jitter(width = 0.1, alpha = 0.6) +
    facet_wrap(~ gene, scales = "free_y") + 
    scale_fill_manual(values = c("Control" = "#4575B4", "EtOH" = "#D73027")) +
    theme_bw() +
    labs(
      title = "Top DIGs: Alteração na Conectividade da Rede (LIONESS)", 
      y = "Grau (Degree)", 
      x = ""
    )
    
  # Força a exibição do gráfico
  print(p_lioness)
}