Влияние домена L7 в белке Nos на транскриптом Drosophila

Визуализация дифференциальной экспрессии при использовании TPM

Для данной задачи будет использована только библиотека tidyverse

library(tidyverse)

Чтение метаданных

mtd <- read_delim("07_marhabaie_2024_metadata.tsv", delim = ";")
mtd <- mtd %>% filter(genotype %in% c("Upf1-NosL7", "Upf1-Nos"))

Чтение данных TPM

tpm_f <- read_tsv("07_marhabaie_2024_cycle9_genes_TPM.tsv")

Переименуем столбцы - теперь их названия совпадают с названиями строк в метаданных

tpm_f <- tpm_f %>%
  rename_at(
    vars(contains("genes.results"))
    , funs(str_replace(., ".genes.results", ""))
  )

Удаление лишнего столбца и сохранение идентификаторов строк

tpm <- tpm_f %>% select(-...1)
rownames(tpm) <- tpm_f$...1

Упорядочивание метаданных

mtd <- mtd %>% filter(sample %in% colnames(tpm))
tpm <- tpm[, mtd$sample]

Разделение на группы по генотипам

group_means <- tpm %>%
  as.data.frame() %>%
  mutate(gene = rownames(.)) %>%
  pivot_longer(-gene, names_to = "sample", values_to = "TPM") %>%
  left_join(mtd, by = c("sample" = "sample")) %>%
  group_by(gene, genotype) %>%
  summarise(mean_TPM = mean(TPM, na.rm = TRUE)) %>%
  pivot_wider(names_from = genotype, values_from = mean_TPM)

Добавление log2 среднего значения и fold change

group_means <- group_means %>%
  mutate(
    avg_log2_TPM = log2((`Upf1-Nos` + `Upf1-NosL7`) / 2 + 1),
    log2_fold_change = log2((`Upf1-Nos` + 1) / (`Upf1-NosL7` + 1))
  )

Добавим данные для графика

group_means <- group_means %>%
  mutate(
    stat_status = case_when(
      log2_fold_change > 1 ~ "Upregulated",
      log2_fold_change < -1 ~ "Downregulated",
      TRUE ~ "Not Significant"
    )
  )

Параметры визуализации

theme_set(theme_bw(base_size = 18))
tmpPointSize <- 0.3
point_color_alpha <- 1
annotationSize <- 5
hlines_col <- "blue"
smearplot_cb_palette <- c("red", "grey", "blue")
group_means %>% 
  ggplot(aes(x = avg_log2_TPM, y = log2_fold_change, color = stat_status, fill = stat_status))+
  geom_point(shape = 16, alpha = point_color_alpha, size = tmpPointSize) +
  geom_hline(yintercept = c(1, -1), col = hlines_col) +
  annotate("text", 15, 8, label = "Upregulated Genes", col = hlines_col, size = annotationSize, hjust = 1) +
  annotate("text", 15, -8, label = "Downregulated Genes", col = hlines_col, size = annotationSize, hjust = 1) +
  labs(
    title = "Upf1-NosL7 vs Upf1-Nos",
    x = expression("Average log"[2] * "(Transcript per Million)"),
    y = expression("log"[2] * "(Fold Change)")
  ) +
  scale_colour_manual(values = smearplot_cb_palette) +
  scale_y_continuous(limits = c(-14, 14), minor_breaks = c(0.0000001)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_line(colour = "black", size = 0.5),
    legend.position = "none",
    legend.title = element_blank()
  ) +
  guides(color = guide_legend(override.aes = list(size = 4)))

Проверим, что экспрессия падает - поставим тест Манна-Уитни.

H0: Медиана log2_fold_change = 0

H1: Медиана log2_fold_change < 0

test_result <- wilcox.test(
  group_means$log2_fold_change,
  mu = 0,
  alternative = "less"
)
test_result
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  group_means$log2_fold_change
## V = 22170164, p-value < 2.2e-16
## alternative hypothesis: true location is less than 0

Анализ дифференциальной экспрессии по каунтам и функциональная аннотация

