04-tiles

library(methylKit)

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: BiocGenerics

## Warning: package 'BiocGenerics' was built under R version 4.0.5

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## Loading required package: GenomeInfoDb

## Warning: package 'GenomeInfoDb' was built under R version 4.0.5

cov.list=list('../data/1_ATCACG.CpG_report.merged_CpG_evidence.cov',
                 '../data/2_CGATGT.CpG_report.merged_CpG_evidence.cov',
                 '../data/3_TTAGGC.CpG_report.merged_CpG_evidence.cov',
                 '../data/4_TGACCA.CpG_report.merged_CpG_evidence.cov',
                 '../data/5_ACAGTG.CpG_report.merged_CpG_evidence.cov',
                 '../data/6_GCCAAT.CpG_report.merged_CpG_evidence.cov',
                 '../data/7_CAGATC.CpG_report.merged_CpG_evidence.cov',
                 '../data/8_ACTTGA.CpG_report.merged_CpG_evidence.cov'
)

myobj_lowCov = methRead(cov.list,
           sample.id = list("1","2","3","4","5","6","7","8"),
           assembly="v081",
           treatment = c(0,0,0,0,1,1,1,1),
           pipeline = "bismarkCoverage",
           context="CpG",
           mincov = 3
           )

## Received list of locations.

## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.

tiles = tileMethylCounts(myobj_lowCov,win.size=1000,step.size=1000,cov.bases = 10)
head(tiles[[1]],3)

#save(tiles, file = "../analyses/tiles_01")

load ("../analyses/tiles_01")

filtered.myobj=filterByCoverage(tiles,lo.count=5,lo.perc=NULL,
                                      hi.count=NULL,hi.perc=99.9)

meth_filter=unite(filtered.myobj, destrand=TRUE)

## uniting...

clusterSamples(meth_filter, dist="correlation", method="ward", plot=TRUE)

## The "ward" method has been renamed to "ward.D"; note new "ward.D2"

## 
## Call:
## hclust(d = d, method = HCLUST.METHODS[hclust.method])
## 
## Cluster method   : ward.D 
## Distance         : pearson 
## Number of objects: 8

PCASamples(meth_filter)

myDiff=calculateDiffMeth(meth_filter,mc.cores=8)

## two groups detected:
##  will calculate methylation difference as the difference of
## treatment (group: 1) - control (group: 0)

# get hyper methylated bases
myDiff25p.hyper=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hyper")

## Warning in max(i): no non-missing arguments to max; returning -Inf

#
# get hypo methylated bases
myDiff25p.hypo=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
#
#
# get all differentially methylated bases
myDiff10p=getMethylDiff(myDiff,difference=10,qvalue=.01)