# The following search string was used in the search field of ArrayExpress:

# ef:"Sex" raw:true samplecount:[30 TO 500] organism:Homo sapiens transcription HG-U133_Plus_2



# Accession numbers of first 39 search hits:

accnumber <- c("E-GEOD-60542", "E-GEOD-71571", "E-GEOD-64415", "E-GEOD-48278", 
  "E-GEOD-58697", "E-GEOD-68720", "E-GEOD-67851", "E-GEOD-64857",
  "E-GEOD-53224", "E-GEOD-49697", "E-GEOD-49695", "E-GEOD-28654",
  "E-GEOD-54543", "E-GEOD-54495", "E-GEOD-63626", "E-GEOD-57065",
  "E-GEOD-61804", "E-GEOD-33532", "E-GEOD-19722", "E-GEOD-59312",
  "E-GEOD-25504", "E-GEOD-39084", "E-GEOD-60028", "E-GEOD-45670",
  "E-GEOD-38662", "E-GEOD-58558", "E-GEOD-49243", "E-GEOD-54837",
  "E-GEOD-50774", "E-GEOD-46474", "E-GEOD-48780", "E-GEOD-48350",
  "E-GEOD-53890", "E-GEOD-51010", "E-GEOD-45216", "E-GEOD-53987",
  "E-GEOD-29623", "E-GEOD-29621", "E-GEOD-45867")



# Links to the corresponding sdrf-files:

sdrffiles <- paste("http://www.ebi.ac.uk/arrayexpress/files/", accnumber, 
  "/", accnumber, ".sdrf.txt", sep="")



# Make a folder in which the phenodata files will be stored:

sdrffolder <- "./ICSVCodeData/DataRetrieval/sdrffiles/"

if(!file.exists(sdrffolder))
  dir.create(sdrffolder)



# Function to read in an sdrf-file and store the content in 
# a data.frame:
  
readsdrf <- function(sdrf) {

adf <- readLines(sdrf)

if(length(unique(sapply(adf, function(x) length(strsplit(x, split="\t")[[1]])))) > 1)
  stop("Number of entries was not the same for all rows.")

datainfo <- data.frame(matrix(nrow=0, ncol=length(strsplit(adf[1], split="\t")[[1]])))
names(datainfo) <- strsplit(adf[1], split="\t")[[1]]

for(i in 1:(length(adf)-1))
  datainfo[i,] <- strsplit(adf[i+1], split="\t")[[1]]
rownames(datainfo) <- NULL

rownames(datainfo) <- datainfo$'Source Name'

return(datainfo)

}



# Download and read in all sdrf-files and each time extract
# the "sex" variable:

# Initialization of list where the "sex" variable for each dataset
# will be stored:
sexlist <- list()


for(i in seq(along=sdrffiles)) {

# Download sdrf-file:
download.file(url=sdrffiles[i], destfile=paste(sdrffolder, accnumber[i], ".sdrf.txt", sep=""))

# Read in sdrf-file:
datainfo <- readsdrf(sdrffiles[i])

if(length(datainfo$'Array Design REF')==0)
  stop("No array design variable.")

# Subset to feature only samples from the HG-U133_Plus_2-microarray:
datainfo <- datainfo[datainfo$'Array Design REF'=="A-AFFY-44",]

if(nrow(datainfo)==0)
  stop("Only samples from wrong array type.")

# Index of the "sex" variable in the pheno data:
sexind <- which(names(datainfo) %in% c("Sex", "Characteristics [Sex]", "Characteristics [sex]"))

if(length(sexind)==0)
  stop("No sex variable.")

# In case there were more than one "sex" variables found, use the first:
sexind <- sexind[1]  

# Store the "sex"-variable in the list initiated above:
eval(parse(text=paste("sexlist[[i]] <- datainfo$'", names(datainfo)[sexind], "'", sep="")))

# Add uniform "sex"-variable to pheno data frame:
datainfo$sexvar <- sexlist[[i]]

# Save subset phenodata as Rda-file:
save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

}




# The "sex" variable in one dataset uses different labels than "male" and "female":

corruptinds <- which(sapply(sexlist, function(x) !(("male" %in% x) | ("female" %in% x))))
sexlist[corruptinds]

# --> Change "46,XX" to "female" and "46,XY" to "male":

