# Те же данные, но через DESeq2 — как считают в настоящих работах.
#
# DESeq2 живёт в Bioconductor и написан на C++; для браузера его никто не
# собирал, поэтому в практикуме он выполняется здесь, на обычном компьютере,
# а на страницу уходит его собственный вывод. Скрипт короткий нарочно: это
# ровно то, что школьник может запустить у себя в RStudio.
#
#   Rscript analysis/deseq.R      (из корня проекта)
#
# Пишет в site/assets/r/:
#   deseq_muha.csv     таблица результатов + оценки разброса по каждому гену
#   deseq_mysh.csv     то же для двух мышиных опытов
#   deseq_faktory.csv  размерные факторы: чем DESeq2 делит вместо CPM
#   deseq_norm.csv     счётчики после нормировки DESeq2 — из них строят plotCounts
#   deseq_vst.csv      те же данные после vst — из них строят PCA и тепловые карты
#
# Последние два нужны, чтобы на странице можно было рисовать свои графики по
# данным DESeq2, а не только смотреть на готовую таблицу результатов.

suppressPackageStartupMessages(library(DESeq2))

HERE <- "site/assets/r"
stopifnot(dir.exists(HERE))

## ---------- муха ----------
counts <- as.matrix(read.csv(file.path(HERE, "muhi.csv"), row.names = 1))

# Таблица образцов: одна строка на столбец счётчиков. Порядок обязан совпадать.
info <- data.frame(row.names = colnames(counts),
                   group = factor(c(rep("norm", 3), rep("mut", 3)), levels = c("norm", "mut")))

# design = ~ group означает «объясняй различия между генами группой».
# Слева от тильды ничего не пишут: что объясняем — это сами счётчики.
dds <- DESeqDataSetFromMatrix(counts, info, design = ~ group)
dds <- dds[rowSums(counts(dds)) >= 10, ]          # тот же фильтр, что в практикуме
dds <- DESeq(dds, quiet = TRUE)                   # три шага сразу, см. ниже

res <- results(dds, contrast = c("group", "mut", "norm"))
disp <- mcols(dds)                                # оценки разброса по каждому гену

out <- data.frame(
  gene = rownames(res),
  baseMean = round(res$baseMean, 2),
  log2FoldChange = round(res$log2FoldChange, 4),
  lfcSE = round(res$lfcSE, 4),
  stat = round(res$stat, 3),
  pvalue = signif(res$pvalue, 5),
  padj = signif(res$padj, 5),
  dispGeneEst = signif(disp$dispGeneEst, 5),
  dispFit = signif(disp$dispFit, 5),
  dispersion = signif(disp$dispersion, 5)
)
write.csv(out, file.path(HERE, "deseq_muha.csv"), row.names = FALSE, na = "NA")

# Сжатые изменения: для слабых генов DESeq2 умеет притягивать log2FoldChange к нулю,
# чтобы шумные гены не выглядели чемпионами. Тип "normal" встроен, пакетов не надо.
shr <- lfcShrink(dds, contrast = c("group", "mut", "norm"), type = "normal", quiet = TRUE)
out$lfcShrink <- round(shr$log2FoldChange, 4)
write.csv(out, file.path(HERE, "deseq_muha.csv"), row.names = FALSE, na = "NA")

# Нормированные счётчики: то же самое, что counts / sizeFactor.
norm <- round(counts(dds, normalized = TRUE), 1)
write.csv(data.frame(gene = rownames(norm), norm), file.path(HERE, "deseq_norm.csv"),
          row.names = FALSE)

# vst — «стабилизация дисперсии»: логарифм, который не раздувает шум у слабых генов.
# Именно на ней DESeq2 строит PCA и тепловые карты, а не на сырых счётчиках.
v <- round(assay(vst(dds, blind = FALSE)), 3)
write.csv(data.frame(gene = rownames(v), v), file.path(HERE, "deseq_vst.csv"), row.names = FALSE)

sig <- subset(out, !is.na(padj) & padj < 0.05 & abs(log2FoldChange) > 1)
cat("муха: генов после фильтра", nrow(out),
    "| padj<0.05 и вдвое:", nrow(sig),
    "| вверх", sum(sig$log2FoldChange > 0), "\n")

faktory <- data.frame(sample = colnames(counts),
                      opyt = "муха",
                      depth = colSums(counts),
                      sizeFactor = round(sizeFactors(dds), 4))

## ---------- мышь ----------
mouse <- as.matrix(read.csv(file.path(HERE, "mysh.csv"), row.names = 1))

run <- function(m, tag) {
  info <- data.frame(row.names = colnames(m),
                     group = factor(c(rep("wt", 3), rep("mut", 3)), levels = c("wt", "mut")))
  d <- DESeqDataSetFromMatrix(m, info, design = ~ group)
  d <- d[rowSums(counts(d)) >= 10, ]
  d <- DESeq(d, quiet = TRUE)
  r <- results(d, contrast = c("group", "mut", "wt"))
  # pvalue не сохраняем: на страницу уходит по сети, а нужен только padj
  res <- data.frame(gene = rownames(r),
                    baseMean = round(r$baseMean, 1),
                    log2FoldChange = round(r$log2FoldChange, 3),
                    padj = signif(r$padj, 4))
  # считаем по округлённым числам — ровно по тем, что увидит страница
  n <- sum(!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1)
  cat("мышь", tag, ": генов", nrow(r), "| padj<0.05 и вдвое:", n, "\n")
  list(res = res,
       sf = data.frame(sample = colnames(m), opyt = tag,
                       depth = colSums(m), sizeFactor = round(sizeFactors(d), 4)))
}

on <- run(mouse[, 1:6], "сенсор цел")
off <- run(mouse[, 7:12], "сенсора нет")

# Оба опыта в одном файле: колонка opyt различает их.
mysh <- rbind(cbind(opyt = "on", on$res), cbind(opyt = "off", off$res))
write.csv(mysh, file.path(HERE, "deseq_mysh.csv"), row.names = FALSE, na = "NA")

write.csv(rbind(faktory, on$sf, off$sf), file.path(HERE, "deseq_faktory.csv"), row.names = FALSE)

for (f in c("deseq_muha.csv", "deseq_mysh.csv", "deseq_faktory.csv",
            "deseq_norm.csv", "deseq_vst.csv")) {
  cat(sprintf("  %-18s %6.0f КБ\n", f, file.size(file.path(HERE, f)) / 1024))
}
cat("DESeq2 версии", as.character(packageVersion("DESeq2")), "\n")