#1. MCPcounter
# library(devtools)
# install_github("ebecht/MCPcounter",ref="master", subdir="Source")
rm(list=ls())
library(tibble)
library(dplyr)
library(tidyr)
library(MCPcounter)
setwd("C:\\Users\\Administrator\\Desktop\\MCPcounter")
cdata <- read.table('uniq.symbol.txt', sep = "\t", header = T,  row.names = 1)

probesets <- read.table("probesets.txt", sep = "\t",
                        stringsAsFactors = FALSE, colClasses = "character")

genes <- read.table("genes.txt", sep = "\t", stringsAsFactors = FALSE,
                    header=TRUE, colClasses="character", check.names=FALSE)

featuresType <- c("affy133P2_probesets","HUGO_symbols","ENTREZ_ID")[2]


results_MCPcounter <- MCPcounter.estimate(cdata, featuresType,
                              probesets = probesets,
                              genes = genes)
results_MCPcounter <- t(results_MCPcounter) 
dim(results_MCPcounter)


results_MCPcounter <- rbind(ID = colnames(results_MCPcounter), results_MCPcounter)
write.table(results_MCPcounter, "results_MCPcounter.txt", sep = "\t", quote = F, col.names = F)

#2. ESTIMATE
#library(utils)
#rforge <- "http://r-forge.r-project.org"
#install.packages("estimate", repos=rforge, dependencies=TRUE)
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.8")
#BiocManager::install("estimate")
library(limma)
library(estimate)
setwd("C:\\Users\\Administrator\\Desktop\\METTL7A_Immune\\2. ESTIMATE")          
inputFile="Merge_Matrix.TPM.txt"                                             

rt=read.table(inputFile,sep="\t",header=T,check.names=F)
dim(rt)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

group_rep <- sapply(strsplit(colnames(data),"\\_"), "[", 1)
length(group_rep) 
group_rep1 <- group_rep[!duplicated(group_rep)]
length(group_rep1) 
data <- data[, group_rep1]
dim(data
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)  
group=gsub("2","1",group)
data=data[,group==0]
dim(data)
data_normal <- data[, group == 1]
data_tumor <- data[, group == 0]
#合并. data <- cbind(data_tumor, data_normal)
out=data[rowMeans(data)>0,]
out=rbind(ID=colnames(out),out)
options(scipen = 200) 
out <- out[-1,]
dim(out)  
str(out) 
colnames_out <- colnames(out)
out1 <- apply(out, 2, as.numeric)  
str(out1)
dim(out1) 
rownames(out1) <- rownames(out)
out2 <- log2(out1 + 1) 
out2 <- rbind(ID=colnames(out2),out2)

write.table(out2,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)


filterCommonGenes(input.f="uniq.symbol.txt", 
                  output.f="commonGenes.gct", 
                  id="GeneSymbol")

estimateScore(input.ds = "commonGenes.gct",
              output.ds="estimateScore.gct", 
              platform="illumina")


scores=read.table("estimateScore.gct",skip = 2,header = T)
rownames(scores)=scores[,1]
scores=t(scores[,3:ncol(scores)])
rownames(scores)=gsub("\\.","\\-",rownames(scores))
out=rbind(ID=colnames(scores),scores)
write.table(out,file="scores.txt",sep="\t",quote=F,col.names=F)

#3.quantiseq

#library(devtools)
#Sys.setlocale(category = "LC_ALL", locale = "us")
#install_github("icbi-lab/immunedeconv")
rm(list = ls())
library(immunedeconv)
library(ggplot2)
library(tidyverse)
setwd("C:\\Users\\Administrator\\Desktop\\")
df <- read.table("uniq.symbol.txt", sep = "\t", header = T, row.names = 1)

res <- deconvolute(df, method = "quantiseq")
write.table(res, "quantiseq.txt", sep = "\t", col.names = T, row.names = F, quote = F)


pdf("quantiseq.pdf")
res %>%
  gather(sample, fraction, -cell_type) %>%

  ggplot(aes(x=sample, y=fraction, fill=cell_type)) +
  geom_bar(stat='identity') +
  coord_flip() +
  scale_fill_brewer(palette="Paired") +
  scale_x_discrete(limits = rev(levels(res)))
dev.off()

#4.EPIC

#library(devtools)
#devtools::install_github("GfellerLab/EPIC", build_vignettes=TRUE)
rm(list = ls())
library(EPIC) 

setwd("C:\\Users\\Administrator\\Desktop\\EPIC")
data <- read.table("uniq.symbol.txt", sep = "\t", header = T, row.names = 1)


EPIC_Results <- EPIC(data)


EPIC_Results_out <- EPIC_Results[["mRNAProportions"]]


EPIC_Results_out <- rbind(ID = colnames(EPIC_Results_out), EPIC_Results_out)
write.table(EPIC_Results_out, "EPIC_Results_out.txt", sep = "\t", quote = F, col.names = F)

#5.K_M analysis
rm(list = ls())
library(survival)
library(survminer)
setwd("C:\\Users\\Administrator\\Desktop\\
res.cut <- surv_cutpoint(svdata, time = "time", event = "event", variables = names(svdata)[3:ncol(svdata)], minprop = 0.3)#
res.cat <- surv_categorize(res.cut)
my.surv <- Surv(res.cat$time, res.cat$event)
pl<-list()          
for (i in colnames(res.cat)[3:ncol(svdata)]) {
  group <- res.cat[,i] 
  survival_dat <- data.frame(group = group)
  fit <- survfit(my.surv ~ group)
  group <- factor(group, levels = c("low", "high"))
  data.survdiff <- survdiff(my.surv ~ group)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
  up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
  low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
  HR <- paste("Hazard Ratio = ", round(HR,2), sep = "")
  CI <- paste("95% CI: ", paste(round(low95,2), round(up95,2), sep = " - "), sep = "")
  svsort <- svdata[order(svdata[,i]),]
  pl[[i]]<-ggsurvplot(fit, data = survival_dat ,conf.int = F, censor = F, palette = c("#E7B800","#00AFBB"), legend.title = i,font.legend = 11,legend.labs=c(paste0(">",round(svsort[fit$n[2],i],2),"(",fit$n[1],")"),paste0("<",round(svsort[fit$n[2],i],2),"(",fit$n[2],")")),pval = paste(pval = ifelse(p.val < 0.001, "p < 0.001", paste("p = ",round(p.val,3), sep = "")),HR, CI, sep = "\n"))}
res <- arrange_ggsurvplots(pl, print = T,ncol = 3, nrow = 4)
ggsave("OS.pdf",res,width = 12,height = 16)

#6.pROC
#install.packages("pROC")
rm(list = ls())
library(pROC)             
inputFile="input.txt"      
outFile="ROC.pdf"         
x="RiskScore"                
setwd("C:\\Users\\Administrator\\Desktop\\")   

rt=read.table(inputFile,header=T,sep="\t",check.names=F,row.names=1)
y=colnames(rt)[1]

rocobj1=roc(rt[,y], as.vector(rt[,x]))
pdf(file=outFile,width=5,height=5)
plot(rocobj1, print.auc=TRUE, col="red")
dev.off()