QC_Pipeline/03_bkgdcor_beta_BMIQ_bySample.R at master · PGC-PTSD-EWAS/QC_Pipeline · GitHub

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
56
57
58
59
60
61
62
63
64
####################################################################################
# Background Corrected: Removing Cross-hybridized probes and running BMIQ
####################################################################################

rm(list=ls())
library(wateRmelon)

####################################################################################
# Step 1: Loading beta matrix, annotation, and CpG sites to remove
####################################################################################

# 450K annotation for probe design
load("fullannotInd.rda")
workspace<-ls()
annot<-fullannot
rm(list=workspace, "workspace")

# CpGassoc QC'ed beta matrix
load("bkgdcor_beta_QC.Rdata")

# Cross-hybridized probes from Chen 2013
chen<-read.table("48639-non-specific-probes-Illumina450k.txt", header=T, sep="\t")

####################################################################################
# Step 2: Removing SNP-associated and Cross-hybridized probes
####################################################################################

# Creating list of probes to remove
probesRmv<-as.character(chen[, "TargetID"]) # list of probes to remove
probesRmv<-probesRmv[!is.na(match(probesRmv, rownames(beta.qc)))]
beta.qc<-beta.qc[-match(probesRmv, rownames(beta.qc)),]
rm(chen, probesRmv)

####################################################################################
# Step 3: Creating design.v. for BMIQ normalization.
####################################################################################

annot<-annot[rownames(beta.qc),]
beta.qc[which(beta.qc==0)]<-0.000001
beta.qc[which(beta.qc==1)]<-0.999999

design.v<-annot[, "INFINIUM_DESIGN_TYPE"]
design.v[design.v=="I"]<-1
design.v[design.v=="II"]<-2
design.v<-as.numeric(design.v)

####################################################################################
# Step 4: Running BMIQ
####################################################################################

# Checking order of beta and design matrices before running BMIQ
for(ii in 1:ncol(beta.qc)){
  beta.v<-beta.qc[, ii]
  sID<-colnames(beta.qc)[ii]
  bmiqTemp<-BMIQ(beta.v=beta.v, design.v=design.v, sampleID=sID)
  beta.qc[,ii]<-bmiqTemp$nbeta
}

bmiq<-beta.qc

# Save list of all 8 outputs from BMIQ function
save(bmiq, design.v, file="bkgd_beta_BMIQ_bySample.Rdata")

rm(list=ls()[-match(c("bmiq", "design.v"), ls())])