i <- corruptinds[1]

# Save subset phenodata:
load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

sexlist[[i]][sexlist[[i]]=="46,XX"] <- "female"
sexlist[[i]][sexlist[[i]]=="46,XY"] <- "male"

datainfo$sexvar <- sexlist[[i]]

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))





# For three datasets information on sex was missing for some patients.
# --> Exclude these patients:


corruptinds <- which(sapply(sexlist, function(x) any(!(x %in% c("male", "female")))))


i <- corruptinds[1]

unique(sexlist[[i]])

# --> Exclude observations where 'sex' is missing:

# Save subset phenodata:
load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

exclind <- which(sexlist[[i]]=="  ")

sexlist[[i]] <- sexlist[[i]][-exclind]
datainfo <- datainfo[-exclind,]

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))



i <- corruptinds[2]

unique(sexlist[[i]])

# --> Exclude observations where 'sex' is missing:

# Save subset phenodata:
load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

exclind <- which(sexlist[[i]]=="  ")

sexlist[[i]] <- sexlist[[i]][-exclind]
datainfo <- datainfo[-exclind,]

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))



i <- corruptinds[3]

unique(sexlist[[i]])

# --> Exclude observations where 'sex' is missing:

# Save subset phenodata:
load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

exclind <- which(sexlist[[i]]=="unknown")

sexlist[[i]] <- sexlist[[i]][-exclind]
datainfo <- datainfo[-exclind,]

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))






# Some datasets contain repeated measurements. In these cases we randomly draw one measurement
# per patient.


# In datasets with multiple measurements ID variables are used to discern the different patients.

# Therefore we iterate through all datasets to spot the ID variables:

IDvars <- list()

for(i in 1:length(sexlist)) {

  IDvarstemp <- c()
  load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

  # Automatically identify candidates for ID variables:
  IDvarstemp <- c(IDvarstemp, grep("atient", names(datainfo), value=TRUE))
  IDvarstemp <- c(IDvarstemp, grep("id", names(datainfo), value=TRUE))
  IDvarstemp <- c(IDvarstemp, grep("ID", names(datainfo), value=TRUE))
  IDvarstemp <- setdiff(IDvarstemp, c(grep("plateid", IDvarstemp, value=TRUE), 
    grep("idhmut", IDvarstemp, value=TRUE), "Hybridization Name"))

  IDvars[[i]] <- list(IDind=i, IDvarstemp=IDvarstemp)

}

IDvars


# Manually inspect candidates for ID variables:

IDvars[sapply(IDvars, function(x) length(x$IDvarstemp)>0)]



# --> Indices of datasets with ID variables:

indwithidvar <- c(4, 6, 13, 14, 23, 31, 32, 35)
IDvars[indwithidvar]

# Names of the ID variables:

namesidvars <- c("Characteristics [subject id]",
"Characteristics [sample_id]",
"Characteristics [subject id]",
"Characteristics [subject id]",
"Characteristics [patient_number]",
"Characteristics [individual]",
"Characteristics [individual]",
"Characteristics [patient number]")



# For each of the datasets with ID variable randomly sample one measurement per patient:


i <- indwithidvar[1]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are repeated measurements.


eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

set.seed(1234)

