-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path02.1_preprocessingSest.R
55 lines (31 loc) · 1.79 KB
/
02.1_preprocessingSest.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
47
48
49
50
51
52
53
54
55
###########################################################################################################################
# Author: Rebekka Koeck
# Lab: Cellular Genomic Medicine, Clinical Genetics, Maastricht University (Medical Centre +)
# script purpose: use minfi package to extract information from .idat files needed for sex prediction using SEst package
# input: data from each cohort / array type processed separately: IVF data (GSE196432), FLEHS 450k (GSE110128)
# output: detection p-values of all sites & samples (.csv)
# raw (unnormalised/unfiltered) beta values for all sites & samples (.csv)
###########################################################################################################################
# load the required packages
suppressPackageStartupMessages(library(minfi))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(data.table))
# set the data directories
base.dir <- "CultMed"
data.dir <- file.path(base.dir, "idat")
sample.annotation <- "annotation.txt"
report.dir <- "cohort1/sEST"
# reaad the raw data files using minfi
targets <- fread(sample.annotation) %>%
mutate(Basename = paste(Sentrix_ID, Sentrix_Position , sep = "_"))
RGSet <- read.metharray.exp(base = data.dir, targets = targets)
# proces the RGSet to an MSet which has beta values
MSet <- preprocessRaw(RGSet)
# calculate beta values by converting this to a ratioSet object
RSet <- ratioConvert(MSet, what = "beta")
# extract beta values and detection p values for all sites
detPval <- detectionP(RGSet)
betas <- getBeta(RSet)
# write the resultant files as csv files for further processing
write.csv(detPval, file = file.path(report.dir, "detPval.csv"), row.names = T)
write.csv(betas, file = file.path(report.dir, "betasRaw.csv"), row.names = T)