# calculate mean intensity
library(GWASdata)
file <- system.file("extdata", "illumina_qxy.gds", package="GWASdata")
gds <- GdsIntensityReader(file)
data(illuminaScanADF)
intenData <- IntensityData(gds, scanAnnot=illuminaScanADF)
meanInten <- meanIntensityByScanChrom(intenData)
intenMatrix <- meanInten$mean.intensity
# find outliers
outliers <- list()
sex <- illuminaScanADF$sex
id <- illuminaScanADF$scanID
allequal(id, rownames(intenMatrix))
for (i in colnames(intenMatrix)) {
if (i != "X") {
imean <- intenMatrix[,i]
imin <- id[imean == min(imean)]
imax <- id[imean == max(imean)]
outliers[[i]] <- c(imin, imax)
} else {
idf <- id[sex == "F"]
fmean <- intenMatrix[sex == "F", i]
fmin <- idf[fmean == min(fmean)]
fmax <- idf[fmean == max(fmean)]
outliers[[i]][["F"]] <- c(fmin, fmax)
idm <- id[sex == "M"]
mmean <- intenMatrix[sex == "M", i]
mmin <- idm[mmean == min(mmean)]
mmax <- idm[mmean == max(mmean)]
outliers[[i]][["M"]] <- c(mmin, mmax)
}
}
par(mfrow=c(2,4))
intensityOutliersPlot(intenMatrix, sex, outliers)
close(intenData)
Run the code above in your browser using DataLab