eval(parse(text=paste("indivar <- datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

subind <- sapply(unique(indivar), function(x) {
  candind <- which(indivar==x)
  if(length(candind)==1)
    sampind <- candind
  else
    sampind <- sample(candind, size=1)
  return(sampind)
})

sexlist[[i]] <- sexlist[[i]][subind]
datainfo <- datainfo[subind,]

eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))





i <- indwithidvar[2]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are no repeated measurements.





i <- indwithidvar[3]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are no repeated measurements.





i <- indwithidvar[4]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are no repeated measurements.




i <- indwithidvar[5]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are repeated measurements.


eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

set.seed(1234)

eval(parse(text=paste("indivar <- datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

subind <- sapply(unique(indivar), function(x) {
  candind <- which(indivar==x)
  if(length(candind)==1)
    sampind <- candind
  else
    sampind <- sample(candind, size=1)
  return(sampind)
})

sexlist[[i]] <- sexlist[[i]][subind]
datainfo <- datainfo[subind,]

eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))




i <- indwithidvar[6]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)

eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

# --> There are no repeated measurements (as verified by looking at
#     the corresponding ArrayExpress-site.





i <- indwithidvar[7]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are repeated measurements.


eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

set.seed(1234)

eval(parse(text=paste("indivar <- datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

subind <- sapply(unique(indivar), function(x) {
  candind <- which(indivar==x)
  if(length(candind)==1)
    sampind <- candind
  else
    sampind <- sample(candind, size=1)
  return(sampind)
})

sexlist[[i]] <- sexlist[[i]][subind]
datainfo <- datainfo[subind,]

eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))





i <- indwithidvar[8]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

eval(parse(text=paste("length(unique(datainfo$'", 
  namesidvars[which(indwithidvar==i)], "'))", sep="")))
nrow(datainfo)
  
# --> There are repeated measurements.


eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

set.seed(1234)

eval(parse(text=paste("indivar <- datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

subind <- sapply(unique(indivar), function(x) {
  candind <- which(indivar==x)
  if(length(candind)==1)
    sampind <- candind
  else
    sampind <- sample(candind, size=1)
  return(sampind)
})

sexlist[[i]] <- sexlist[[i]][subind]
datainfo <- datainfo[subind,]

eval(parse(text=paste("datainfo$'", namesidvars[which(indwithidvar==i)], "'", sep="")))

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))






# In the following we spot those datasets for which the CEL-files are not given for each patient
# and subset these to feature only patients for which the CEL-files are given:


corruptarraydatafile <- FALSE

for(i in 1:length(sexlist)) {

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

corruptarraydatafile[i] <- !all((1:nrow(datainfo)) %in% c(grep("\\.cel", datainfo$'Array Data File'), grep("\\.CEL", datainfo$'Array Data File')))

}


corruptind <- which(corruptarraydatafile)

# --> There is only one dataset where there is not a CEL-file for each patient.


i <- corruptind[1]

accnumber[i]

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

datainfo$'Array Data File'

includeind <- which(datainfo$'Array Data File'!="  ")

sexlist[[i]] <- sexlist[[i]][includeind]
datainfo <- datainfo[includeind,]

save(datainfo, file=paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))






# Exclude datasets which feature any samples from other datasets:

allrownames <- c()
datasetind <- c()

for(i in 1:length(sexlist)) {

load(paste("./ICSVCodeData/DataRetrieval/sdrffiles/", accnumber[i], "_datainfo.Rda", sep=""))

allrownames <- c(allrownames, rownames(datainfo))
datasetind <- c(datasetind, rep(i, nrow(datainfo)))

}

datasetswithduplicates <- as.numeric(names(table(datasetind[duplicated(allrownames)])))

includeind <- setdiff(1:length(sexlist), datasetswithduplicates)
sexlist <- sexlist[includeind]
accnumber <- accnumber[includeind]




# Exclude datasets with less than 20 obsevations:

sort(sapply(sexlist, length))

includeind <- which(sapply(sexlist, length) >= 20)
sexlist <- sexlist[includeind]
accnumber <- accnumber[includeind]


accnumber



# ---> Accession numbers of remaining datasets to be used in the analysis:

datasets <- c("E-GEOD-60542", 
 "E-GEOD-71571", 
 "E-GEOD-64415", 
 "E-GEOD-48278", 
 "E-GEOD-58697", 
 "E-GEOD-68720", 
 "E-GEOD-67851", 
 "E-GEOD-64857", 
 "E-GEOD-53224", 
 "E-GEOD-49697", 
 "E-GEOD-28654", 
 "E-GEOD-54543", 
 "E-GEOD-63626", 
 "E-GEOD-57065", 
 "E-GEOD-61804", 
 "E-GEOD-33532", 
 "E-GEOD-19722", 
 "E-GEOD-59312", 
 "E-GEOD-39084", 
 "E-GEOD-60028", 
 "E-GEOD-45670", 
 "E-GEOD-38662", 
 "E-GEOD-58558", 
 "E-GEOD-49243", 
 "E-GEOD-54837", 
 "E-GEOD-50774", 
 "E-GEOD-46474", 
 "E-GEOD-48780", 
 "E-GEOD-48350", 
 "E-GEOD-53890", 
 "E-GEOD-45216", 
 "E-GEOD-53987", 
 "E-GEOD-29623", 
 "E-GEOD-45867")

