library(minfi)
library(minfiData)
library(shiny)
source("/data/minfiLab/extractFromRGSet450k.R")


#extractedDataMinfiLab <- extractFromRGSet450k(RGSet, file = "extractedDataMinfiLab.Rda")



### To read the data in from the IDAT files:
baseDir <- system.file("extdata",package="minfiData")
targets <- read.450k.sheet(baseDir)
head(targets)
RGSet <- read.450k.exp(base = baseDir, targets = targets)
RGSet

### To extract the phenotype data:
phenoData <- pData(RGSet)
phenoData[,1:6]


### Manifest
manifest <- getManifest(RGSet)
manifest
head(getProbeInfo(manifest))

### To get probe SNP probes:
snpProbesI <- manifest@data$TypeSnpI$Name
snpProbesII <- manifest@data$TypeSnpII$Name

### 
typeIProbes <- getProbeInfo(manifest, type = "I")$Name
typeIProbesGrn <- getProbeInfo(manifest, type = "I-Green")$Name
snpProbesI <- getProbeInfo(manifest, type = "SnpI")$Name
snpProbesII <- getProbeInfo(manifest, type = "SnpI")$Name
head(typeIProbes)
head(snpProbesII)

### Quality control - Detection p-values:
detP <- detectionP(RGSet) 
failed <- detP > 0.01 ### Failed if the detection pvalue is > 0.01
head(failed, n=3)

colMeans(failed) # Fraction of failed positions per sample
sum(rowMeans(failed)>0.5) # How many positions failed in more than 50% of the samples

### Converting the RGChannelSet to a MethylSet:
MSet <- preprocessRaw(RGSet) 
MSet
head(getMeth(MSet)[,1:3])
head(getUnmeth(MSet)[,1:3])

### To extract the quality control information:
qc <- getQC(MSet)
head(qc)
plotQC(qc)


## Density plots:
densityPlot(RGSet, sampGroups = phenoData$Sample_Group, main= "Beta", xlab = "Beta")
densityBeanPlot(RGSet, sampGroups = phenoData$Sample_Group)

## Control probe density plot:
controlStripPlot(RGSet, controls="BISULFITE CONVERSION II")

### Report:
qcReport(RGSet, pdf= "qcReport.pdf")

##########################################
######  shinyMethyl

source("/data/minfiLab/extractFromRGSet450k.R")
extractedDataMinfiLab <- extractFromRGSet450k(RGSet,
file = "extractedDataMinfiLab.Rda")

### For the Amazon machine:
source("/data/minfiLab/myRunApp.R")

### To run the app:
myRunApp(appDir = "/data/minfiLab/shinyMethyl")


### Preprocess Illumina: 
MSet.illumina <- preprocessIllumina(RGSet, bg.correct = TRUE,
                               normalize = "controls")
                               
### Preprocess SWAN:
MSet.swan <- preprocessSWAN(RGSet)


### To access Beta-values and M-values of a MethylSet:
getBeta(MSet, type = "Illumina")[1:4,1:3]
getM(MSet)[1:4,1:3]

### To create a RatioSet
ratioSet <- ratioConvert(MSet, what = "both", keepCN = TRUE)
ratioSet

### Map to genome

gRatioSet <- mapToGenome(ratioSet, genomeBuild = "hg19", mergeManifest = TRUE)
gRatioSet
showClass("GenomicRatioSet")


### Different accessors for the genomicRatioSet

getBeta(gRatioSet)
getM(gRatioSet)
getCN(gRatioSet)


sampleNames <- sampleNames(gRatioSet)
probeNames <- featureNames(gRatioSet)
pheno <- pData(gRatioSet)
gRanges <- rowData(gRatioSet)
head(gRanges, n= 2)

### Sex prediction:
predictedSex <- getSex(gRatioSet, cutoff = -2)$predictedSex
head(predictedSex)

### To plot the predicted set
plotSex(getSex(gRatioSet, cutoff = -2))

### preprocessQuantile

gRatioSet.quantile <- preprocessQuantile(RGSet, fixOutliers = TRUE,
removeBadSamples = TRUE, badSampleCutoff = 10.5,
quantileNormalize = TRUE, stratified = TRUE, 
mergeManifest = FALSE, sex = NULL)

### dmpFinder
beta <- getBeta(gRatioSet.quantile)
predictedSex <- pData(gRatioSet.quantile)$predictedSex
dmp <- dmpFinder(beta, pheno = predictedSex  , type = "categorical")
head(dmp)


### bumphunter
pheno <- pData(gRatioSet.quantile)$age
designMatrix <- model.matrix(~ pheno)
dmrs <- bumphunter(gRatioSet.quantile, design = designMatrix, cutoff = 0.05, B=1)
names(dmrs)
head(dmrs$table, n =3)




### Exercices:

# 1. Before processing a RGChannelSet further, coudl you remove the probes which failed more than 50% of the samples in the example dataset?

# 2. For the top loci that we found differentially methylated for the predicted sex, could you tell if those loci are mostly mapped to the X and Y chromosomes? 

# 3 It is known that the Beta-value distribution of the Type I probes is different from the Beta value distribution of the Type II probes. Can you verify this with by plotting the Beta-value distribution density for each type separately?