sir3-8-M.ECOGII 25C vs 37C aggregate methylation with nanopore sequencing

# LOESS PLOT FUNCTION
# d = data, ch = chromosome, b = beginning, e = end
plot_loess <- function(d0, d1, ch, b, e){
  region_methyl_0 <- d0[chrom == ch & start > b & start < e]
  lo_methyl_0 <- loess(percentage ~ start, data = region_methyl_0, weights = region_methyl_0$coverage, enp.target = 100)
  
  region_methyl_1 <- d1[chrom == ch & start > b & start < e]
  lo_methyl_1 <- loess(percentage ~ start, data = region_methyl_1, weights = region_methyl_1$coverage, enp.target = 100)
  
  plot(x=lo_methyl_0$x[order(lo_methyl_0$x)], y=lo_methyl_0$fitted[order(lo_methyl_0$x)], 
     type = 'l', lwd=2, ylim = c(0,100), col = "darkturquoise", frame.plot = FALSE,
     xlab = "", ylab = "")
  lines(x=lo_methyl_1$x[order(lo_methyl_1$x)], y=lo_methyl_1$fitted[order(lo_methyl_1$x)],
      type = 'l', lwd=2, ylim = c(0,100), col = "mediumpurple4")
  box(bty="l")
}

# DATA
columns <- c("chrom", "start", "end", "name", "score", 
             "strand", "startCodon", "stopCodon", "color", 
             "coverage", "percentage")
select_cols <- c("chrom", "start", "coverage", "percentage")

# sir3-8-ecoGII 37C
methyl_37 <- fread("/home/mbrothers/nanopore/210205_Ocular/modified_bases.aggregate11.6mA.bed",
                col.names = columns)
methyl_filtered_37 <- methyl_37[coverage > 10, ..select_cols]

# sir3-8-ecoGII 25C
methyl_25 <- fread("/home/mbrothers/nanopore/210205_Ocular/modified_bases.aggregate08.6mA.bed",
                col.names = columns)
methyl_filtered_25 <- methyl_25[coverage > 10, ..select_cols]

# PLOT CONTROL REGION
plot_loess(methyl_filtered_37, methyl_filtered_25, "III", 80e3, 105e3)

# with genes:
# abline(v = c(79162, 82275, 83620, 83997, 91324, 92418, 92777, 94270, 94621, 95763, 96281, 101191, 
#              101317, 101788, 102075, 103358, 103571, 104350))

# PLOT HML
plot_loess(methyl_filtered_37, methyl_filtered_25, "III", 0, 25e3)

# with genes:
# abline(v = c(11146, 14849, 9706, 11082, 6479, 8326, 15798, 16880, 17290, 18561, 18816, 22106, 22429, 23379,
#              23584, 23925, 23523, 23981, 24032, 24325))
# # with telomere:
# abline(v = 1098)

# PLOT HMR
plot_loess(methyl_filtered_37, methyl_filtered_25, "III", 280e3, 305e3)

# with genes:
# abline(v = c(292674, 292769, 294805, 294864, 288170, 289258, 297049, 298605, 286762, 287937, 280117, 286443,
#              301271, 302221, 304361, 305467))
# with Ty elements:
# abline(v = c(291922, 292167, 291373, 291712, 295958, 296190, 295003, 295330), col = "red")
# with tRNA:
# abline(v = c(295484, 295556), col = "green")

# PLOT TELOMERES
chromosomes <- c("I", "II", "III", "IV", "V", "VI", "VII", "VIII",
                 "IX", "X", "XI", "XII", "XIII", "XIV", "XV", "XVI")
right_ends <- c(230218, 813184, 316620, 1531933, 576874, 270161, 1090940, 562643,
                439888, 745751, 666816, 1078177, 917000, 765000, 1091291, 948010)
for(i in 1:16){
  plot_loess(methyl_filtered_37, methyl_filtered_25, chromosomes[i], 0, 10e3)
}

for(i in 1:16){
  plot_loess(methyl_filtered_37, methyl_filtered_25, chromosomes[i], right_ends[i]-10e3, right_ends[i])
}

sir3-8-M.ECOGII time course aggregate methylation with nanopore sequencing (biological replicate 1)

