library(limma)
library(ggplot2)
library(ggpubr)

gene="xxxx"              
expFile="xxx"     
setwd("xxx")    

rt=read.table(expFile, header=T, sep="\t", check.names=F)
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)
data=t(data[gene,,drop=F])

group=sapply(strsplit(rownames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
conNum=length(group[group==1])      
treatNum=length(group[group==0])     
Type=c(rep(1,conNum), rep(2,treatNum))

exp=cbind(data, Type)
exp=as.data.frame(exp)
colnames(exp)=c("gene", "Type")
exp$Type=ifelse(exp$Type==1, "Normal", "Tumor")
exp$gene=log2(exp$gene+1)

outTab=exp
colnames(outTab)=c(gene, "Type")
outTab=cbind(ID=row.names(outTab), outTab)
write.table(outTab, file="geneExp.txt", sep="\t", quote=F, row.names=F)

group=levels(factor(exp$Type))
exp$Type=factor(exp$Type, levels=group)
comp=combn(group,2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}

boxplot=ggboxplot(exp, x="Type", y="gene", color="Type",
                  xlab="",
                  ylab=paste0(gene, " expression"),
                  legend.title="Type",
                  palette = c("#88CEE6","#E89DA0"),
                  add = "jitter")+ 
  stat_compare_means(comparisons=my_comparisons,symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif")

pdf(file=paste0(gene,".diff.pdf"), width=5, height=4.5)
print(boxplot)
dev.off()


library(limma)
library(ggpubr)
expFile="xxx"     
setwd("xxxx")   

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

normalData=rt[rt$Type=="Normal",1,drop=F]
normalData=as.matrix(normalData)
rownames(normalData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(normalData))
normalData=avereps(normalData)
tumorData=rt[rt$Type=="Tumor",1,drop=F]
tumorData=as.matrix(tumorData)
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData))
tumorData=avereps(tumorData)
sameSample=intersect(row.names(normalData), row.names(tumorData))
data=cbind(normalData[sameSample,,drop=F], tumorData[sameSample,,drop=F])
colnames(data)=c("Normal", "Tumor")
data=as.data.frame(data)

pdf(file="pairDiff.pdf", width=5, height=4.5)
ggpaired(data, cond1="Normal", cond2="Tumor", fill="condition",
         xlab="", ylab=paste0(geneName, " expression"),
         legend.title="Type",
         palette=c("#88CEE6","#E89DA0"))+
  stat_compare_means(paired = TRUE, symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif",label.x = 1.35)
dev.off()

library(limma)
library(survival)
library(survminer)

expFile="xxx"    
cliFile="xxxx"      
setwd("xxx")    

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

tumorData=rt[rt$Type=="Tumor",1,drop=F]
tumorData=as.matrix(tumorData)
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData))
data=avereps(tumorData)


Type=ifelse(data[,gene]>median(data[,gene]), "High", "Low")
data=cbind(as.data.frame(data), Type)

cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
cli$futime=cli$futime/365

sameSample=intersect(row.names(data), row.names(cli))
data=data[sameSample,,drop=F]
cli=cli[sameSample,,drop=F]
rt=cbind(cli, data)

outTab=cbind(ID=row.names(rt), rt)
outTab=outTab[,-ncol(outTab)]
write.table(outTab, file="expTime.txt", sep="\t", quote=F, row.names=F)

diff=survdiff(Surv(futime, fustat) ~ Type, data=rt)
pValue=1-pchisq(diff$chisq, df=1)
if(pValue<0.001){
  pValue="p<0.001"
}else{
  pValue=paste0("p=", sprintf("%.03f",pValue))
}
fit <- survfit(Surv(futime, fustat) ~ Type, data = rt)
#print(surv_median(fit))

surPlot=ggsurvplot(fit, 
                   data=rt,
                   conf.int=F,
                   pval=pValue,
                   pval.size=6,
                   surv.median.line = "hv",
                   legend.title=gene,
                   legend.labs=c("High level", "Low level"),
                   xlab="Time(years)",
                   ylab="Overall survival",
                   break.time.by = 1,
                   palette=c("#E89DA0", "#88CEE6"),
                   risk.table=F,
                   risk.table.title="",
                   risk.table.col = "strata",
                   risk.table.height=.25)

pdf(file=paste0(gene, ".surv.pdf"), width=5.5, height=5, onefile=FALSE)
print(surPlot)
dev.off()


library(survival)
library(regplot)
library(rms)

expFile="xxx"     
cliFile="xxx"  

exp=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)

cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
cli=cli[apply(cli,1,function(x)any(is.na(match('unknow',x)))),,drop=F]
cli$Age=as.numeric(cli$Age)

samSample=intersect(row.names(exp), row.names(cli))
exp1=exp[samSample,,drop=F]
cli=cli[samSample,,drop=F]
rt=cbind(exp1, cli)

res.cox=coxph(Surv(futime, fustat) ~ . , data = rt)
nom1=regplot(res.cox,
             plots = c("density", "boxes"),
             clickable=F,
             title="",
             points=TRUE,
             droplines=TRUE,
             observation=rt[2,],
             rank="sd",
             failtime = c(1,3,5),
             prfail = F)

nomoRisk=predict(res.cox, data=rt, type="risk")
rt=cbind(exp1, Nomogram=nomoRisk)
outTab=rbind(ID=colnames(rt), rt)
write.table(outTab, file="nomoRisk.txt", sep="\t", col.names=F, quote=F)


pdf(file="calibration.pdf", width=5, height=5)
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=1)
cal <- calibrate(f, cmethod="KM", method="boot", u=1, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1),
     xlab="Nomogram-predicted OS (%)", ylab="Observed OS (%)", lwd=1.5, col="green", sub=F)
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=3)
cal <- calibrate(f, cmethod="KM", method="boot", u=3, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", lwd=1.5, col="blue", sub=F, add=T)
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=5)
cal <- calibrate(f, cmethod="KM", method="boot", u=5, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="",  lwd=1.5, col="red", sub=F, add=T)
legend('bottomright', c('1-year', '3-year', '5-year'),
       col=c("green","blue","red"), lwd=1.5, bty = 'n')
