R/SXTpls.R

SXTpls <- function(sample=NULL,
                   qc=NULL,
                   info=NULL,
                  #used data
                  scalemethod="auto",
                  plsmethod=c("plsreg","plsr"),
                  width=7,
                  height=7,
                  QC=FALSE,
                  text=FALSE,
                  ellipse=FALSE,
                  color=c("palegreen","firebrick1","royalblue","yellow","black","cyan","gray48"),
                  shape=c(17,19,15,18,2,8,11),
                  cexa=1,
                  cexlab = 1.5,
                  cexaxis = 1.8,
                  cexlegend = 1.8)
  #parameter setting
{
   # browser()
  options(warn=-1)
  if (is.null(sample))  stop("sample is NULL")
  if (!is.null(qc))
  {if (ncol(sample)!=ncol(qc)) stop("the column number of sample and qc must be same")}
  if (is.null(qc)&QC) stop("QC shoud be FALSE because qc is NULL")
  if (is.null(info)) stop("info must not be NULL")


  #load needed packages
  need.packages1<-c("plspm","plsdepot","scatterplot3d","pls","lars",
                    "ROCR","ggplot2","pROC","ellipse")

  packages <- library()[[2]][,1]
  for (i in need.packages1) {
    if (!any(packages==i)) {install.packages(i)}
  }


  need.packages2<-c("SXTscale","SXTvip")
  for (i in need.packages2) {
    if (!any(packages==i)) stop(paste("Please install",i,"package form Japse Shen!"))
  }

  require(plspm); require(plsdepot); require(scatterplot3d)
  require(pls); require(lars); require(ROCR); require(ggplot2)
  require(pROC); require(ellipse); require(SXTscale); require(SXTvip)

  int <- sample
  index <- NULL
  for (i in 1:length(info)) {
    index1<-as.character(info[[i]])
    index<-c(index,index1)
  }
  if (length(which(index==""))!=0)  {index<-index[-which(index=="")]}


  index <- index[!is.na(index)]
  index <- match(index, rownames(int))
  index <- index[!is.na(index)]
  int <- int[index, ]

  ifelse(QC,int<-rbind(int,qc),int<-int)
  #######
  q<-grep("QC",rownames(int))
  name<-rownames(int)
  # browser()
  Y <- NULL
  label<-list()
  for (i in 1:length( info )) {
    label[[i]]<-match(as.character(info[[i]]),name)
    label[[i]]<-label[[i]][!is.na(label[[i]])]
    Y[label[[i]]]<-i-1
  }
  if (QC) {Y[q]<- length(info)}

  int.scale <- SXTscale(int,method=scalemethod)
  # int.Y<-SXTscale(Y,method=scalemethod)
  int.Y <- Y
  ncompa <- nrow(int) - 1

  if (plsmethod=="plsr") {
    pls1<-plsr(int.Y~int.scale,scale = FALSE,validation = "CV",ncomp = ncompa,method = "oscorespls")
    save(pls1,file="pls1")
    #########select the number of compents#################

    msep<-MSEP(pls1)
    save(msep,file="msep")
    msep<-msep$val[,,]

    yesnot<-"y"
    while (yesnot=="y"|yesnot=="") {
      comps.number<-readline("How many comps do you want to see? ")
      while (!exists("comps.number")|comps.number=="")
      {cat("You must give a comps number to continute!!!\n")
        comps.number<-readline("How many comps do you want to see? ")}
      comps.number<-as.numeric(comps.number)
      plot(x=c(1:comps.number),y=msep[1,2:(comps.number+1)],type="b",col="firebrick1",pch=20,
           xlab="ncomp",ylab="MSEP",cex.lab=1.3,cex.axis=1.3)
      points(x=c(1:comps.number),y=msep[2,2:(comps.number+1)],type="b",pch=2)
      legend("top",legend = c("CV","adjCV"),col = c("firebrick1","black"),pch=c(20,2),
             bty = "n",cex = 1.3,pt.cex = 1.3)
      yesnot<-readline("Do you want to see the next plot? (y/n)")
    }

    pdf("MSEP plot.pdf")
    par(mar=c(5,5,4,2))
    plot(x=c(1:comps.number),y=msep[1,2:(comps.number+1)],type="b",col="firebrick1",pch=20,
         xlab="ncomp",ylab="MSEP",cex.lab=1.3,cex.axis=1.3)
    points(x=c(1:comps.number),y=msep[2,2:(comps.number+1)],type="b",pch=2)
    legend("top",legend = c("CV","adjCV"),col = c("firebrick1","black"),pch=c(20,2),
           bty = "n",cex = 1.3,pt.cex = 1.3)
    dev.off()

    number<-readline("Please type number and press Enter  to continute:  ")
    while (!exists("number")|number=="") {cat("You must give a number to continute!!!\n")
      number<-readline("Please type comps number value and press Enter  to continute: ")}
    number<-as.numeric(number)

    ##################construct final pls model###################
    pls2 <- plsr(int.Y~int.scale, scale = FALSE,validation = "CV",ncomp = number,method = "oscorespls")
    save(pls2, file = "pls2")
    vip <- SXTvip(pls2)
    save(vip,file = "vip")

    scores <- scores(pls2)
    x <- scores[,1]
    y <- scores[,2]
    if (number > 2) {z <- scores[,3];zmin <- 1.2*min(z);zmax <- 1.2*max(z)}

    xmin <- 1.2 * min(x)
    xmax <- 1.2 * max(x)
    ymin <- 1.2 * min(y)
    ymax <- 1.2 * max(y)
  }

  else {
    # browser()
    require(SXTdummy)
    dummy <- SXTdummy(Y)
    # int.dummy<-SXTscale(dummy,method=scalemethod)
    int.dummy <- dummy
    # ncompa = nrow(int.scale) - 1
    ncompa <- min(nrow(int), ncol(int))
    pls1 <- plsreg1(int.scale,Y, comps = ncompa)
    save(pls1,file = "pls1")
    #########select the number of compents#################
    Q2cum <- pls1$Q2[,5]
    Q2cum[is.nan(Q2cum)] <- 1
    yesnot<-"y"
    while (yesnot=="y"|yesnot=="") {
      comps.number<-readline("How many comps do you want to see? ")
      while (!exists("comps.number")|comps.number=="") {cat("You must give a comps number to continute!!!\n")
        comps.number<-readline("How many comps do you want to see? ")}
      comps.number<-as.numeric(comps.number)
      barplot(Q2cum[1:comps.number],xlab="ncomp",ylab="Q2cum",cex.lab=1.3,cex.axis=1.3)
      a <- barplot(Q2cum[1:comps.number],xlab="ncomp",ylab="Q2cum",cex.lab=1.3,cex.axis=1.3)
      abline(h=0)
      points(a,Q2cum[1:comps.number],type="b",col="red",pch=20,cex=2)
      yesnot <- readline("Do you want to see the next plot? (y/n)")
    }
    pdf("Q2cum plot.pdf",width=7,height=7)
    par(mar=c(5,5,4,2))
    barplot(Q2cum[1:comps.number],xlab="ncomp",ylab="Q2cum",cex.lab=1.3,cex.axis=1.3)
    a<-barplot(Q2cum[1:comps.number],xlab="ncomp",ylab="Q2cum",cex.lab=1.3,cex.axis=1.3)
    abline(h=0)
    points(a,Q2cum[1:comps.number],type="b",col="red",pch=20,cex=2)
    dev.off()

    number<-readline("Please type number and press Enter  to continute:  ")
    while (!exists("number")|number=="") {cat("You must give a number to continute!!!\n")
      number<-readline("Please type comps number value and press Enter  to continute: ")}
    number<-as.numeric(number)

    ##################construct final pls model###################
    cat(paste("Construct PLS model with all peaks using",number,"comps ...","\n"))
    pls2 <- plsreg1(int.scale,Y,comps = number)
    save(pls2,file = "pls2")
    pls.temp <- plsreg2(int.scale,int.dummy, comps = number)
    vip <- pls.temp$VIP
    Q2cum <- pls2$Q2[,5]
    R2cum <- cumsum(pls2$R2)

    write.csv(cbind(R2cum,Q2cum),"R2Q2.csv",row.names = F)

    ##draw barplot of Q2cum, R2Xcum and R2Ycum
    Q2R2 <- cbind( R2cum, Q2cum)
    colnames(Q2R2) <- c( "R2cum","Q2cum")
    pdf("Q2R2cum.pdf",width = 8,height = 6)
    par(mar=c(5,5,4,2))
    barplot( t(Q2R2), beside = T, col = c( "palegreen", "royalblue"),
             cex.lab = 1.3, cex.axis=1.3, cex.names = 1.3)
    abline( h = 0)
    legend( "topleft", legend = c("R2Ycum", "Q2cum"),pch=15,
            col = c( "palegreen", "royalblue"), cex = 1.3, pt.cex = 1.3, bty = "n")
    dev.off()

    save(vip,file="vip")
    save(Q2cum,file = "Q2cum")
    save(R2cum,file = "R2cum")

    x <- pls2$x.scores[,1]
    y <- pls2$x.scores[,2]
    if (number>2) {z<-pls2$x.scores[,3];zmin<-1.2*min(z);zmax<-1.2*max(z)}

    xmin <- 1.2 * min(x)
    xmax <- 1.2 * max(x)
    ymin <- 1.2 * min(y)
    ymax <- 1.2 * max(y)
  }

  legend<-NULL
  for (i in 1:length(label)) {
    legend[label[[i]]]<- names(info)[i]
  }
  if (QC) {legend[q]<-"QC"}

  colour<-NULL
  colourlist<-color
  for (i in 1:length(label)) {
    colour[label[[i]]]<-colourlist[i]
  }
  if (QC) {colour[q]<-colourlist[length(info)+1]}

  pcha<-NULL
  pchalist<-shape
  for (i in 1:length(label)) {
    pcha[label[[i]]]<-pchalist[i]
  }
  if (QC) {pcha[q]<-pchalist[length(info)+1]}

  #PLS 2D
  pdf("plsplot 2d t1 vs t2.pdf",width=width,height=height)
  par(mar=c(5,5,4,2))
  plot(x,y,
       xlim = c(xmin,xmax),
       ylim = c(ymin,ymax),
       col = colour,
       pch = pcha,
       xlab = "t[1]",
       ylab = "t[2]",
       cex = cexa,
       cex.axis = cexaxis,
       cex.lab = cexlab)
  abline(v = 0, lty = 2)
  abline(h = 0, lty = 2)
  if (text) {text(x,y,rownames(int),pos=4)}
  if (ellipse) {lines(ellipse::ellipse(0,scale=c(sd(x),sd(y)),centre=c(mean(x),mean(y))),lty=2)}

  if (QC) {legend("topleft",
                  c(names(info),"QC"),
                  pch = pchalist[1:(length(info)+1)],
                  col = colourlist[1:(length(info)+1)],
                  bty="n",
                  cex=cexlegend,
                  pt.cex = cexa)
  }else{legend("topleft",names(info),
               pch=pchalist[1: length(info)],
               col=colourlist[1:length(info)],
               bty="n",
               cex=cexlegend,
               pt.cex = cexa)}
  dev.off()

  #PLS 3D
  if (number>2) {
    pdf("plsplot 3d.pdf",width=width,height=height)
    par(mar=c(5,5,4,2))
    scatterplot3d(x,y,z,color=colour,xlab="t[1]",ylab="t[2]",zlab="t[3]",angle=50,
                  pch=pcha,box=FALSE,cex.symbol=cexa,cex.lab=1.3,cex.axis=1.3,
                  xlim=c(xmin,xmax),ylim=c(ymin,ymax),zlim=c(zmin,zmax))
    if (QC) {
      legend("topleft",c( names(info),"QC"),
             pch = pchalist[1:(length(info)+1)],
             col = colourlist[1:(length(info)+1)],
             bty="n",
             cex=cexlegend,
             pt.cex = cexa)
    }else{
      legend("topleft",
             names(info),
             pch = pchalist[1:length(info)],
             col = colourlist[1:length(info)],
             bty="n",
             cex=cexlegend,
             pt.cex = cexa)
    }
    dev.off()
  }

  cat("PLS analysis is done\n")
}
jaspershen/SXTpls documentation built on May 18, 2019, 5:56 p.m.