Description Usage Arguments Details Value Author(s) Examples
This function is for conveniently specifying a pipeline to preprocess methylation array data using common parameters from the minfi package and from common manipulations of experiment summary objects. It is recommended this be used to supplement, but not replace, effective analysis workflows, ie. to compare the efficacies of different normalization methods. Methylation experiment objects that have already been preprocessed can still be conveniently run through this workflow, ie. to map probes to the genome, convert to a new experiment object class, or conveniently filter on recognized lists of cross-reactive probes.
1 2 3 4 | methyPre(meExpObj, workflow = c("norm", "pfilt", "map", "minfifilt", "ufilt", "crxcgfilt", "batch"),
normfun = c("noob", "SWAN", "illumina", "funnorm"), pfilt = TRUE, detPcutoff = 0.05,
minfiFilt = TRUE, uFilt = NULL, crxcgFilt = c(NULL, "hm450", "epic"), batchCorrect = FALSE,
bcCovariate = NULL, bcBatch = NULL, returnObj = c("preObj1", "preObj2", "preObj3", "gexp1", "gexp2", "gexp3"))
|
meExpObj |
Experiment summary object containing methylation, annotation, and sample data. |
workflow |
User-specified options (as character strings) outlining workflow steps to be taken. All or part of the possible options may be used, with the implicit order of operations being: 1. norm (normalize, see normfun arg); 2. pfilt (intensity p-val filter, see pfilt and detPcutoff args); 3. map (map data to genome); 4. minfifilt (use minfi filtering functions and annotation, see minfiFilt arg); 5. ufilt (user-specified CpG IDs to filter, see uFilt); 6. crxcgfilt(filter low-quality probes, see crxcgFilt arg); 7. batch (perform ComBat batch correction, see batchCorrect, bcCovariate, and bcBatch variables). |
normfun |
Normalization function(s) to be used. User can specify any combination of options: Illumina; SWAN; noob; funnorm. The implicit order of operations is background (Illumina or noob) followed by quantile (SWAN or funnorm) normalizations. |
pfilt |
Whether to obtain intensity p-values for filtering. Requires that meExpObj be an RGChannelSet. |
detPcutoff |
Mean p-value cutoff for intensity p-value. ie. CpGs removed if mean of p-value across samples is above the cutoff. |
minfiFilt |
Whether to perform a series of filters using minfi functions and annotation: CpGs are filtered if they map to or are associated with SNPs, are CH/non-CpG probes, or map to X or Y chromosomes. |
uFilt |
User-specified list of CpGs to filter out. Should be a vector or list of characters strings. |
crxcgFilt |
Filter cross-reactive or poor-quality probes. Options are "hm450" or "epic", which filter either the HM450 CpGs specified in Chen et al 2013 (see ?chen_crxcg for info) or EPIC CpGs identified by Illumina (?illumina_crxcg) and Pidsley et al 2016 (?pidsley_crxcg) as being cross reactive, sensitive to genetic variation, or of otherwise poor quality. |
batchCorrect |
Should batch correction with ComBat be performed? |
bcCovariate |
Covariate variable (as vector or list) to be included in batch correction (ie. tissue type). |
bcBatch |
Variable (vector or list) used to specify batches for correction (ie. array plate ID). |
returnObj |
Objects to return, where gcombat is automatically returned if batch correction is specified. Options are: 1. preObj1: Illumina, noob, or unnormalized experiment object. 2. preObj2: preObj1 after SWAN, funnorm, or no additional normalization. 3. preObj3: preObj2 after detection-p value or no filter applied. 4. gexp1: preObj3 after mapping or no mapping performed. 5. gexp2: gexp1 after minfi CpG filters applied or not applied. 6. gexp3: gexp2 after user-defined and/or crxcg filters are applied or aren't applied. |
This is a wrapper function for common workflows to preprocess Illumina methylation arrays using pipelines in the minfi package.
Note, this function is intended enable not just apply a pipeline efficiently, but a convenient way to compare different workflows.
List of experiment summary objects containing methylation, annotation, and sample data.
(preObj1|preObj2|preObj3|gexp1|gexp2|gexp3|gcombat) |
Methylation experiment summary object. |
...
Sean Maden, smaden@fredhutch.org
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 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 | ## The function is currently defined as
function(meExpObj, workflow=c("norm","pfilt","map","minfifilt","ufilt","crxcgfilt","batch"),
normfun=c("noob","SWAN","illumina","funnorm"),
pfilt=TRUE, detPcutoff=0.05,minfiFilt=TRUE,uFilt=NULL,
crxcgFilt=c(NULL,"hm450","epic"), batchCorrect=FALSE,
bcCovariate=NULL,bcBatch=NULL,
returnObj=c("preObj1","preObj2","preObj3","gexp1","gexp2","gexp3")){
require(minfi)
return.list <- list()
if(length(workflow[grepl("norm",workflow)])>0){
if(length(normfun[grepl("illumina",normfun)])>0){
message("starting Illumina normalization..")
preObj1 <- preprocessIllumina(meExpObj)
message("finished Illumina normalization!")
} else if(length(normfun[grepl("noob",normfun)])>0){
message("starting Noob normalization..")
preObj1 <- preprocessNoob(meExpObj)
message("finished Noob normalization!")
} else{
message("no background normalizations performed. Continuing..")
preObj1 <- meExpObj
}
if(length(returnObj[grepl("preObj1",returnObj)])>0){
return.list <- append(return.list,list(preObj1))
}
if(length(normfun[grepl("swan",normfun)])>0){
message("beginning SWAN normalization..")
preObj2 <- preprocessSWAN(meExpObj,mSet=preObj1)
message("finished SWAN normalization!")
} else if(length(normfun[grepl("funnorm",normfun)])>0){
message("beginning Functional Normalization..")
preObj2 <- preprocessFunnorm(meExpObj,bgCorr=FALSE)
message("finished Functional Normalization!")
} else{
preObj2 <- preObj1
}
} else{
preObj2 <- meExpObj
}
if(length(returnObj[grepl("preObj2",returnObj)])>0){
return.list <- append(return.list,list(preObj2))
}
if(length(workflow[grepl("pfilt",workflow)])>0 & class(meExpObj)=="RGChannelSet"){
message("getting intensity p-values from RG set...")
detP <- detectionP(meExpObj)
message(paste0("filtering on mean detection cutoff =",detPcutoff))
failed <-rowMeans(detP) > detPcutoff
preObj3 <- preObj2[!failed,]
} else{
preObj3 <- preObj2
}
if(length(returnObj[grepl("preObj3",returnObj)])>0){
return.list <- append(return.list,list(preObj3))
}
if(length(workflow[grepl("map",workflow)])>0){
gexp1 <- mapToGenome(preObj3)
} else{
gexp1 <- preObj3
}
if(length(returnObj[grepl("gexp1",returnObj)])>0){
return.list <- append(return.list,list(gexp1))
}
if(length(workflow[grepl("minfifilt",workflow)])>0){
message("Applying minfi filters..")
message("Removing cg probes containing SNPs..")
gexp2 <- dropLociWithSnps(gexp1, snps=c("SBE", "CpG"))
message("Removing CH and SNP-assoc. probes...")
gexp2 <- dropMethylationLoci(gexp2,dropRS=TRUE,dropCH=TRUE)
message("Removing chrX and chrY-assoc. probes..")
gexp2 <- gexp2[!getAnnotation(gexp2)$chr
message("After applying minfi filters, ",nrow(gexp2)," CpGs remain.")
} else{
gexp2 <- gexp1
}
if(length(returnObj[grepl("gexp2",returnObj)])>0){
return.list <- append(return.list,list(gexp2))
}
if(length(workflow[grepl("ufilt",workflow)])>0|
length(workflow[grepl("crxcgfilt",workflow)])>0){
if(crxcgFilt=="hm450"){
message("Filtering cross-reactive HM450 probes...")
data(chen_crxcg)
gexp3 <- gexp2[!rownames(gexp2)
}
if(crxcgFilt=="epic"){
message("Filtering cross-reactive EPIC probes")
data(illumina_crxcg);data(pidsley_crxcg)
gexp3 <- gexp2[!rownames(gexp2)
}
if(!is.null(uFilt)){
message("Filtering user-specified cg list...")
gexp3 <- gexp2[!rownames(gexp2)
}
} else{
gexp3 <- gexp2
}
if(length(returnObj[grepl("gexp3",returnObj)])>0){
return.list <- append(return.list,list(gexp3))
}
if(batchCorrect){
bc.check <- readline(message("You selected to batch correct. Do you want to save now or continue? Enter 'save' or 'continue'"))
if(length(bc.check[grepl("save",bc.check)])>0){
return(gexp3)
}
if(length(bc.check[grepl("continue",bc.check)])>0){
require("sva")
if(!is.null(bcCovariate)){
message("Applying ComBat correction to M-values..")
combat_M <- ComBat(getM(gexp3),mod=model.matrix(~bcCovariate),batch=bcBatch)
message("Assembling corrected experiment object..")
gcombat <- GenomicRatioSet(gr=granges(gexp3),
Beta=ilogit2(combat_M),
M=combat_M,
annotation=annotation(gexp3),
preprocessMethod=preprocessMethod(gexp3),
pData=pData(gexp3))
return.list <- append(return.list,gcombat)
message("Completed all specified steps! Returning..")
return(return.list)
} else{
message("Applying ComBat correction to M-values..")
combat_M <- ComBat(getM(gexp3),batch=bcBatch)
message("Assembling corrected experiment object...")
gcombat <- GenomicRatioSet(gr=granges(gexp3),
Beta=ilogit2(combat_M),
M=combat_M,
annotation=annotation(gexp3),
preprocessMethod=preprocessMethod(gexp3),
pData=pData(gexp3))
return.list <- append(return.list,gcombat)
message("Completed all specified steps! Returning..")
return(return.list)
}
}
} else{
message("Completed all specified steps! Returning...")
return(return.list)
}
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.