dev.off()


library(limma)
library(ggplot2)
library(ggpubr)
library(ggExtra)

gene="xxx"             
corFilter=0.6           
pFilter=0.001            
expFile="xxx"    
setwd("xx")      


rt=read.table(expFile, header=T, sep="\t", check.names=F)
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)
data=data[rowMeans(data)>1,]

group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=log2(data+1)

x=as.numeric(data[gene,])
outTab=data.frame()
for(j in rownames(data)){
  if(gene==j){next}
  y=as.numeric(data[j,])
  corT=cor.test(x, y, method = 'pearson')
  cor=corT$estimate
  pvalue=corT$p.value
  outTab=rbind(outTab, cbind(Query=gene, Gene=j, cor, pvalue))
  if((abs(cor)>corFilter) & (pvalue<pFilter)){
    df1=as.data.frame(cbind(x,y))
    p1=ggplot(df1, aes(x, y)) + 
      xlab(paste0(gene, " expression"))+ ylab(paste0(j, " expression"))+
      geom_point()+ geom_smooth(method="lm", formula=y~x) + theme_bw()+
      stat_cor(method = 'pearson', aes(x =x, y =y))
    pdf(file=paste0("cor.", j, ".pdf"), width=5, height=4.6)
    print(p1)
    dev.off()
  }
}

write.table(file="corResult.txt", outTab, sep="\t", quote=F, row.names=F)
outTab=outTab[abs(as.numeric(outTab$cor))>corFilter & as.numeric(outTab$pvalue)<pFilter,]
write.table(file="corSig.txt", outTab, sep="\t", quote=F, row.names=F)



library(limma)
library(ggplot2)
library(ggpubr)


gene <- "xxx"                   
expFile <- "TCGA_tpm_mRNA.txt"   

plotGenes <- c(
  "TMEM38A", "ATP6V1G3", "BSND", "FOXI1", "HEPACAM2",
  "MYO9B", "STAC3", "CARD9", "FMNL1", "TRPM2", "GMIP"
)

rt <- read.table(expFile, header = TRUE, sep = "\t", check.names = FALSE)
rt <- as.matrix(rt)
rownames(rt) <- rt[, 1]
exp <- rt[, 2:ncol(rt)]

dimnames_list <- list(rownames(exp), colnames(exp))
data <- matrix(as.numeric(as.matrix(exp)), nrow = nrow(exp), dimnames = dimnames_list)
data <- avereps(data)
data <- data[rowMeans(data) > 1, ]
group <- sapply(strsplit(colnames(data), "\\-"), "[", 4)
group <- sapply(strsplit(group, ""), "[", 1)
group <- gsub("2", "1", group)
data <- data[, group == 0]
data <- log2(data + 1)
allGenes <- c(gene, plotGenes)
missGenes <- setdiff(allGenes, rownames(data))
if(length(missGenes) > 0){
  cat("以下基因在表达矩阵中未找到：\n")
  print(missGenes)
}

plotGenes <- intersect(plotGenes, rownames(data))

if(!(gene %in% rownames(data))){
  stop(paste("目标基因", gene, "在表达矩阵中不存在，无法作图。"))
}

x <- as.numeric(data[gene, ])
outTab <- data.frame()
for(j in plotGenes){
  
  y <- as.numeric(data[j, ])
  
  corT <- cor.test(x, y, method = "pearson")
  cor_value <- as.numeric(corT$estimate)
  pvalue <- corT$p.value
  
  outTab <- rbind(
    outTab,
    data.frame(
      Query = gene,
      Gene = j,
      cor = cor_value,
      pvalue = pvalue
    )
  )
  df1 <- data.frame(
    x = x,
    y = y
  )
  p1 <- ggplot(df1, aes(x = x, y = y)) +
    geom_point(color = "#4C72B0", size = 1.8, alpha = 0.75) +
    geom_smooth(method = "lm", formula = y ~ x, color = "#D55E00", fill = "#F4A582") +
    xlab(paste0(gene, " expression")) +
    ylab(paste0(j, " expression")) +
    ggtitle(paste0(gene, " vs ", j)) +
    theme_bw(base_size = 13) +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold"),
      axis.title = element_text(face = "bold"),
      panel.grid = element_blank()
    ) +
    stat_cor(
      method = "pearson",
      aes(label = paste(..r.label.., ..p.label.., sep = "~`,`~")),
      label.x.npc = "left",
      label.y.npc = "top",
      size = 4
    )
  pdf(file = paste0("cor.", j, ".pdf"), width = 5, height = 4.6)
  print(p1)
  dev.off()
}
write.table(
  outTab,
  file = "corResult_selected11.txt",
  sep = "\t",
  quote = FALSE,
  row.names = FALSE
)

cat("完成：已输出 11 个基因的散点图和相关性结果表。\n")

library(limma)
library(pheatmap)

gene="xxx"             
expFile="xxx"    
logFCfilter=1          
fdrFilter=0.05          
rt=read.table(expFile, header=T, sep="\t", check.names=F)
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=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(data))
rt=t(avereps(data))

low=rt[gene,]<=median(rt[gene,])
high=rt[gene,]>median(rt[gene,])
lowRT=rt[,low]
highRT=rt[,high]
conNum=ncol(lowRT)        
treatNum=ncol(highRT)     
data=cbind(lowRT,highRT)
Type=c(rep(1,conNum), rep(2,treatNum))

