# Мухи без Adar против обычных мух. 3 повтора на группу.
# Данные: GSE128332 (головы, Adar5G1 null vs w1118), сырые счётчики featureCounts.
# Только базовый R. Запуск: Rscript urok.R  (из этого каталога)

## ШАГ 1. Что вообще посчитано?
d <- read.csv("muhi.csv")
d[d$gene %in% c("Adar", "DptA", "Act5C"), ]

## ШАГ 2. Можно ли сравнивать столбцы напрямую?
glub <- colSums(d[, 2:7])                       # сколько всего прочтений в каждом образце
cpm  <- t((t(d[, 2:7]) + 1) / glub * 1e6)       # доля гена на миллион прочтений
rownames(cpm) <- d$gene
rbind(syrye = as.numeric(d[d$gene == "ninaE", 2:7]), CPM = round(cpm["ninaE", ]))

## ШАГ 3. Про какие гены вообще можно что-то сказать?
porog <- 1                                      # ползунок: 0.5 / 1 / 2 / 5
est <- rowSums(cpm >= porog) >= 3
cpm <- cpm[est, ]
c(ostalos = nrow(cpm), vybrosili = sum(!est))

## ШАГ 4. Во сколько раз изменился ген?
lg <- log2(cpm)
fc <- rowMeans(lg[, 4:6]) - rowMeans(lg[, 1:3])
round(head(sort(fc, decreasing = TRUE), 10), 2)
c(log2 = fc["DptA"], v_razah = 2^fc["DptA"])

## ШАГ 5. А вдруг это случайность?
t.test(lg["DptA", 4:6], lg["DptA", 1:3], var.equal = TRUE)
p <- apply(lg, 1, function(x) t.test(x[4:6], x[1:3], var.equal = TRUE)$p.value)
c(znachimyh = sum(p < 0.05), zhdem_sluchayno = round(0.05 * length(p)))

## ШАГ 6. Мы проверили 10 000 генов — сколько «значимых» получилось само собой?
q <- p.adjust(p, method = "BH")
sapply(c(0.05, 0.1, 0.2), function(Q) sum(q < Q))
geny <- names(q)[q < 0.1 & abs(fc) > 1]
c(vsego = length(geny), vverh = sum(fc[geny] > 0), vniz = sum(fc[geny] < 0))

## ШАГ 7. Что это за гены, если смотреть на них вместе?
go  <- read.csv("go.csv"); go <- go[go$gene %in% rownames(lg), ]
vse <- unique(go$gene); n <- table(go$process)
k   <- table(go$process[go$gene %in% geny]); k <- k[k > 0]
pgo <- phyper(k - 1, n[names(k)], length(vse) - n[names(k)], sum(vse %in% geny), lower.tail = FALSE)
head(data.frame(process = names(k), nashli = as.integer(k), vsego = as.integer(n[names(k)]),
                p = signif(as.numeric(pgo), 2))[order(pgo), ], 10)

## ШАГ 8. Как эти гены связаны друг с другом?
sub <- go[go$gene %in% geny, ]
M <- table(sub$gene, sub$process)               # ген x процесс
A <- M %*% t(M); diag(A) <- 0                   # сколько общих процессов у пары генов
ord <- hclust(dist(M))$order; A <- A[ord, ord]  # похожие гены — рядом на круге
ug <- rownames(A); K <- length(ug); u <- 2 * pi * seq_len(K) / K
x <- cos(u); y <- sin(u)
sviaz <- 5                                      # ползунок: сколько общих процессов считаем связью
e <- which(A >= sviaz & upper.tri(A), arr.ind = TRUE)
png("karta.png", 1100, 1100, res = 120); par(mar = c(0, 0, 0, 0))
plot(x, y, type = "n", axes = FALSE, xlab = "", ylab = "", asp = 1, xlim = c(-1.35, 1.35))
segments(x[e[, 1]], y[e[, 1]], x[e[, 2]], y[e[, 2]], col = "grey70")
points(x, y, pch = 19, cex = 1.3, col = ifelse(fc[ug] > 0, "firebrick", "steelblue"))
text(x * 1.14, y * 1.14, ug, cex = 0.62)
dev.off()
c(uzlov = K, reber = nrow(e))