# Module F: the files of an RNA-seq experiment, opened by hand.  Videos: F1 to F8.
# Working directory: the "workshop" folder. Base R only, nothing to install.
# Data (real): the first 10,000 read pairs of sample SRR6068402, and Pkm from Ensembl release 99 (GRCm38).
d <- "data/foundations"

# ---- A. FASTQ: four lines per read ---------------------------------------------------------------
fq1 <- readLines(file.path(d, "SRR6068402_1.first10k.fastq"))
fq2 <- readLines(file.path(d, "SRR6068402_2.first10k.fastq"))
length(fq1) / 4                       # number of reads
head(fq1, 4)                          # the first record
name1 <- fq1[seq(1, length(fq1), 4)]; seq1 <- fq1[seq(2, length(fq1), 4)]; qual1 <- fq1[seq(4, length(fq1), 4)]
name2 <- fq2[seq(1, length(fq2), 4)]; seq2 <- fq2[seq(2, length(fq2), 4)]; qual2 <- fq2[seq(4, length(fq2), 4)]
table(nchar(seq1))                    # read lengths
# are the two files in the same order? compare the read names without the /1 and /2
all(sub(" .*", "", name1) == sub(" .*", "", name2))

# ---- B. quality characters -> Phred scores -> error probabilities ----------------------------------
phred <- function(q) utf8ToInt(q) - 33
phred(qual1[1])[1:10]
10^(-phred(qual1[1])[1:10] / 10)      # error probability of the first 10 letters
table(strsplit(paste(qual1, collapse = ""), "")[[1]])   # which quality characters occur?
Q1 <- t(sapply(qual1, phred, USE.NAMES = FALSE))         # reads x positions
Q2 <- t(sapply(qual2, phred, USE.NAMES = FALSE))
plot(colMeans(Q1), type = "l", col = "#3f7cac", lwd = 2, ylim = c(0, 41), xlab = "position in read", ylab = "mean Q",
     main = "mean quality per position (first 10,000 pairs)")
lines(colMeans(Q2), col = "#c2185b", lwd = 2)
abline(h = c(20, 28, 30), lty = c(3, 3, 2), col = "grey50")
legend("bottomright", c("read 1", "read 2"), col = c("#3f7cac", "#c2185b"), lwd = 2)
c(read1_Q30 = mean(Q1 >= 30), read2_Q30 = mean(Q2 >= 30))
mean(grepl("N", seq1))                # share of reads with at least one N

# ---- C. paired-end: do the mates overlap? -----------------------------------------------------------
revcomp <- function(s) chartr("ACGTN", "TGCAN", paste(rev(strsplit(s, "")[[1]]), collapse = ""))
overlap <- function(a, b, min_ov = 15, max_mm = 2) {
  b_rc <- revcomp(b)
  for (k in nchar(a):min_ov) {
    x <- strsplit(substr(a, nchar(a) - k + 1, nchar(a)), "")[[1]]
    y <- strsplit(substr(b_rc, 1, k), "")[[1]]
    if (sum(x != y) <= max_mm) return(k)
  }
  NA
}
ov <- overlap(seq1[1], seq2[1])
c(overlap = ov, fragment_length = 2 * nchar(seq1[1]) - ov)
# EXERCISE C1: run overlap() on the first 500 pairs. What share of pairs overlap? What does that tell you
#              about the fragment lengths in this library?

# ---- D. FASTA: one entry per transcript -------------------------------------------------------------
fa <- readLines(file.path(d, "Pkm_transcripts.fa"))
hdr <- grep("^>", fa)
ids <- sub("^>(\\S+).*", "\\1", fa[hdr])
seqs <- sapply(seq_along(hdr), function(i) {
  end <- if (i < length(hdr)) hdr[i + 1] - 1 else length(fa)
  paste(fa[(hdr[i] + 1):end], collapse = "")
})
names(seqs) <- ids
sort(nchar(seqs), decreasing = TRUE)
fa[hdr[1]]                            # read the header: ID, type, location, gene, biotype, symbol
gc <- function(s) mean(strsplit(s, "")[[1]] %in% c("G", "C"))
round(sapply(seqs, gc), 3)

# ---- E. GTF: where the exons are ------------------------------------------------------------------------
gtf <- read.delim(file.path(d, "Pkm.gtf"), header = FALSE, comment.char = "#", quote = "",
                  col.names = c("chr", "source", "feature", "start", "end", "score", "strand", "frame", "attributes"))
table(gtf$feature)
attr_of <- function(x, key) sub(paste0('.*', key, ' "([^"]+)".*'), "\\1", x)
ex <- subset(gtf, feature == "exon")
ex$transcript <- attr_of(ex$attributes, "transcript_id")
ex$biotype <- attr_of(ex$attributes, "transcript_biotype")
sort(table(ex$transcript), decreasing = TRUE)          # exons per transcript
unique(ex[, c("transcript", "biotype")])
# draw the transcripts
tx <- unique(ex$transcript)
plot(NA, xlim = range(c(ex$start, ex$end)), ylim = c(0.5, length(tx) + 0.5), yaxt = "n", xlab = "chromosome 9 position (GRCm38)",
     ylab = "", main = "Pkm transcripts (Ensembl 99)")
for (i in seq_along(tx)) {
  e <- ex[ex$transcript == tx[i], ]
  segments(min(e$start), i, max(e$end), i, col = "grey70")
  rect(e$start, i - 0.3, e$end, i + 0.3, col = "#3f7cac", border = NA)
}
axis(2, at = seq_along(tx), labels = substr(tx, 13, 18), las = 1, cex.axis = 0.7)

# ---- F. SAM fields by hand -------------------------------------------------------------------------------
decode_flag <- function(f) {
  bits <- c(paired = 1, proper_pair = 2, unmapped = 4, mate_unmapped = 8, reverse = 16, mate_reverse = 32,
            read1 = 64, read2 = 128, secondary = 256, qc_fail = 512, duplicate = 1024, supplementary = 2048)
  names(bits)[bitwAnd(f, bits) > 0]
}
decode_flag(99)
decode_flag(147)
decode_flag(83)            # the flag on the RNA-seq read in the 2025 NGS lecture slides
cigar_ops <- function(cg) {
  n <- as.integer(regmatches(cg, gregexpr("[0-9]+", cg))[[1]])
  op <- regmatches(cg, gregexpr("[MIDNSHP=X]", cg))[[1]]
  c(read_length = sum(n[op %in% c("M", "I", "S", "=", "X")]),
    reference_span = sum(n[op %in% c("M", "D", "N", "=", "X")]))
}
cigar_ops("150M")
cigar_ops("75M1I74M")
cigar_ops("82M406N19M")    # the Myh6 read from video F7
# EXERCISE F1: a read has flag 147 and CIGAR 75M1I74M. In words: which mate is it, which strand, and what
#              does its alignment look like? How many reference bases does it cover?