outTab=data.frame()
for(i in row.names(data)){
  rt=data.frame(expression=data[i,], Type=Type)
  wilcoxTest=wilcox.test(expression ~ Type, data=rt)
  conGeneMeans=mean(data[i,1:conNum])
  treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
  logFC=log2(treatGeneMeans)-log2(conGeneMeans)
  pvalue=wilcoxTest$p.value
  conMed=median(data[i,1:conNum])
  treatMed=median(data[i,(conNum+1):ncol(data)])
  diffMed=treatMed-conMed
  if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){  
    outTab=rbind(outTab,cbind(gene=i,lowMean=conGeneMeans,highMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
  }
}
pValue=outTab[,"pValue"]
fdr=p.adjust(as.numeric(as.vector(pValue)), method="fdr")
outTab=cbind(outTab, fdr=fdr)

outDiff=outTab[( abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$fdr))<fdrFilter),]
write.table(outDiff, file="diff.txt", sep="\t", row.names=F, quote=F)

geneNum=50   
outDiff=outDiff[order(as.numeric(as.vector(outDiff$logFC))),]
diffGeneName=as.vector(outDiff[,1])
diffLength=length(diffGeneName)
hmGene=c()
if(diffLength>(2*geneNum)){
  hmGene=diffGeneName[c(1:geneNum,(diffLength-geneNum+1):diffLength)]
}else{
  hmGene=diffGeneName
}
hmExp=log2(data[hmGene,]+0.01)
Type=c(rep("Low",conNum),rep("High",treatNum))
Type=factor(Type, levels=c("Low","High"))
names(Type)=colnames(data)
Type=as.data.frame(Type)
pdf(file="heatmap.pdf", width=10, height=7)
pheatmap(hmExp, 
         annotation=Type, 
         color = colorRampPalette(c(rep("#88CEE6",5), "white", rep("#E89DA0",5)))(50),
         cluster_cols =F,
         show_colnames = F,
         scale="row",
         fontsize = 8,
         fontsize_row=5,
         fontsize_col=8)
dev.off()
library(ggplot2)
library(ggrepel)
volcanoData <- outTab[outTab$gene != gene, ]

volcanoData$logFC <- as.numeric(as.character(volcanoData$logFC))
volcanoData$fdr <- as.numeric(as.character(volcanoData$fdr))
volcanoData$group <- "Not Sig"
volcanoData$group[volcanoData$logFC > logFCfilter & volcanoData$fdr < fdrFilter] <- "Up"
volcanoData$group[volcanoData$logFC < -logFCfilter & volcanoData$fdr < fdrFilter] <- "Down"
volcanoData$negLogFDR <- -log10(volcanoData$fdr + 1e-300)

label_up <- subset(volcanoData, group == "Up")
label_up <- label_up[order(label_up$fdr, -label_up$logFC), ]
label_up <- head(label_up, 4)

label_down <- subset(volcanoData, group == "Down")
label_down <- label_down[order(label_down$fdr, label_down$logFC), ]
label_down <- head(label_down, 3)

labelData <- rbind(label_up, label_down)

up_num <- sum(volcanoData$group == "Up")
down_num <- sum(volcanoData$group == "Down")

x_lim <- quantile(abs(volcanoData$logFC), 0.99, na.rm = TRUE)
y_max <- max(volcanoData$negLogFDR, na.rm = TRUE)

p_volcano <- ggplot() +
  geom_point(
    data = subset(volcanoData, group == "Not Sig"),
    aes(x = logFC, y = negLogFDR),
    color = "grey88",
    size = 0.65,
    alpha = 0.22
  ) +
  geom_point(
    data = subset(volcanoData, group == "Down"),
    aes(x = logFC, y = negLogFDR),
    color = "#88CEE6",
    size = 1.6,
    alpha = 0.9
  ) +
  geom_point(
    data = subset(volcanoData, group == "Up"),
    aes(x = logFC, y = negLogFDR),
    color = "#E89DA0",
    size = 1.6,
    alpha = 0.9
  ) +
  geom_vline(
    xintercept = c(-logFCfilter, logFCfilter),
    linetype = "dashed",
    linewidth = 0.5,
    color = "grey68"
  ) +
  geom_hline(
    yintercept = -log10(fdrFilter),
    linetype = "dashed",
    linewidth = 0.5,
    color = "grey68"
  ) +
  geom_text_repel(
    data = labelData,
    aes(x = logFC, y = negLogFDR, label = gene),
    size = 3.1,
    box.padding = 0.3,
    point.padding = 0.2,
    segment.color = "grey60",
    segment.size = 0.35,
    max.overlaps = 100
  ) +
  annotate(
    "text",
    x = -x_lim * 0.78,
    y = y_max * 0.98,
    label = paste0("Down: ", down_num),
    color = "#88CEE6",
    size = 4.4,
    fontface = "bold"
  ) +
  annotate(
    "text",
    x = x_lim * 0.78,
    y = y_max * 0.98,
    label = paste0("Up: ", up_num),
    color = "#E89DA0",
    size = 4.4,
    fontface = "bold"
  ) +
  coord_cartesian(xlim = c(-x_lim, x_lim)) +
  labs(
    title = paste0(gene, " high vs low expression"),
    x = expression(log[2]("Fold Change")),
    y = expression(-log[10]("FDR"))
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
    axis.title = element_text(face = "bold", color = "black", size = 14),
    axis.text = element_text(color = "black", size = 12),
    legend.position = "none",
    axis.line = element_line(linewidth = 1),
    axis.ticks = element_line(linewidth = 1),
    plot.margin = margin(10, 12, 10, 10)
  )

pdf("volcano_final.pdf", width = 7, height = 6)
print(p_volcano)
dev.off()

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
library(circlize)
library(RColorBrewer)
library(dplyr)
library(ComplexHeatmap)

pvalueFilter=0.05      
qvalueFilter=0.05       

colorSel="qvalue"
if(qvalueFilter>0.05){
  colorSel="pvalue"
}
ontology.col=c("#00AFBB", "#E7B800", "#90EE90")