# PLOT FUNCTION
# t = time point data, ch = chromosome, b = beginning, e = end
plot_timecourse <- function(t0, t1, t2, t3, t4, t5, tss, ch, b, e){
  region_methyl_0 <- t0[chrom == ch & start > b & start < e]
  lo_methyl_0 <- loess(percentage ~ start, data=region_methyl_0, weights=region_methyl_0$coverage, enp.target = 100)

  region_methyl_1 <- t1[chrom == ch & start > b & start < e]
  lo_methyl_1 <- loess(percentage ~ start, data=region_methyl_1, weights=region_methyl_1$coverage, enp.target = 100)

  region_methyl_2 <- t2[chrom == ch & start > b & start < e]
  lo_methyl_2 <- loess(percentage ~ start, data=region_methyl_2, weights=region_methyl_2$coverage, enp.target = 100)

  region_methyl_3 <- t3[chrom == ch & start > b & start < e]
  lo_methyl_3 <- loess(percentage ~ start, data=region_methyl_3, weights=region_methyl_3$coverage, enp.target = 100)

  region_methyl_4 <- t4[chrom == ch & start > b & start < e]
  lo_methyl_4 <- loess(percentage ~ start, data=region_methyl_4, weights=region_methyl_4$coverage, enp.target = 100)

  region_methyl_5 <- t5[chrom == ch & start > b & start < e]
  lo_methyl_5 <- loess(percentage ~ start, data=region_methyl_5, weights=region_methyl_5$coverage, enp.target = 100)
  
  region_methyl_ss <- tss[chrom == ch & start > b & start < e]
  lo_methyl_ss <- loess(percentage ~ start, data=region_methyl_ss, weights=region_methyl_ss$coverage, enp.target = 100)

  plot(x=lo_methyl_0$x[order(lo_methyl_0$x)], y=lo_methyl_0$fitted[order(lo_methyl_0$x)],
     type = 'l', lwd=2, ylim = c(0,100), col = "black", frame.plot = FALSE,
     xlab = "", ylab = "")
  lines(x=lo_methyl_1$x[order(lo_methyl_1$x)], y=lo_methyl_1$fitted[order(lo_methyl_1$x)],
        type = 'l', lwd=2, ylim = c(0,100), col = "forestgreen")
  # lines(x=lo_methyl_2$x[order(lo_methyl_2$x)], y=lo_methyl_2$fitted[order(lo_methyl_2$x)],
  #       type = 'l', lwd=2, ylim = c(0,100), col = "black")
  lines(x=lo_methyl_3$x[order(lo_methyl_3$x)], y=lo_methyl_3$fitted[order(lo_methyl_3$x)],
        type = 'l', lwd=2, ylim = c(0,100), col = "mediumpurple3")
  # lines(x=lo_methyl_4$x[order(lo_methyl_4$x)], y=lo_methyl_4$fitted[order(lo_methyl_4$x)],
  #       type = 'l', lwd=2, ylim = c(0,100), col = "darkturquoise")
  lines(x=lo_methyl_5$x[order(lo_methyl_5$x)], y=lo_methyl_5$fitted[order(lo_methyl_5$x)],
        type = 'l', lwd=2, ylim = c(0,100), col = "deeppink")
  lines(x=lo_methyl_ss$x[order(lo_methyl_ss$x)], y=lo_methyl_ss$fitted[order(lo_methyl_ss$x)],
        type = 'l', lwd=2, lty = 3, ylim = c(0,100), col = "gray50")
  box(bty="l")
}

# DATA
columns <- c("chrom", "start", "end", "name", "score", 
             "strand", "startCodon", "stopCodon", "color", 
             "coverage", "percentage")
select_cols <- c("chrom", "start", "coverage", "percentage")
my_pal <- c("gray50", "forestgreen", "darkturquoise", "mediumpurple3", "deeppink", "black")

# change for different samples
methyl_0 <- fread("/home/mbrothers/nanopore/210304_Amanita/modified_bases.aggregate07.6mA.bed",
                col.names = columns)
methyl_filtered_0 <- methyl_0[coverage > 10, ..select_cols]

methyl_15 <- fread("/home/mbrothers/nanopore/210304_Amanita/modified_bases.aggregate08.6mA.bed",
                  col.names = columns)
methyl_filtered_15 <- methyl_15[coverage > 10, ..select_cols]

methyl_30 <- fread("/home/mbrothers/nanopore/210304_Amanita/modified_bases.aggregate09.6mA.bed",
                  col.names = columns)
methyl_filtered_30 <- methyl_30[coverage > 10, ..select_cols]

