###NASTAVENI R
#nacteni knihoven
library(oligo)
library(affycoretools)
library(pd.clariom.d.human)
library(clariomdhumanprobeset.db)
#nastaveni cesty
cel_path = "Projekt/"

### NACTENI DAT
cel_files = list.celfiles(cel_path, full.names=TRUE)
raw = read.celfiles(cel_files)
str(raw)
s = sampleNames(raw)
s = substring(s, 1, nchar(s) - 4)
ns = length(s)

#ulozeni informaci o vzorcich
a = function(x) { r = strsplit(x, "_")[[1]]; r }
info = as.data.frame(t(data.frame(lapply(s, a))))
colnames(info) = c() ## nutno vlozit vlastni colnames!
rownames(info) = sampleNames(raw)
info = info[,1:3]
save(info, file="sample_info.rdata")
phenoData(raw) = new('AnnotatedDataFrame', 
                     data=info, 
                     varMetadata=data.frame(labelDescription = c("PatID", "blk", "morph"), 
                                           channel=factor('_ALL_')))

#ulozenie raw dat
save(raw, file="./raw_cels_to_use.rdata")

###KONTROLA KVALITY
#histogram
pdf('Histogram.pdf')
par(mfrow=c(1,1))
hist(raw, transfo=log2, which=c("all"), main='Expresie')
dev.off()

#MA plot
pdf('MAplot.pdf')
par(mfrow=c(1,1))
MAplot(raw)
dev.off()

#PLM
qc.plm = fitProbeLevelModel(raw)
save(qc.plm, file="./PLM_cels_to_use.rdata")

#RLE
qc.rle = RLE(qc.plm)
pdf('RLE.pdf')
RLE(qc.plm)
dev.off()

#NUSE
qc.nuse = NUSE(qc.plm)
abline(a=1.02,b=0, col='red')
median_nuse_values <- apply(qc.nuse, 2, function(x) median(x, na.rm = TRUE))
median_nuse_values[median_nuse_values > 1.02]

save(qc.nuse, qc.rle, file="./QC-metrics_cels_to_use.rdata")
pdf('NUSE.pdf')
NUSE(qc.plm)
abline(a=1.02,b=0, col='red')
dev.off()

#vykresleni obrazku cipu
par(mfrow=c(1,1))
for (i in c(seq(1,39))) {
  file_name <- paste(qc.plm@phenoData@data$PatID[i], ".pdf", sep="")
  pdf(file_name, width=10, height=10)
  oligo::image(qc.plm, which=i, type="residuals")
  dev.off()
}

###NORMALIZACE
#pomoci rma
detach("package:affycoretools", unload = TRUE)
normRMA = rma(raw, target = "probeset", background = T, normalize = T)
save(normRMA, file="./rma_cels_to_use.rdata")

#MAplot po normalizaci
pdf('MAplot_afterRMA.pdf')
par(mfrow=c(1,1))
MAplot(normRMA)
dev.off()

#extrakce expresi
X = t(exprs(normRMA)) 
dim(X)
X[1:5,1:5]

#nazvy genu z probeset IDs
library(affycoretools)
eset <- annotateEset(normRMA, pd.clariom.d.human, type = "probeset")
matched_names <- fData(eset)$SYMBOL[match(colnames(X), fData(eset)$PROBEID)]
colnames(X) <- matched_names
X[1:5,1:5]
na_columns <- which(is.na(colnames(X)))
X <- X[,-na_columns]
dim(X)