С помощью DESeq2 сравним экспрессию у Epf1-NosL7 и Epf1-Nos.

Читаем и фильтруем метаданные

mtd <- read_delim("07_marhabaie_2024_metadata.tsv", delim = ";")
mtd <- mtd %>% filter(developmental_stage == "cycle 9")

Чтение и преобразование данных о каунтах

cnt_9_f <- read_tsv("07_marhabaie_2024_cycle9_genes_counts.tsv")
cnt9 <- cnt_9_f %>% select(-...1)
cnt9 <- round(cnt9)
cnt9 <- cnt9 %>%
  rename_at(
    vars(contains("genes.results"))
    , funs(str_replace(., ".genes.results", ""))
  )
cnt9 <- as.matrix(cnt9)
rownames(cnt9) <- cnt_9_f$...1

Сортируем строки и столбцы в двух таблицах шли в одном порядке

mtd <- mtd[match(colnames(cnt9), mtd$sample), ]

Проводим анализ диф экспрессии

library(DESeq2)
dds <- DESeqDataSetFromMatrix(
  countData = cnt9,
  colData = mtd,
  design = ~genotype
)
dds <- DESeq(dds)
res <- results(dds, contrast = c("genotype", "Upf1-NosL7", "Upf1-Nos"))

Добавление log10(p-value) и преобразование в data.frame

res_df <- as.data.frame(res) %>%
  mutate(
    gene = rownames(cnt9),
    log10_pvalue = -log10(pvalue),
    significant = ifelse(!is.na(padj) & padj < 0.05 & abs(log2FoldChange) > 1, "Significant", "Not Significant")
  )

Построение Volcano-Plot

volcano_plot <- ggplot(res_df, aes(x = log2FoldChange, y = log10_pvalue, color = significant)) +
  geom_point(alpha = 0.7, size = 2.0) +  # Прозрачность и размер точек
  scale_color_manual(values = c("red", "blue")) +
  theme_bw(base_size = 14) +  # Базовый стиль и размер шрифта
  theme(
    panel.grid.major = element_line(color = "grey80", linetype = "dotted"),
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold", hjust = 0.5, size = 12), # Уменьшен размер заголовка
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14)
  ) +
  labs(
    title = "Volcano Plot: Upf1-NosL7 vs Upf1-Nos",
    x = expression(Log[2] ~ "Fold Change"),
    y = expression(-Log[10] ~ "(P-value)"),
    color = "Significance"
  ) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black", size = 1) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black", size = 1) +
  xlim(-10, 10)

volcano_plot

Фильтрация значимых генов

significant_genes <- res_df %>%
  filter(abs(log2FoldChange) > 1 & !is.na(padj) & padj < 0.05) %>%
  pull(gene)
library(AnnotationDbi)
library(biomaRt) 
library(clusterProfiler)
library(org.Dm.eg.db)

Преобразование Ensembl в Entrez ID

gene_entrez <- bitr(significant_genes, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Dm.eg.db)

Обогащение по Gene Ontology (биологические процессы)

ego <- enrichGO(
  gene          = gene_entrez$ENTREZID,
  OrgDb         = org.Dm.eg.db,
  ont           = "BP",  # "BP" - биологические процессы, можно также использовать "MF" и "CC"
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.2,
  readable      = TRUE   # Преобразование Entrez ID обратно в символы генов
)
barplot(ego, showCategory = 10, title = "Top 10 GO Categories\nfor enrichment between\nUpf1-NosL7 and Upf1-Nos")

dotplot(ego, showCategory = 10, title = "Top 10 GO Categories\nfor enrichment between\nUpf1-NosL7 and Upf1-Nos")

Сделаем то же самое для белков Upf1 и Upf1-Nos - найдем и у них обогащенные категории. Сравнивать лучше с Upf1, а не wt, поскольку белок Upf1 сам может влиять на транскриптом.

## Warning: Removed 6101 rows containing missing values or values outside the scale range
## (`geom_point()`).

### Выводы