rt=read.table("diff.txt", header=T, sep="\t", check.names=F)   

genes=unique(as.vector(rt[,1]))
entrezIDs=mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
entrezIDs=as.character(entrezIDs)
gene=entrezIDs[entrezIDs!="NA"]        

kk=enrichGO(gene=gene, OrgDb=org.Hs.eg.db, pvalueCutoff=1, qvalueCutoff=1, ont="all", readable=T)
GO=as.data.frame(kk)
GO=GO[(GO$pvalue<pvalueFilter & GO$qvalue<qvalueFilter),]

write.table(GO, file="GO.txt", sep="\t", quote=F, row.names = F)


showNum=10
if(nrow(GO)<30){
  showNum=nrow(GO)
}


pdf(file="barplot.pdf", width=8, height=10)
bar=barplot(kk, drop=TRUE, showCategory=showNum, label_format=30, split="ONTOLOGY", color=colorSel) + facet_grid(ONTOLOGY~., scale='free')
print(bar)
dev.off()


pdf(file="bubble.pdf", width=8, height=10)
bub=dotplot(kk, showCategory=showNum, orderBy="GeneRatio", label_format=30, split="ONTOLOGY", color=colorSel) + facet_grid(ONTOLOGY~., scale='free')
print(bub)
dev.off()

data=GO[order(GO$p.adjust),]
datasig=data[data$p.adjust<0.05,,drop=F]
BP = datasig[datasig$ONTOLOGY=="BP",,drop=F]
CC = datasig[datasig$ONTOLOGY=="CC",,drop=F]
MF = datasig[datasig$ONTOLOGY=="MF",,drop=F]
BP = head(BP,6)
CC = head(CC,6)
MF = head(MF,6)
data = rbind(BP,CC,MF)
main.col = ontology.col[as.numeric(as.factor(data$ONTOLOGY))]

BgGene = as.numeric(sapply(strsplit(data$BgRatio,"/"),'[',1))
Gene = as.numeric(sapply(strsplit(data$GeneRatio,'/'),'[',1))
ratio = Gene/BgGene
logpvalue = -log(data$pvalue,10)
logpvalue.col = brewer.pal(n = 8, name = "Reds")
f = colorRamp2(breaks = c(0,2,4,6,8,10,15,20), colors = logpvalue.col)
BgGene.col = f(logpvalue)
df = data.frame(GO=data$ID,start=1,end=max(BgGene))
rownames(df) = df$GO
bed2 = data.frame(GO=data$ID,start=1,end=BgGene,BgGene=BgGene,BgGene.col=BgGene.col)
bed3 = data.frame(GO=data$ID,start=1,end=Gene,BgGene=Gene)
bed4 = data.frame(GO=data$ID,start=1,end=max(BgGene),ratio=ratio,col=main.col)
bed4$ratio = bed4$ratio/max(bed4$ratio)*9.5

pdf("GO.circlize.pdf",width=15,height=15)
par(omi=c(0.1,0.1,0.1,1.5))
circos.par(track.margin=c(0.01,0.01))
circos.genomicInitialize(df,plotType="none")
circos.trackPlotRegion(ylim = c(0, 1), panel.fun = function(x, y) {
  sector.index = get.cell.meta.data("sector.index")
  xlim = get.cell.meta.data("xlim")
  ylim = get.cell.meta.data("ylim")
  circos.text(mean(xlim), mean(ylim), sector.index, cex = 0.8, facing = "bending.inside", niceFacing = TRUE)
}, track.height = 0.08, bg.border = NA,bg.col = main.col)

for(si in get.all.sector.index()) {
  circos.axis(h = "top", labels.cex = 0.6, sector.index = si,track.index = 1,
              major.at=seq(0,max(BgGene),by=100),labels.facing = "clockwise")
}
f = colorRamp2(breaks = c(-1, 0, 1), colors = c("green", "black", "red"))
circos.genomicTrack(bed2, ylim = c(0, 1),track.height = 0.1,bg.border="white",
                    panel.fun = function(region, value, ...) {
                      i = getI(...)
                      circos.genomicRect(region, value, ytop = 0, ybottom = 1, col = value[,2], 
                                         border = NA, ...)
                      circos.genomicText(region, value, y = 0.4, labels = value[,1], adj=0,cex=0.8,...)
                    })
circos.genomicTrack(bed3, ylim = c(0, 1),track.height = 0.1,bg.border="white",
                    panel.fun = function(region, value, ...) {
                      i = getI(...)
                      circos.genomicRect(region, value, ytop = 0, ybottom = 1, col = '#BA55D3', 
                                         border = NA, ...)
                      circos.genomicText(region, value, y = 0.4, labels = value[,1], cex=0.9,adj=0,...)
                    })
circos.genomicTrack(bed4, ylim = c(0, 10),track.height = 0.35,bg.border="white",bg.col="grey90",
                    panel.fun = function(region, value, ...) {
                      cell.xlim = get.cell.meta.data("cell.xlim")
                      cell.ylim = get.cell.meta.data("cell.ylim")
                      for(j in 1:9) {
                        y = cell.ylim[1] + (cell.ylim[2]-cell.ylim[1])/10*j
                        circos.lines(cell.xlim, c(y, y), col = "#FFFFFF", lwd = 0.3)
                      }
                      circos.genomicRect(region, value, ytop = 0, ybottom = value[,1], col = value[,2], 
                                         border = NA, ...)
                      #circos.genomicText(region, value, y = 0.3, labels = value[,1], ...)
                    })
circos.clear()

middle.legend = Legend(
  labels = c('Number of Genes','Number of Select','Rich Factor(0-1)'),
  type="points",pch=c(15,15,17),legend_gp = gpar(col=c('pink','#BA55D3',ontology.col[1])),
  title="",nrow=3,size= unit(3, "mm")
)
circle_size = unit(1, "snpc")
draw(middle.legend,x=circle_size*0.42)

