# Module 2: CPM, TPM and size factors, by hand, and checked against the tools.
# Videos: T3 (TPM) and T5 (size factors).  Working directory: the "workshop" folder.
suppressPackageStartupMessages(library(DESeq2))

# ---- A. TPM from kallisto's own numbers ------------------------------------------------------
ab <- read.delim("data/kallisto/SRR6068402/abundance.tsv")
head(ab)
# kallisto divides by the *effective* length: the number of positions a fragment can start from
rpk <- ab$est_counts / (ab$eff_length / 1000)
my_tpm <- rpk / sum(rpk) * 1e6
cat("sum of my TPM:", sum(my_tpm), "\n")
cat("max difference from kallisto's TPM:", max(abs(my_tpm - ab$tpm)), "\n")
# EXERCISE A1: repeat with ab$length instead of ab$eff_length. How different is it, and for which transcripts?

# ---- B. why do gene-level TPM columns not sum to a million? -----------------------------------
tpm <- read.delim("data/tpm_gene.txt", row.names = 1, check.names = FALSE)
colSums(tpm)
# EXERCISE B1: what share of transcript-level TPM in SRR6068402 belongs to transcripts that are
# NOT in data/mart_export.txt?  (Hint: %in%)
mart <- read.delim("data/mart_export.txt")
cat("TPM share of non protein coding transcripts:",
    sum(ab$tpm[!ab$target_id %in% mart$Transcript.stable.ID.version]) / 1e6, "\n")

# ---- C. CPM ---------------------------------------------------------------------------------
meta <- read.delim("data/metadata.txt", row.names = 1)
count <- round(as.matrix(read.delim("data/count_gene.txt", row.names = 1, check.names = FALSE)[, rownames(meta)]))
lib <- colSums(count)
lib
cpm <- t(t(count) / lib * 1e6)
colSums(cpm)

# ---- D. size factors by the median of ratios -------------------------------------------------
my_size_factors <- function(m) {
  m <- m[rowSums(m == 0) == 0, ]                 # skip genes with a zero in any sample
  log_geo_mean <- rowMeans(log(m))               # geometric mean on the log scale
  apply(m, 2, function(x) exp(median(log(x) - log_geo_mean)))
}
sf_mine <- my_size_factors(count)
dds <- DESeqDataSetFromMatrix(count, data.frame(row.names = rownames(meta), conds = factor(meta$condition)), ~conds)
dds <- estimateSizeFactors(dds)
round(cbind(mine = sf_mine, DESeq2 = sizeFactors(dds), libsize_ratio = lib / mean(lib)), 3)

# EXERCISE D1: the worked example from the video uses 5 genes and 2 samples. Reproduce it:
ex <- count[c("Actb", "Gapdh", "Myh6", "Nppa", "Postn"), c("SRR6068402", "SRR6068404")]
ex
my_size_factors(ex)
# EXERCISE D2: which samples have the biggest gap between size factor and library-size ratio?
#              Look at the share of reads taken by Myh6 in each sample. Can you explain the gap?
round(100 * count["Myh6", ] / lib, 2)