
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.12")


geneName="FTL"                                                           #
 
  #【A】#############################
  ###  ---different expression  ---### 
  ####################################---#-----------------------------------------------------------------
  
  
  library(limma)
  
  filename.singlegeneName.txt <- paste("The data for Figure S1A and B.txt",sep = "")
  rt=read.table(filename.singlegeneName.txt,sep="\t",header=T,row.names=1,check.names=F)#row.names=F,
  #duplicated(rt$ID)
  colnames(rt)=c("expression","Type")
  rt$expression=log2(rt$expression+1)
  
  ylab = paste(geneName," log2 (expression)",sep="")
  
  
##ggplot------
  
  library(plyr)
  library(ggplot2)
  library(ggpubr)
  
  options(stringsAsFactors=F)
  
  
 #p-value
  clinical= colnames(rt)[2]                                                      
  gene=colnames(rt)[1]                                                           
  
  xlabel=vector()
  tab1=table(rt[,clinical])
  labelNum=length(tab1)
  dotCol=MyPal1
  if(labelNum==3){
    dotCol=c(2,3,4)
  }
  if(labelNum==4){
    dotCol=c(2,3,4,5)
  }
  if(labelNum>4){
    dotCol=rainbow(labelNum)
  }
  for(i in 1:labelNum){
    xlabel=c(xlabel,names(tab1[i]) )
  }
  
  outTab=data.frame()
  
  rt1=rbind(expression=rt[,gene],clinical=rt[,clinical])
  rt1=data.frame(t(rt1))
  rt1[,1]=as.numeric(as.character(rt1[,1]))
  
  class(rt1[,1])
  if(labelNum==2){
    wilcoxTest<-wilcox.test(expression ~ clinical, data=rt1)
  }else{
    wilcoxTest<-kruskal.test(expression ~ clinical, data = rt1)}
  pValue=wilcoxTest$p.value
  outTab=rbind(outTab,cbind(gene=gene,pVal=pValue))
  pval=0
  if(pValue<0.001){
    pval=signif(pValue,4)
    pval=format(pval, scientific = TRUE)
  }else{
    pval=round(pValue,3)
  }
  
  plotType="dotbPlot"
  
  clinical= colnames(rt)[2]                                                     
  gene=colnames(rt)[1]                                                        
  
  data=rt[c(gene,clinical)]
  data=data[(data[,gene] > -4),]
  
  colnames(data)=c("gene","clinical")
  data=data[(data[,"clinical"]!="unknow"),]
  data=na.omit(data)
  
  #group
  group=levels(factor(data$clinical))
  data$clinical=factor(data$clinical, levels=group)
  comp=combn(group,2)
  my_comparisons=list()
  
  for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
  
  #boxplot
  plotcli<-ggboxplot(data, x="clinical", y="gene", color="clinical",
                     xlab=clinical,
                     ylab=ylab,
                     legend.title=clinical,
                     add = "jitter")+ 
    stat_compare_means(comparisons = my_comparisons)
  
  width=5
  height=6
  outFile=paste("Figure 2A.pdf",sep="")
  pdf(file=outFile,width=width,height=height)
  print(plotcli)
  dev.off()
  
  
  #【B】#############################
  ###---    paired expression   ---### 
  ####################################---#-------------------------------------------------------------------
  
  MyPal <- c("#E41A1C","#377EB8","#4DAF4A","#FF7F00","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")#（"#0571b0"-蓝,"#ca0020"-红）
  
  rt.singleGene=read.table(filename.singlegeneName.txt,sep="\t",header=T,check.names=F)             #读取文件
  colnames(rt.singleGene)=c("ID","expression","Type")
  rt.singleGene$expression=log2(rt.singleGene$expression+1)
  colnames(rt.singleGene)=c("ID",geneName,"Type")
  ylab = paste(geneName," log2 (expression)",sep="")
  
  add.ID=rt.singleGene[,1]
  
  #ID
  rt.singleGene.addID=cbind(rt.singleGene,add.ID)
  rt.singleGene.addID[,1]=substr(rt.singleGene.addID[,1],1,12)
  
  #normal-tumor
  rt.singleGene.addID2 = rt.singleGene.addID
  rt.singleGene.Normal=rt.singleGene.addID[grep(pattern="Normal",rt.singleGene.addID[,3]),]
  rt.singleGene.Tumor=rt.singleGene.addID[grep(pattern="Tumor",rt.singleGene.addID[,3]),]
  
  rt.singleGene.Normal0=rt.singleGene.Normal
  rt.singleGene.Tumor0=rt.singleGene.Tumor
  rownames(rt.singleGene.Normal0)=rt.singleGene.Normal0[,1]
  rownames(rt.singleGene.Tumor0)=rt.singleGene.Tumor0[,1]
  #same name
  rt.singleGene.Normal=rt.singleGene.addID[grep(pattern="Normal",rt.singleGene.addID[,3]),]
  rownames(rt.singleGene.Normal)=rt.singleGene.Normal[,1]
  selename.normal = rownames(rt.singleGene.Normal)
  
  rownames(rt.singleGene.Tumor)=rt.singleGene.Tumor[,1]
  rt.singleGene.Tumor=rt.singleGene.Tumor[selename.normal,]
  rt.singleGene.Tumor=na.omit(rt.singleGene.Tumor)
  
  rt.singleGene.Pairprepare=merge(rt.singleGene.Tumor,rt.singleGene.Tumor,all=FALSE)
  
  selename.Pair= rownames(rt.singleGene.Pairprepare)=rt.singleGene.Pairprepare[,1]
  
  
  rt.singleGene.Normal.Pair=rt.singleGene.Normal0[selename.Pair,]
  rt.singleGene.Tumor.Pair=rt.singleGene.Tumor0[selename.Pair,]
  
  rownames(rt.singleGene.Normal.Pair)<-NULL
  rownames(rt.singleGene.Tumor.Pair)<-NULL
  
  rt.singleGene.NormalTumor.Pair=rbind(rt.singleGene.Normal.Pair,rt.singleGene.Tumor.Pair)
  
  rt.singleGene.NormalTumor.Pair=rt.singleGene.NormalTumor.Pair[,c(4,2,3,1)]
  colnames(rt.singleGene.NormalTumor.Pair)[1] ="ID"
  rt.singleGene.geneExp=rt.singleGene.NormalTumor.Pair[,-4]
  rt.singleGene.Group=rt.singleGene.NormalTumor.Pair[,-2]
  
  #filename.paired.geneName.txt <- paste(geneName,"-","Exp.Paired.singleGene.txt",sep = "")
  #write.table(rt.singleGene.geneExp,file=filename.paired.geneName.txt,sep="\t",quote=F,row.names=F)
  
  #filename.paired.Group.txt <- paste(geneName,"-","Paired.samplegroup.txt",sep = "")
  #write.table(rt.singleGene.Group,file=filename.paired.Group.txt,sep="\t",quote=F,row.names=F,col.names=F)
  
  #---paired 
  group = rt.singleGene.Group#
  #group = read.table(filename.paired.Group.txt,header=F,sep="\t")                               #读取分组信息文件
  df = rt.singleGene.geneExp#
  #df = read.table(filename.paired.geneName.txt,row.names=1,header=T,sep="\t",check.names=F)        #读取输入文件
  row.names(group) = group[,1]
  row.names(df) = df[,1]
  df=df[,-1]
  m = match(group[,1],rownames(df))
  df = df[m,]
  
  