main.legend = Legend(
  labels = c("Biological Process", "Molecular Function","Cellular Component"),  type="points",pch=15,
  legend_gp = gpar(col=ontology.col), title_position = "topcenter",
  title = "ONTOLOGY", nrow = 3,size = unit(3, "mm"),grid_height = unit(5, "mm"),
  grid_width = unit(5, "mm")
)

logp.legend = Legend(
  labels=c('(0,2]','(2,4]','(4,6]','(6,8]','(8,10]','(10,15]','(15,20]','>=20'),
  type="points",pch=16,legend_gp=gpar(col=logpvalue.col),title="-log10(Pvalue)",
  title_position = "topcenter",grid_height = unit(5, "mm"),grid_width = unit(5, "mm"),
  size = unit(3, "mm")
)
lgd = packLegend(main.legend,logp.legend)
circle_size = unit(1, "snpc")
print(circle_size)
draw(lgd, x = circle_size*0.85, y=circle_size*0.55,just = "left")
dev.off()
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
pvalueFilter=0.05      
qvalueFilter=0.05      


colorSel="qvalue"
if(qvalueFilter>0.05){
  colorSel="pvalue"
}
rt=read.table("diff.txt", header=T, sep="\t", check.names=F)     
genes=unique(as.vector(rt[,1]))
entrezIDs=mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
entrezIDs=as.character(entrezIDs)
rt=data.frame(genes, entrezID=entrezIDs)
gene=entrezIDs[entrezIDs!="NA"]        
kk <- enrichKEGG(gene=gene, organism="hsa", pvalueCutoff=1, qvalueCutoff=1)
KEGG=as.data.frame(kk)
KEGG$geneID=as.character(sapply(KEGG$geneID,function(x)paste(rt$genes[match(strsplit(x,"/")[[1]],as.character(rt$entrezID))],collapse="/")))
KEGG=KEGG[(KEGG$pvalue<pvalueFilter & KEGG$qvalue<qvalueFilter),]
write.table(KEGG, file="KEGG.txt", sep="\t", quote=F, row.names = F)
showNum=30
if(nrow(KEGG)<showNum){
  showNum=nrow(KEGG)
}
pdf(file="barplot.pdf", width=9, height=7)
barplot(kk, drop=TRUE, showCategory=showNum, label_format=30, color=colorSel)
dev.off()
pdf(file="bubble.pdf", width = 9, height = 7)
dotplot(kk, showCategory=showNum, orderBy="GeneRatio", label_format=30, color=colorSel)
dev.off()
library(limma)
library(org.Hs.eg.db)
library(clusterProfiler)
library(enrichplot)
gene="xxx"              
expFile="xxx"     
gmtFile="c2.cp.kegg.v7.4.symbols.gmt"    
rt=read.table(expFile, header=T, sep="\t", check.names=F)
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)
data=data[rowMeans(data)>0,]
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(data))
data=t(avereps(data))
dataL=data[,data[gene,]<=median(data[gene,]),drop=F]
dataH=data[,data[gene,]>median(data[gene,]),drop=F]
meanL=rowMeans(dataL)
meanH=rowMeans(dataH)
meanL[meanL<0.00001]=0.00001
meanH[meanH<0.00001]=0.00001
logFC=log2(meanH)-log2(meanL)
logFC=sort(logFC,decreasing=T)
genes=names(logFC)
gmt=read.gmt(gmtFile)
kk=GSEA(logFC, TERM2GENE=gmt, pvalueCutoff = 1)
kkTab=as.data.frame(kk)
kkTab=kkTab[kkTab$pvalue<0.05,]
write.table(kkTab,file="GSEA.result.txt",sep="\t",quote=F,row.names = F)
termNum=5     
if(nrow(kkTab)>=termNum){
  showTerm=row.names(kkTab)[1:termNum]
  gseaplot=gseaplot2(kk, showTerm, base_size=8, title="")
  pdf(file="GSEA.pdf", width=7.5, height=5.6)
  print(gseaplot)
  dev.off()
}
library(limma)
library(estimate)
expFile="xx"      
rt=read.table(expFile, header=T, sep="\t", check.names=F)
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=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
out=rbind(ID=colnames(data),data)
write.table(out,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")
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="TMEscores.txt", sep="\t", quote=F, col.names=F)
library(limma)
library(reshape2)
library(ggpubr)
expFile="xx"          
scoreFile="xx"     
setwd("xxx")  
exp=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)
gene=colnames(exp)[1]
exp=exp[exp$Type=="Tumor",1,drop=F]
exp$Type=ifelse(exp[,gene]>median(exp[,gene]), "High", "Low")
score=read.table(scoreFile, header=T, sep="\t", check.names=F, row.names=1)
score=score[,1:3]
sameSample=intersect(row.names(exp), row.names(score))
exp=exp[sameSample,"Type",drop=F]
score=score[sameSample,,drop=F]
rt=cbind(score, exp)
rt$Type=factor(rt$Type, levels=c("Low", "High"))
data=melt(rt, id.vars=c("Type"))
colnames(data)=c("Type", "scoreType", "Score")
p=ggviolin(data, x="scoreType", y="Score", fill = "Type",
           xlab="",
           ylab="TME score",
           legend.title=gene,
           add = "boxplot", add.params = list(color="white"),
           palette = c("#88CEE6","#E89DA0"), width=1)
p=p+rotate_x_text(45)
p1=p+stat_compare_means(aes(group=Type),
                        method="wilcox.test",
                        symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", " ")),
                        label = "p.signif")
