library(affy)
library(limma)
##import phenotype data
phenoData = read.AnnotatedDataFrame('target.txt')
pheno = pData(phenoData)
View(pheno)
#import Annotaion
anno = read.csv("HG-U133_Plus_2-na36-annot.csv",head=T)
View(head(anno))
##RMA normalization
#eset.rma = RMA(Data)
eset.rma <- justRMA(filenames=paste(rownames(pheno),'.CEL',sep=''), celfile.path='./RAW')
datExpr = exprs(eset.rma)

library(impute)

imputed_gene_exp = impute.knn(datExpr,k=10,rowmax = 0.5,
                              colmax=0.8,maxp =3000, rng.seed=362436069)
datExpr2 = imputed_gene_exp$data

write.table(datExpr2,file="Expdata.txt",sep="\t")
library(sva)
batch <- pheno[ ,c('batch')]
batch
ComBatmodel <- model.matrix(~group, data = pheno)   
dataCombat <- ComBat(dat = datExpr2,batch = batch, 
                     mod = ComBatmodel, par.prior = T, prior.plot = T)
ExpMatrix <- dataCombat 


Group = factor(pheno$group,levels=c('Normal','Tumour'))
design = model.matrix(~0+Group)
colnames(design) <- c('Normal','Tumour')
design

fit <- lmFit(datExpr2, design)

contrast.matrix <- makeContrasts(Tumour-Normal,
                                 levels=design)

library(xlsx)
library(futile.logger)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.05,
                lfc=log(2,2),number=5000,sort.by = 'logFC')
DEG <- topTable(fit2, coef = 1, n = Inf, adjust.method = "fdr", 
                sort.by = 'logFC')     

 
DEG$change <- as.factor(ifelse(DEG$adj.P.Val < 0.05 & abs(DEG$logFC) > 1, 
                               ifelse(DEG$logFC > 1 ,'UP', 'DOWN'), 
                               'NOT'))
DEG$Gene = anno$Gene.Symbol[match(rownames(DEG),anno$ID_REF)]
DEG$ID_REF = rownames(DEG)
DEG = DEG[,c(8,7,1:6)]
DEG = DEG[DEG$Gene != '---',]
DEG <- DEG[!duplicated(DEG$Gene), ]
write.xlsx2(DEG,'DEG3.xlsx',sheetName=colnames(contrast.matrix)[1],col.names=T,row.names=F,append=T)

diff$Gene = anno$Gene.Symbol[match(rownames(diff),anno$ID_REF)]
diff$ID_REF = rownames(diff)
diff = diff[,c(8,7,1:6)]
diff = diff[diff$Gene != '---',]
diff <- diff[!duplicated(diff$Gene), ]

#output
write.xlsx2(diff,'DEG3.xlsx',sheetName=colnames(contrast.matrix)[1],col.names=T,row.names=F,append=T)
write.table(diff,file = 'DEG.txt', quote = F, na = "", row.names = F, sep = "\t")
 DEG$Gene = anno$Gene.Symbol[match(rownames(DEG),anno$ID_REF)]
 DEG$ID_REF = rownames(DEG)
 DEG = DEG[,c(8,7,1:6)]
 DEG = DEG[DEG$Gene != '---',]
 DEG <- DEG[!duplicated(DEG$Gene), ]
 write.table(DEG,file = 'DEG3.txt', quote = F, na = "", row.names = F, sep = "\t")
 library(ggplot2)
 titl <- paste0('Cutoff for logFC is ', '1',
                '\nThe number of up gene is ', 
                nrow(DEG[DEG$change == 'UP',]), 
                '\nThe number of down gene is ', 
                nrow(DEG[DEG$change == 'DOWN',]))
 ggplot(data = DEG, aes(x = logFC, y = -log10(P.Value), color = change)) +
   geom_point(alpha = 0.5, size = 1) +
   theme_set(theme_set(theme_bw(base_size = 20)))+
   xlab("log2(FC)") + 
   ylab("-log10(P value)") +
   ggtitle(titl) + 
   theme(plot.title = element_text(size = 15,hjust = 0.5))+
   scale_colour_manual(values = c('blue', 'grey10', 'red')) 