methyl_45 <- fread("/home/mbrothers/nanopore/210304_Amanita/modified_bases.aggregate10.6mA.bed",
                  col.names = columns)
methyl_filtered_45 <- methyl_45[coverage > 10, ..select_cols]

methyl_60 <- fread("/home/mbrothers/nanopore/210304_Amanita/modified_bases.aggregate11.6mA.bed",
                   col.names = columns)
methyl_filtered_60 <- methyl_60[coverage > 10, ..select_cols]

methyl_90 <- fread("/home/mbrothers/nanopore/210304_Amanita/modified_bases.aggregate12.6mA.bed",
                    col.names = columns)
methyl_filtered_90 <- methyl_90[coverage > 10, ..select_cols]

methyl_ss <- fread("/home/mbrothers/nanopore/210205_Ocular/modified_bases.aggregate08.6mA.bed",
                    col.names = columns)
methyl_filtered_ss <- methyl_ss[coverage > 10, ..select_cols]

# PLOT CONTROL
plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, "III", 90e3, 95e3)

# PLOT HML
plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, "III", 10e3, 15.5e3)
# with genes:
abline(v = c(11146, 14849, 11082, 13282, 13809, 12386, 13018))

# PLOT HMR
plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, "III", 292e3, 295.5e3)
# with genes:
abline(v = c(292674, 292769, 294805, 294864, 293179, 293538, 293835, 294321))

# PLOT TELOMERES
chromosomes <- c("I", "II", "III", "IV", "V", "VI", "VII", "VIII",
                 "IX", "X", "XI", "XII", "XIII", "XIV", "XV", "XVI")
right_ends <- c(230218, 813184, 316620, 1531933, 576874, 270161, 1090940, 562643,
                439888, 745751, 666816, 1078177, 917000, 765000, 1091291, 948010)
for(i in 1:16){
  plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, chromosomes[i], 0, 10e3)
}

for(i in 1:16){
  plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, chromosomes[i], right_ends[i]-10e3, right_ends[i])
}

## sir3-8-M.ECOGII time course aggregate methylation with nanopore sequencing (biological replicate 2)

# PLOT FUNCTION
# t = time point data, ch = chromosome, b = beginning, e = end
plot_timecourse <- function(t0, t1, t2, t3, t4, t5, tss, ch, b, e){
  region_methyl_0 <- t0[chrom == ch & start > b & start < e]
  lo_methyl_0 <- loess(percentage ~ start, data=region_methyl_0, weights=region_methyl_0$coverage, enp.target = 100)

  region_methyl_1 <- t1[chrom == ch & start > b & start < e]
  lo_methyl_1 <- loess(percentage ~ start, data=region_methyl_1, weights=region_methyl_1$coverage, enp.target = 100)

  region_methyl_2 <- t2[chrom == ch & start > b & start < e]
  lo_methyl_2 <- loess(percentage ~ start, data=region_methyl_2, weights=region_methyl_2$coverage, enp.target = 100)

  region_methyl_3 <- t3[chrom == ch & start > b & start < e]
  lo_methyl_3 <- loess(percentage ~ start, data=region_methyl_3, weights=region_methyl_3$coverage, enp.target = 100)

  region_methyl_4 <- t4[chrom == ch & start > b & start < e]
  lo_methyl_4 <- loess(percentage ~ start, data=region_methyl_4, weights=region_methyl_4$coverage, enp.target = 100)

  region_methyl_5 <- t5[chrom == ch & start > b & start < e]
  lo_methyl_5 <- loess(percentage ~ start, data=region_methyl_5, weights=region_methyl_5$coverage, enp.target = 100)
  
  region_methyl_ss <- tss[chrom == ch & start > b & start < e]
  lo_methyl_ss <- loess(percentage ~ start, data=region_methyl_ss, weights=region_methyl_ss$coverage, enp.target = 100)

  plot(x=lo_methyl_0$x[order(lo_methyl_0$x)], y=lo_methyl_0$fitted[order(lo_methyl_0$x)],
     type = 'l', lwd=2, ylim = c(0,100), col = "black", frame.plot = FALSE,
     xlab = "", ylab = "")
  lines(x=lo_methyl_1$x[order(lo_methyl_1$x)], y=lo_methyl_1$fitted[order(lo_methyl_1$x)],
        type = 'l', lwd=2, ylim = c(0,100), col = "forestgreen")
  # lines(x=lo_methyl_2$x[order(lo_methyl_2$x)], y=lo_methyl_2$fitted[order(lo_methyl_2$x)],
  #       type = 'l', lwd=2, ylim = c(0,100), col = "black")
  lines(x=lo_methyl_3$x[order(lo_methyl_3$x)], y=lo_methyl_3$fitted[order(lo_methyl_3$x)],
        type = 'l', lwd=2, ylim = c(0,100), col = "mediumpurple3")
  # lines(x=lo_methyl_4$x[order(lo_methyl_4$x)], y=lo_methyl_4$fitted[order(lo_methyl_4$x)],
  #       type = 'l', lwd=2, ylim = c(0,100), col = "darkturquoise")
  lines(x=lo_methyl_5$x[order(lo_methyl_5$x)], y=lo_methyl_5$fitted[order(lo_methyl_5$x)],
        type = 'l', lwd=2, ylim = c(0,100), col = "deeppink")
  lines(x=lo_methyl_ss$x[order(lo_methyl_ss$x)], y=lo_methyl_ss$fitted[order(lo_methyl_ss$x)],
        type = 'l', lwd=2, lty = 3, ylim = c(0,100), col = "gray50")
  box(bty="l")
}