pdf(file="vioplot.pdf", width=6, height=5)
print(p1)
dev.off()
#' CIBERSORT R script v1.03
#' Note: Signature matrix construction is not currently available; use java version for full functionality.
#' Author: Aaron M. Newman, Stanford University (amnewman@stanford.edu)
#' Requirements:
#'       R v3.0 or later. (dependencies below might not work properly with earlier versions)
#'       install.packages('e1071')
#'       install.pacakges('parallel')
#'       install.packages('preprocessCore')
#'       if preprocessCore is not available in the repositories you have selected, run the following:
#'           source("http://bioconductor.org/biocLite.R")
#'           biocLite("preprocessCore")
#' Windows users using the R GUI may need to Run as Administrator to install or update packages.
#' This script uses 3 parallel processes.  Since Windows does not support forking, this script will run
#' single-threaded in Windows.
#'
#' Usage:
#'       Navigate to directory containing R script
#'
#'   In R:
#'       source('CIBERSORT.R')
#'       results <- CIBERSORT('sig_matrix_file.txt','mixture_file.txt', perm, QN)
#'
#'       Options:
#'       i)  perm = No. permutations; set to >=100 to calculate p-values (default = 0)
#'       ii) QN = Quantile normalization of input mixture (default = TRUE)
#'
#' Input: signature matrix and mixture file, formatted as specified at http://cibersort.stanford.edu/tutorial.php
#' Output: matrix object containing all results and tabular data written to disk 'CIBERSORT-Results.txt'
#' License: http://cibersort.stanford.edu/CIBERSORT_License.txt
#' Core algorithm
#' @param X cell-specific gene expression
#' @param y mixed expression per sample
#' @export
CoreAlg <- function(X, y){
  
  #try different values of nu
  svn_itor <- 3
  
  res <- function(i){
    if(i==1){nus <- 0.25}
    if(i==2){nus <- 0.5}
    if(i==3){nus <- 0.75}
    model<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
    model
  }
  
  if(Sys.info()['sysname'] == 'Windows') out <- mclapply(1:svn_itor, res, mc.cores=1) else
    out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)
  
  nusvm <- rep(0,svn_itor)
  corrv <- rep(0,svn_itor)
  
  #do cibersort
  t <- 1
  while(t <= svn_itor) {
    weights = t(out[[t]]$coefs) %*% out[[t]]$SV
    weights[which(weights<0)]<-0
    w<-weights/sum(weights)
    u <- sweep(X,MARGIN=2,w,'*')
    k <- apply(u, 1, sum)
    nusvm[t] <- sqrt((mean((k - y)^2)))
    corrv[t] <- cor(k, y)
    t <- t + 1
  }
  
  #pick best model
  rmses <- nusvm
  mn <- which.min(rmses)
  model <- out[[mn]]
  
  #get and normalize coefficients
  q <- t(model$coefs) %*% model$SV
  q[which(q<0)]<-0
  w <- (q/sum(q))
  
  mix_rmse <- rmses[mn]
  mix_r <- corrv[mn]
  
  newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
  
}

#' do permutations
#' @param perm Number of permutations
#' @param X cell-specific gene expression
#' @param y mixed expression per sample
#' @export
doPerm <- function(perm, X, Y){
  itor <- 1
  Ylist <- as.list(data.matrix(Y))
  dist <- matrix()
  
  while(itor <= perm){
    #print(itor)
    
    #random mixture
    yr <- as.numeric(Ylist[sample(length(Ylist),dim(X)[1])])
    
    #standardize mixture
    yr <- (yr - mean(yr)) / sd(yr)
    
    #run CIBERSORT core algorithm
    result <- CoreAlg(X, yr)
    
    mix_r <- result$mix_r
    
    #store correlation
    if(itor == 1) {dist <- mix_r}
    else {dist <- rbind(dist, mix_r)}
    
    itor <- itor + 1
  }
  newList <- list("dist" = dist)
}

#' Main functions
#' @param sig_matrix file path to gene expression from isolated cells
#' @param mixture_file heterogenous mixed expression
#' @param perm Number of permutations
#' @param QN Perform quantile normalization or not (TRUE/FALSE)
#' @export
CIBERSORT <- function(sig_matrix, mixture_file, perm=0, QN=TRUE){
  library(e1071)
  library(parallel)
  library(preprocessCore)
  
  #read in data
  X <- read.table(sig_matrix,header=T,sep="\t",row.names=1,check.names=F)
  Y <- read.table(mixture_file, header=T, sep="\t", row.names=1,check.names=F)
  
  X <- data.matrix(X)
  Y <- data.matrix(Y)
  
  #order
  X <- X[order(rownames(X)),]
  Y <- Y[order(rownames(Y)),]
  
  P <- perm #number of permutations
  
  #anti-log if max < 50 in mixture file
  if(max(Y) < 50) {Y <- 2^Y}
  
  #quantile normalization of mixture file
  if(QN == TRUE){
    tmpc <- colnames(Y)
    tmpr <- rownames(Y)
    Y <- normalize.quantiles(Y)
    colnames(Y) <- tmpc
    rownames(Y) <- tmpr
  }
  
  #intersect genes
  Xgns <- row.names(X)
  Ygns <- row.names(Y)
  YintX <- Ygns %in% Xgns
  Y <- Y[YintX,]
  XintY <- Xgns %in% row.names(Y)
  X <- X[XintY,]
  #if(substr(Sys.Date(),1,4)>(2000+22)){next}
  #standardize sig matrix
  X <- (X - mean(X)) / sd(as.vector(X))
  
  #empirical null distribution of correlation coefficients
  if(P > 0) {nulldist <- sort(doPerm(P, X, Y)$dist)}
  
  #print(nulldist)
  
  header <- c('Mixture',colnames(X),"P-value","Correlation","RMSE")
  #print(header)
  
  output <- matrix()
  itor <- 1
  mixtures <- dim(Y)[2]
  pval <- 9999
  
  #iterate through mixtures
  while(itor <= mixtures){
    
    y <- Y[,itor]
    
    #standardize mixture
    y <- (y - mean(y)) / sd(y)
    
    #run SVR core algorithm
    result <- CoreAlg(X, y)
    
    #get results
    w <- result$w
    mix_r <- result$mix_r
    mix_rmse <- result$mix_rmse
    
    #calculate p-value
    if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))}
    
    #print output
    out <- c(colnames(Y)[itor],w,pval,mix_r,mix_rmse)
    if(itor == 1) {output <- out}
    else {output <- rbind(output, out)}
    
    itor <- itor + 1
    
  }
  
  #save results
  write.table(rbind(header,output), file="CIBERSORT-Results.txt", sep="\t", row.names=F, col.names=F, quote=F)
  
  #return matrix object containing all results
  obj <- rbind(header,output)
  obj <- obj[,-1]
  obj <- obj[-1,]
  obj <- matrix(as.numeric(unlist(obj)),nrow=nrow(obj))
  rownames(obj) <- colnames(Y)
  colnames(obj) <- c(colnames(X),"P-value","Correlation","RMSE")
  obj
}


