library(methylKit)

setwd("~/Projects/BullsRRBS_Heat_Stress_AUG2016/bismark_aligned/")
list.files()

# Get list of files. Here we will read from coverage file output by bismark
files <- list(
                "sample_7307_C1_bismark_out/7307_C1_R1_val_1_bismark_pe.bismark.cov",
                "sample_7334_C1_bismark_out/7334_C1_R1_val_1_bismark_pe.bismark.cov",
                "sample_8062_C1_bismark_out/8062_C1_R1_val_1_bismark_pe.bismark.cov",
                "sample_7307_T8_bismark_out/7307_T8_R1_val_1_bismark_pe.bismark.cov",
                "sample_7334_T8_bismark_out/7334_T8_R1_val_1_bismark_pe.bismark.cov",
                "sample_8062_T8_bismark_out/8062_T8_R1_val_1_bismark_pe.bismark.cov"
              )

myobj <- methylKit::methRead(files,
               sample.id=list("S7307_C1", "S7334_C1", "S8062_C1",
                              "S7307_T8", "S7334_T8", "S8062_T8"),
               treatment=c(0, 0, 0, 1, 1, 1),
               context="CpG",
               pipeline = "bismarkCoverage",
               assembly = "UMD3.1")

# Plot methylation percentage histograms for each sample
getMethylationStats(myobj[[1]],plot=TRUE,both.strands=FALSE)
getMethylationStats(myobj[[2]],plot=TRUE,both.strands=FALSE)
getMethylationStats(myobj[[3]],plot=TRUE,both.strands=FALSE)
getMethylationStats(myobj[[4]],plot=TRUE,both.strands=FALSE)
getMethylationStats(myobj[[5]],plot=TRUE,both.strands=FALSE)
getMethylationStats(myobj[[6]],plot=TRUE,both.strands=FALSE)

# Plot read coverage for each of the samples
getCoverageStats(myobj[[1]],plot=TRUE,both.strands=FALSE)
getCoverageStats(myobj[[2]],plot=TRUE,both.strands=FALSE)
getCoverageStats(myobj[[3]],plot=TRUE,both.strands=FALSE)
getCoverageStats(myobj[[4]],plot=TRUE,both.strands=FALSE)
getCoverageStats(myobj[[5]],plot=TRUE,both.strands=FALSE)
getCoverageStats(myobj[[6]],plot=TRUE,both.strands=FALSE)

# Remove CpG sites with low coverage (below 10 reads) and extremely high coverage > 99.9 percentile
filtered.myobj <- filterByCoverage(myobj,lo.count=10,lo.perc=NULL,
                                hi.count=NULL,hi.perc=99.9)
norm.myobj <- normalizeCoverage(filtered.myobj, method="median")
# Get CpGs covered in all samples
meth <- unite(norm.myobj, destrand=FALSE)
meth
# We are only comparing 7705 CpGs between Control and Day 8

# Plot correlations and scatter plots between samples
getCorrelation(meth,plot=TRUE)

# Cluster samples based on correlation. Method "Ward"
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)

# PCA plot
PCASamples(meth)

#####################################################################################################################
## Get differentially methylated CpGs
myDiff <- calculateDiffMeth(meth, num.cores = 4)

# Get hyper- and hypo-methylated CpGs
# get hyper methylated bases
myDiff25p.hyper <- getMethylDiff(myDiff,difference=20,qvalue=0.05,type="hyper") # 17 Hyper CpGs
myDiff25p.hyper
# get hypo methylated bases
myDiff25p.hypo <- getMethylDiff(myDiff,difference=20,qvalue=0.05,type="hypo") # 3 Hypo CpGs
myDiff25p.hypo
# get all differentially methylated bases
myDiff25p <- getMethylDiff(myDiff,difference=20,qvalue=0.05) # 20 differentially methylated
myDiff25p