#dot----
  plotType="dot"
  data = df[,1]
  cell = colnames(df)[1]
  cell = gsub(' ','_',cell)
  #Plot = function(data,group,outpdf){
  xfactors = as.factor(group[,2])
  xnumsample = as.numeric(xfactors)
  xaxis = levels(xfactors)
  link = group[,3]
  links = unique(group[,3])
  x1data = data[xnumsample==1]
  x2data = data[xnumsample==2]
  wilcoxP=wilcox.test(x1data,x2data,exact=FALSE)$p.value#)$p.value
  pvalue.Pairdot=signif(wilcoxP,3)
  if(pvalue.Pairdot<0.001){
    pval=signif(pvalue.Pairdot,4)
    pval=format(pvalue.Pairdot, scientific = TRUE)
  }else{
    pval=round(pvalue.Pairdot,3)
  }
  
  width=5
  height=6
  outpdf=paste("Figure 2B.pdf",sep="")
  pdf(outpdf,width=width,height=height)
  par(las=1)
  
  yMin=min(data)-max(data)/10                     #纵轴最小值                      #纵轴最小值
  yMax=max(data)*1.1                     #纵轴最大值
  
  plot(1,xlim=c(0.5,2.5),ylim=c(yMin,yMax),type="n",xlab="",ylab = ylab,xaxt="n",cex.main=1.5, cex.lab=1.3, cex.axis=1.2)
  points(rep(1,length(x1data)),x1data,pch=16,cex=2,col=MyPal[1])
  points(rep(2,length(x2data)),x2data,pch=15,cex=2,col=MyPal[2])
  axis(1,1:2,xaxis)
  for(i in links){
    w = which(link==i)
    x1 = xnumsample[w[1]]
    y1 = data[w[1]]
    x2 = xnumsample[w[2]]
    y2 = data[w[2]]
    segments(x1,y1,x2,y2)
  }
  par(xpd=T)
  arrows(1,max(data)*1.1,2,max(data)*1.1,angle=90,code=3,length=0.1)
  text(1.5,max(data)*1.05,paste("p =",pval),pos=3,cex=1)
  dev.off()
  #}
  
  #data = df[,1]
  #cell = colnames(df)[1]
  #cell = gsub(' ','_',cell)
  #outpdf=paste(geneName,"-diff.Paired.",plotType,"-h",height,"w",width,"-",pvalue.Pairdot,".pdf",sep="")
  #Plot(data,group,outpdf)
 