library("limma")     
expFile="xx"     
rt=read.table(expFile, header=T, sep="\t", check.names=F)
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)
data=data[rowMeans(data)>0,]
v=voom(data, plot=F, save.plot=F)
out=v$E
out=rbind(ID=colnames(out), out)
write.table(out,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)      
source("Gene25.CIBERSORT.R")
results=CIBERSORT("ref.txt", "uniq.symbol.txt", perm=1000, QN=TRUE)

library(limma)
library(reshape2)
library(ggpubr)
library(vioplot)
library(ggExtra)

expFile="xxxx"          
immFile="CIBERSORT-Results.txt"  
pFilter=0.05                     
fdrFilter=0.05                 

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

tumorData=rt[rt$Type=="Tumor",1,drop=F]
tumorData=as.matrix(tumorData)
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData))
data=avereps(tumorData)

data=as.data.frame(data)
data$gene=ifelse(data[,gene]>median(data[,gene]), "High", "Low")

immune=read.table(immFile, header=T, sep="\t", check.names=F, row.names=1)
immune=immune[immune[,"P-value"]<pFilter,]
immune=as.matrix(immune[,1:(ncol(immune)-3)])


group=sapply(strsplit(row.names(immune),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
immune=immune[group==0,]
row.names(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", row.names(immune))
immune=avereps(immune)

sameSample=intersect(row.names(immune), row.names(data))
rt=cbind(immune[sameSample,,drop=F], data[sameSample,,drop=F])

data=rt[,-(ncol(rt)-1)]
data=melt(data,id.vars=c("gene"))
colnames(data)=c("gene", "Immune", "Expression")


group=levels(factor(data$gene))
data$gene=factor(data$gene, levels=c("Low","High"))
data$Immune=factor(data$Immune, levels=unique(data$Immune))

diffTab=data.frame()

for(i in levels(data$Immune)){
  tmp=data[data$Immune==i,]
  test=wilcox.test(Expression~gene, data=tmp)
  
  diffTab=rbind(
    diffTab,
    data.frame(
      Immune=i,
      pvalue=test$p.value
    )
  )
}

diffTab$FDR=p.adjust(diffTab$pvalue, method="BH")

diffTab$signif=ifelse(
  diffTab$FDR<0.001, "***",
  ifelse(diffTab$FDR<0.01, "**",
         ifelse(diffTab$FDR<0.05, "*", ""))
)

write.table(diffTab,
            file="immune.diff.FDR.result.txt",
            sep="\t",
            row.names=F,
            quote=F)


bioCol=c("#88CEE6","#E89DA0","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length(group)]


boxplot=ggboxplot(data, x="Immune", y="Expression", fill="gene",
                  xlab="",
                  ylab="Fraction",
                  legend.title=gene,
                  width=0.8,
                  palette=bioCol)+
  rotate_x_text(50)+
  scale_y_continuous(expand=expansion(mult=c(0.05, 0.20)))

sigDiff=diffTab[diffTab$FDR<fdrFilter,]

if(nrow(sigDiff)>0){
  
  yMax=max(data$Expression, na.rm=T)
  yMin=min(data$Expression, na.rm=T)
  yRange=yMax-yMin
  
  lineY=yMax+0.08*yRange
  textY=yMax+0.11*yRange
  
  for(i in 1:nrow(sigDiff)){
    
    cell=sigDiff$Immune[i]
    xpos=which(levels(data$Immune)==cell)
    
    boxplot=boxplot+
      annotate("segment",
               x=xpos-0.25,
               xend=xpos+0.25,
               y=lineY,
               yend=lineY,
               linewidth=0.5)
    
    boxplot=boxplot+
      annotate("segment",
               x=xpos-0.25,
               xend=xpos-0.25,
               y=lineY-0.01*yRange,
               yend=lineY,
               linewidth=0.5)
    
    boxplot=boxplot+
      annotate("segment",
               x=xpos+0.25,
               xend=xpos+0.25,
               y=lineY-0.01*yRange,
               yend=lineY,
               linewidth=0.5)
    
    boxplot=boxplot+
      annotate("text",
               x=xpos,
               y=textY,
               label=sigDiff$signif[i],
               size=6,
               fontface="bold")
  }
}

pdf(file="immune.diff.pdf", width=7, height=6)
print(boxplot)
dev.off()


inputFile="cor.result.txt"      
data = read.table(inputFile, header=T, sep="\t", check.names=F)

if(!("FDR" %in% colnames(data))){
  data$FDR = p.adjust(data$pvalue, method="BH")
}

p.col = c('gold','pink','orange','LimeGreen','darkgreen')
fcolor = function(x,p.col){
  color = ifelse(x>0.8,p.col[1],ifelse(x>0.6,p.col[2],ifelse(x>0.4,p.col[3],
                                                             ifelse(x>0.2,p.col[4], p.col[5])
  )))
  return(color)
}

p.cex = seq(2.5, 5.5, length=5)
fcex = function(x){
  x=abs(x)
  cex = ifelse(x<0.1,p.cex[1],ifelse(x<0.2,p.cex[2],ifelse(x<0.3,p.cex[3],
                                                           ifelse(x<0.4,p.cex[4],p.cex[5]))))
  return(cex)
}

points.color = fcolor(x=data$FDR,p.col=p.col)
data$points.color = points.color

points.cex = fcex(x=data$cor)
data$points.cex = points.cex
data=data[order(data$cor),]

xlim = ceiling(max(abs(data$cor))*10)/10        

pdf(file="Lollipop.FDR.pdf", width=9, height=7)    
layout(mat=matrix(c(1,1,1,1,1,0,2,0,3,0),nc=2),width=c(8,2.2),heights=c(1,2,1,2,1))
par(bg="white",las=1,mar=c(5,18,2,4),cex.axis=1.5,cex.lab=2)

plot(1,type="n",xlim=c(-xlim,xlim),ylim=c(0.5,nrow(data)+0.5),
     xlab="Correlation Coefficient",ylab="",yaxt="n",yaxs="i",axes=F)

rect(par('usr')[1],par('usr')[3],par('usr')[2],par('usr')[4],
     col="#F5F5F5",border="#F5F5F5")

grid(ny=nrow(data),col="white",lty=1,lwd=2)
segments(x0=data$cor,y0=1:nrow(data),x1=0,y1=1:nrow(data),lwd=4)
points(x=data$cor,y = 1:nrow(data),col = data$points.color,pch=16,cex=data$points.cex)
text(par('usr')[1],1:nrow(data),data$Cell,adj=1,xpd=T,cex=1.5)
fdr.text=ifelse(data$FDR<0.001,'<0.001',sprintf("%.03f",data$FDR))
redcutoff_cor=0
redcutoff_FDR=0.05
text(par('usr')[2],1:nrow(data),fdr.text,adj=0,xpd=T,
     col=ifelse(abs(data$cor)>redcutoff_cor & data$FDR<redcutoff_FDR,"red","black"),
     cex=1.5)

axis(1,tick=F)
par(mar=c(0,4,3,4))
plot(1,type="n",axes=F,xlab="",ylab="")
legend("left",legend=c(0.1,0.2,0.3,0.4,0.5),col="black",
       pt.cex=p.cex,pch=16,bty="n",cex=2,title="abs(cor)")

par(mar=c(0,6,4,6),cex.axis=1.5,cex.main=2)
barplot(rep(1,5),horiz=T,space=0,border=NA,col=p.col,
        xaxt="n",yaxt="n",xlab="",ylab="",main="FDR")

axis(4,at=0:5,c(1,0.8,0.6,0.4,0.2,0),tick=F)

dev.off()
library(limma)
library(reshape2)
library(ggplot2)
library(ggpubr)
library(corrplot)

fdrFilter=0.05                
corCutoff=0.3                
geneName="ARHGAP22"            
expFile="xxx"   
geneFile="gene.txt"           
groupCol <- c("Low" = "#88CEE6", "High" = "#E89DA0")

corCol <- colorRampPalette(c("#88CEE6", "white", "#E89DA0"))(100)
rt <- read.table(expFile, header = TRUE, sep = "\t", check.names = FALSE)
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)
gene <- read.table(geneFile, header = FALSE, sep = "\t", check.names = FALSE)
sameGene <- intersect(rownames(data), as.vector(gene[, 1]))

if (!(geneName %in% rownames(data))) {
  stop(paste("目标基因不在表达矩阵中:", geneName))
}

data <- t(data[c(geneName, sameGene), ])
data <- log2(data + 1)
group <- sapply(strsplit(rownames(data), "\\-"), "[", 4)
group <- sapply(strsplit(group, ""), "[", 1)
group <- gsub("2", "1", group)
data <- data[group == 0, ]
rownames(data) <- gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1-\\2-\\3", rownames(data))
data <- t(avereps(data))

x <- as.numeric(data[geneName, ])
allTab <- data.frame()

for (i in sameGene) {
  if (i == geneName) next
  
  y <- as.numeric(data[i, ])
  corT <- cor.test(x, y, method = "pearson")
  corValue <- as.numeric(corT$estimate)
  pvalue <- corT$p.value
  
  if (!is.na(corValue)) {
    allTab <- rbind(allTab,
                    data.frame(Query = geneName,
                               Gene = i,
                               cor = corValue,
                               pvalue = pvalue))
  }
}

if (nrow(allTab) > 0) {
  allTab$FDR <- p.adjust(allTab$pvalue, method = "BH")
  
  allTab <- allTab[order(abs(allTab$cor), decreasing = TRUE), ]
}

write.table(allTab, file = "corResult.all.txt", sep = "\t", quote = FALSE, row.names = FALSE)
outTab <- allTab[allTab$FDR < fdrFilter & abs(allTab$cor) >= corCutoff, ]

write.table(outTab, file = "corResult.txt", sep = "\t", quote = FALSE, row.names = FALSE)
if (nrow(outTab) > 1) {
  plotData <- t(data[c(geneName, as.vector(outTab[, "Gene"])), ])
  M <- cor(plotData, use = "pairwise.complete.obs")
  
  pdf("corplot_circle.pdf", width = 7, height = 7)
  corrplot(M,
           method = "circle",
           order = "original",
           type = "upper",
           col = corCol,
           tl.col = "black")
  dev.off()
  
  pdf("corplot_number.pdf", width = 8, height = 8)
  corrplot(M,
           order = "original",
           method = "color",
           number.cex = 0.7,
           addCoef.col = "black",
           diag = TRUE,
           tl.col = "black",
           col = corCol)
  dev.off()
}