# DATA
columns <- c("chrom", "start", "end", "name", "score", 
             "strand", "startCodon", "stopCodon", "color", 
             "coverage", "percentage")
select_cols <- c("chrom", "start", "coverage", "percentage")
my_pal <- c("gray50", "forestgreen", "darkturquoise", "mediumpurple3", "deeppink", "black")

# change for different samples
methyl_0 <- fread("/home/mbrothers/nanopore/210310_Russula/modified_bases.aggregate01.6mA.bed",
                col.names = columns)
methyl_filtered_0 <- methyl_0[coverage > 10, ..select_cols]

methyl_15 <- fread("/home/mbrothers/nanopore/210310_Russula/modified_bases.aggregate02.6mA.bed",
                  col.names = columns)
methyl_filtered_15 <- methyl_15[coverage > 10, ..select_cols]

methyl_30 <- fread("/home/mbrothers/nanopore/210310_Russula/modified_bases.aggregate03.6mA.bed",
                  col.names = columns)
methyl_filtered_30 <- methyl_30[coverage > 10, ..select_cols]

methyl_45 <- fread("/home/mbrothers/nanopore/210310_Russula/modified_bases.aggregate04.6mA.bed",
                  col.names = columns)
methyl_filtered_45 <- methyl_45[coverage > 10, ..select_cols]

methyl_60 <- fread("/home/mbrothers/nanopore/210310_Russula/modified_bases.aggregate05.6mA.bed",
                   col.names = columns)
methyl_filtered_60 <- methyl_60[coverage > 10, ..select_cols]

methyl_90 <- fread("/home/mbrothers/nanopore/210310_Russula/modified_bases.aggregate06.6mA.bed",
                    col.names = columns)
methyl_filtered_90 <- methyl_90[coverage > 10, ..select_cols]

methyl_ss <- fread("/home/mbrothers/nanopore/210205_Ocular/modified_bases.aggregate08.6mA.bed",
                    col.names = columns)
methyl_filtered_ss <- methyl_ss[coverage > 10, ..select_cols]

# PLOT CONTROL
plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, "III", 90e3, 95e3)

# PLOT HML
plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, "III", 10e3, 15.5e3)

# with genes:
# abline(v = c(11146, 14849, 11082, 13282, 13809, 12386, 13018))

# PLOT HMR
plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, "III", 292e3, 295.5e3)

# with genes:
# abline(v = c(292674, 292769, 294805, 294864, 293179, 293538, 293835, 294321))

# PLOT TELOMERES
chromosomes <- c("I", "II", "III", "IV", "V", "VI", "VII", "VIII",
                 "IX", "X", "XI", "XII", "XIII", "XIV", "XV", "XVI")
right_ends <- c(230218, 813184, 316620, 1531933, 576874, 270161, 1090940, 562643,
                439888, 745751, 666816, 1078177, 917000, 765000, 1091291, 948010)
for(i in 1:16){
  plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, chromosomes[i], 0, 10e3)
}

for(i in 1:16){
  plot_timecourse(methyl_0, methyl_15, methyl_30, methyl_45, methyl_60, methyl_90, methyl_ss, chromosomes[i], right_ends[i]-10e3, right_ends[i])
}