-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathquality_control.R
46 lines (38 loc) · 1.9 KB
/
quality_control.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
############################################################################
############################################################################
########### ###########
########### Script for data Quality control for IDATs (RAW 450K) ###########
########### Author: Diego Montiel Gonzalez ###########
########### ###########
########### Erasmus MC University Medical Centre ###########
########### Rotterdam, The Netherlands ###########
########### ###########
########### [email protected] ###########
########### ###########
############################################################################
############################################################################
## Load libraries ##
library(minfi)
library(data.table)
library(wateRmelon)
library(ENmix)
library(RPMM)
setwd("PATH/Y-CpG/IDAT/")
## IDATs QC
# 0) Read IDATs in a rgSet object
rgSet <- read.metharray.exp(getwd(), extended = T)
qcinfo <- ENmix::QCinfo(rgSet)
# 1) Detect failed probes and samples with SQN w/rgSet
badsample <- unlist(lapply(X = strsplit(qcinfo$badsample, split = '[_]'), FUN = function(X) {X[1]})) # split the header names
length(qcinfo$badCpG)
write.csv(qcinfo$badsample, file="PATH/badsamples.txt", quote = F, row.names = F)
write.csv(qcinfo$badCpG, file="PATH/badCpGs.txt", quote = F, row.names = F)
# 2) Predict Sex on IDAT samples with Qantile Normalization
um.sqn <- preprocessQuantile(rgSet)
pdf(file = "sex_estimation_raw.pdf", width = 10, height = 7)
plotSex(um.sqn)
dev.off()
# 3) Plot density beta methylation values
pdf(file = "densityplot_raw.pdf", width = 10, height = 7)
densityPlot(rgSet)
dev.off()