Edit November 28, 2020:
Further reproducible code: A: GPL6883_HumanRef-8_V3_0_R0_11282963_A (illumina expression beadchip)
--
Most Illumina 'chip' studies that I have seen on GEO do not contain the raw data IDAT files. You can start with the tab-delimited file, but will also require the annotation file (contained in the *_RAW.tar file), and the usual targets file for limma.
The targets file may look like this:
IDATfile Group
raw/EX249_001.idat Peripheral_Blood_cDC1
raw/EX249_005.idat Peripheral_Blood_cDC2
raw/EX249_007.idat Peripheral_Blood_early_pre_DC
...
Please! always ensure that this is aligned correctly with the data with which you are working. Some of these functions have no internal checks to ensure that the targets file data is aligned to the expression data. Check it regularly.
# general config
baseDir <- '/home/kblighe/Escritorio/GSE20159/'
bgxfile <- 'Annot/GPL6947_HumanHT-12_V3_0_R1_11283641_A.bgx.gz'
targetsfile <- 'targets.txt'
setwd(baseDir)
options(scipen = 99)
require(limma)
# read in the data and convert the data to an EListRaw object
x <- read.table(
paste0(baseDir, 'raw/GSE80171_non_normalized.txt'),
header = TRUE, sep = '\t', stringsAsFactors = FALSE, skip = 0)
# extract detection p-value columns
detectionpvalues <- x[,grep('Detection.Pval', colnames(x))]
x <- x[,-grep('Detection.Pval', colnames(x))]
# set rownames and tidy up final expression matrix
probes <- x$ID_REF
x <- data.matrix(x[,2:ncol(x)])
rownames(x) <- probes
colnames(x) <- gsub('^X', '',
gsub('\\.AVG_Signal', '', colnames(x)))
# rename samples manually (optional)
...
# read in annotation and align it with the expression data
annot <- illuminaio::readBGX(bgxfile)$probes
annot <- annot[,which(colnames(annot) %in% c('Source','Symbol','Transcript','ILMN_Gene','RefSeq_ID',
'Entrez_Gene_ID','Symbol','Protein_Product','Probe_Id','Probe_Type',
'Probe_Start','Chromosome','Probe_Chr_Orientation','Probe_Coordinates',
'Cytoband', 'Definition', 'Ontology_Component', 'Ontology_Process',
'Ontology_Function', 'Synonyms'))]
annot <- annot[which(annot$Probe_Id %in% rownames(x)),]
annot <- annot[match(rownames(x), annot$Probe_Id),]
# update the target info
targetinfo <- readTargets(targetsfile, sep = '\t')
rownames(targetinfo) <- gsub('\\.idat$', '',
gsub('^raw/', '', targetinfo$IDATfile))
x <- x[,match(rownames(targetinfo), colnames(x))]
if (!all(colnames(x) == rownames(targetinfo)))
stop('Target info is not aligned to expression data - they must be in the same order')
# create a custom EListRaw object
project <- new('EListRaw')
project@.Data[[1]] <- 'illumina'
project@.Data[[2]] <- targetinfo
project@.Data[[3]] <- NULL
project@.Data[[4]] <- x
project@.Data[[5]] <- NULL
project$E <- x
project$targets <- targetinfo
project$genes <- NULL
project$other$Detection <- detectionpvalues
# generate QC plots of raw intensities
dir.create('QC/')
...
# for BeadArrays, background correction and normalisation are handled by a single function: neqc()
project.bgcorrect.norm <- neqc(project, offset = 16)
# filter out control probes, those with no symbol, and those that failed
annot <- annot[which(annot$Probe_Id %in% rownames(project.bgcorrect.norm)),]
project.bgcorrect.norm <- project.bgcorrect.norm[which(rownames(project.bgcorrect.norm) %in% annot$Probe_Id),]
annot <- annot[match(rownames(project.bgcorrect.norm), annot$Probe_Id),]
project.bgcorrect.norm@.Data[[3]] <- annot
project.bgcorrect.norm$genes <- annot
Control <- project.bgcorrect.norm$genes$Source=="ILMN_Controls"
NoSymbol <- project.bgcorrect.norm$genes$Symbol == ""
isexpr <- rowSums(project.bgcorrect.norm$other$Detection <= 0.05) >= 3
project.bgcorrect.norm.filt <- project.bgcorrect.norm[!Control & !NoSymbol & isexpr, ]
dim(project.bgcorrect.norm)
dim(project.bgcorrect.norm.filt)
# remove annotation columns we no longer need
project.bgcorrect.norm.filt$genes <- project.bgcorrect.norm.filt$genes[,c(
'Probe_Id',
'Definition','Ontology_Component','Ontology_Process','Ontology_Function',
'Chromosome','Probe_Coordinates','Cytoband','Probe_Chr_Orientation',
'RefSeq_ID','Entrez_Gene_ID','Symbol')]
head(project.bgcorrect.norm.filt$genes)
# summarise across genes by mean - use ID
project.bgcorrect.norm.filt.mean <- avereps(project.bgcorrect.norm.filt,
ID = project.bgcorrect.norm.filt$genes$Symbol)
