R/multipleRegion_plot.R

In neuroconductor-devel-releases/brainKCCA: Region-Level Connectivity Network Construction via Kernel Canonical Correlation Analysis

#' Multiple Regions Connection Plot
#'
#' This function can draw connection line(s) among multiple brain regions,
#' based on the result from funciton permkCCA_multipleRegion(...).
#' @author Xubo Yue, Chia-Wei Hsu (tester), Jian Kang (maintainer)
#' @param input it can accept two types of argument:
#' (1) output generated by "permkCCA_multipleRegion(...)" (with p value).
#' (2) output generated by summary_kcca() or meanConnection() (without p value).
#' @param significance significance level of hypothesis testing.
#' If there is no p value in your data (i.e. only region name and
#' region index, then set significance as NA)
#' @param regionCodeProvided use TRUE if region code was provided in nii2RData()
#' funciton and FALSE if not.
#' @param view view of brain to see. One can choose coronal, axial, SL, or SR.
#' @param color color of dots and lines. User can specify multiple colors in this
#' argument and the function will use these colors periodically.
#' @param screenShot whether to save figure as png. If it is "None", then the
#' figure will not be saved; If it is a name, for example, "plot.png",
#' then it will be saved as "plot.png".
#' @return 1-4 view(s) (depends on the number of views user would like to see)
#' of brain graphic network. All significant regions are connected by colored lines.
#' @details If you provide an input with multiple patients, the function will ask
#' you to specify which patients data you would like to use.
#' @importFrom rgl rgl.snapshot
#' @export
#' @examples
#' \donttest{
#' #It will take more than 3 min to run
#' #write data into temp file
#' filePath <- tempdir()
#' #the nii.gz fMRI imaging file is created (toy example)
#' oro.nifti::writeNIfTI(brainKCCA::input_img, paste(filePath, "/",  "temp", sep=""))
#' #read fMRI data
#' testcase1 <- nii2RData(niiFile1 = "temp", resolution = "3mm", imgPath = filePath)
#' result1<-permkCCA_multipleRegion(imageDat = testcase1, region = c(1,5,10))
#' multipleRegion_plot(result1, view="axial")
#' }

multipleRegion_plot = function(input, significance=0.05, regionCodeProvided=FALSE,
                               view="coronal", color="blue", screenShot="None"){

  template <- oro.nifti::readNIfTI(system.file("MNI152_T1_1mm_brain.nii.gz", package = "brainR"),reorient = FALSE)

  if(length(view)<1) {warning("length of your vector should be greater than 0, please check your input vector!", call. = FALSE); return(0)}
  if(length(view)>4) {warning("length of your vector should be less than 5, please check your input vector!", call. = FALSE); return(0)}
  if(sum(table(view)>1)) {warning("Same view appears more than once, please check your input vector!", call. = FALSE); return(0)}

  #check which patients' data is of interest
  input_length <- length(input)-1
  if(input_length>=2){
    cat("\n", "You have ", input_length, " patients' data", "\n")
    index <- readline(prompt="Which data you would like to see? ")
    index <- as.numeric(index)
    if(index>input_length | index<0){warning("length of your input is out of bound.", call. = FALSE); return(0)}
    input<-input[[index]]
  }
  else input<-input[[1]]

  #check whether region code is provided
  if(regionCodeProvided!=FALSE){
    node <- input[[length(input)]]
  }

  #read brain template
  template <- oro.nifti::readNIfTI(system.file("MNI152_T1_1mm_brain.nii.gz", package = "brainR"),reorient = FALSE)


  ##############################################################################################
  #only significant region names and region index w/o p value
  ##############################################################################################
  if(is.na(significance)){

    if(length(view)==1) rgl::mfrow3d(nr = 1, nc = 1, sharedMouse = TRUE)
    else if(length(view)==2) rgl::mfrow3d(nr = 1, nc = 2, sharedMouse = TRUE)
    else if(length(view)==3) rgl::mfrow3d(nr = 1, nc = 3, sharedMouse = TRUE)
    else if(length(view)==4) rgl::mfrow3d(nr = 2, nc = 2, sharedMouse = TRUE)

    if("coronal" %in% view) {rgl::par3d(userMatrix=coronal_view);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "coronal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("axial" %in% view) {rgl::next3d();rgl::par3d(userMatrix=axial_view);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "axial", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("SL" %in% view) {rgl::next3d();rgl::par3d(userMatrix=sagittal_view_L);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "Left Sagittal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("SR" %in% view) {rgl::next3d();rgl::par3d(userMatrix=sagittal_view_R);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "Right Sagittal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }

    for(i in 1:dim(input[[1]])[1]){

      for(j in 1:length(view)){
        rgl::next3d(clear = FALSE)
        rgl::plot3d(c(node[input[[1]][i,1],1]+90,node[input[[1]][i,2],1]+90),
                    c(node[input[[1]][i,1],2]+126,node[input[[1]][i,2],2]+126),
                    c(node[input[[1]][i,1],3]+72,node[input[[1]][i,2],3]+72),
                    add=TRUE,col=color[(i%%length(color))+1],size=5)

        rgl::segments3d(c(node[input[[1]][i,1],1]+90,node[input[[1]][i,2],1]+90),
                        c(node[input[[1]][i,1],2]+126,node[input[[1]][i,2],2]+126),
                        c(node[input[[1]][i,1],3]+72,node[input[[1]][i,2],3]+72)
                        ,col=color[(i%%length(color))+1],lwd=3)
      }

    }

    interest<-as.vector(as.numeric(names(table(input[[1]]))))
    for(j in 1:length(view)){
      rgl::next3d(clear = FALSE)
      for(i in 1:length(interest))
        rgl::text3d(node[interest[i],1]+90,node[interest[i],2]+126,node[interest[i],3]+72,node[interest[i],6])
    }
    rgl::par3d("windowRect"= c(0,50,1200,1200))
    if(screenShot==TRUE) rgl::rgl.snapshot("plot.png", fmt = "png", top = TRUE )
    return()
  }



  ###############################################
  #two region plot
  if(is.null(dim(input[[1]]))){

    if(length(view)==1) rgl::mfrow3d(nr = 1, nc = 1, sharedMouse = TRUE)
    else if(length(view)==2) rgl::mfrow3d(nr = 1, nc = 2, sharedMouse = TRUE)
    else if(length(view)==3) rgl::mfrow3d(nr = 1, nc = 3, sharedMouse = TRUE)
    else if(length(view)==4) rgl::mfrow3d(nr = 2, nc = 2, sharedMouse = TRUE)

    if("coronal" %in% view) {rgl::par3d(userMatrix=coronal_view);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::plot3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                  c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                  c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72),
                  add=TRUE,col=color,size=5)
      rgl::decorate3d(main = "coronal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("axial" %in% view) {rgl::next3d();rgl::par3d(userMatrix=axial_view);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::plot3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                  c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                  c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72),
                  add=TRUE,col=color,size=5)
      rgl::decorate3d(main = "axial", xlab="", ylab="", zlab="",  box=FALSE,axes=FALSE)
    }
    if("SL" %in% view) {rgl::next3d();rgl::par3d(userMatrix=sagittal_view_L);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::plot3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                  c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                  c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72),
                  add=TRUE,col=color,size=5)
      rgl::decorate3d(main = "Left Sagittal", xlab="", ylab="", zlab="",  box=FALSE,axes=FALSE)
    }
    if("SR" %in% view) {rgl::next3d();rgl::par3d(userMatrix=sagittal_view_R);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::plot3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                  c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                  c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72),
                  add=TRUE,col=color,size=5,xlab = "Right Sagittal")
      rgl::decorate3d(main = "Right Sagittal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }

    if(input[[2]]<significance){
      cat("p value of kccc between ", input[[4]][1], " and ", input[[4]][2], " is: ", input[[2]])
      for(i in 1:length(view)){
        rgl::next3d(clear = FALSE)
        rgl::text3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                    c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                    c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72),
                    node[input[[1]],6])
        rgl::segments3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                        c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                        c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72)
                        ,col=color[(i%%length(color))+1],lwd=2)
      }
    }
    else{
      #cat("p value of kccc between ", input[[4]][1], " and ", input[[4]][2], " is: ", input[[2]])
      for(i in 1:length(view)){
        rgl::next3d(clear = FALSE)
        rgl::text3d(c(node[input[[1]][1],1]+90,node[input[[1]][2],1]+90),
                    c(node[input[[1]][1],2]+126,node[input[[1]][2],2]+126),
                    c(node[input[[1]][1],3]+72,node[input[[1]][2],3]+72),
                    node[input[[1]],6])
      }
    }
    rgl::par3d("windowRect"= c(0,50,1200,1200))
    if(screenShot!="None") rgl.snapshot(screenShot, fmt = "png", top = TRUE )
  }

  #multiple region plot
  else{

    if(length(view)==1) rgl::mfrow3d(nr = 1, nc = 1, sharedMouse = TRUE)
    else if(length(view)==2) rgl::mfrow3d(nr = 1, nc = 2, sharedMouse = TRUE)
    else if(length(view)==3) rgl::mfrow3d(nr = 1, nc = 3, sharedMouse = TRUE)
    else if(length(view)==4) rgl::mfrow3d(nr = 2, nc = 2, sharedMouse = TRUE)

    if("coronal" %in% view) {rgl::par3d(userMatrix=coronal_view);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "coronal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("axial" %in% view) {rgl::next3d();rgl::par3d(userMatrix=axial_view);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "axial", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("SL" %in% view) {rgl::next3d();rgl::par3d(userMatrix=sagittal_view_L);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "Left Sagittal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }
    if("SR" %in% view) {rgl::next3d();rgl::par3d(userMatrix=sagittal_view_R);
      misc3d::contour3d(template, level = 4500,alpha = 0.2, draw = TRUE,add=TRUE)
      rgl::decorate3d(main = "Right Sagittal", xlab="", ylab="", zlab="", box=FALSE,axes=FALSE)
    }

    for(i in 1:length(input[[2]])){

      if(input[[2]][[i]]<0.05){
        cat("p value of kccc between ", input[[4]][i,1], " and ", input[[4]][i,2], " is: ", input[[2]][[i]], "\n")

        for(j in 1:length(view)){
          rgl::next3d(clear = FALSE)
          rgl::plot3d(c(node[input[[1]][i,1],1]+90,node[input[[1]][i,2],1]+90),
                      c(node[input[[1]][i,1],2]+126,node[input[[1]][i,2],2]+126),
                      c(node[input[[1]][i,1],3]+72,node[input[[1]][i,2],3]+72),
                      add=TRUE,col=color[(i%%length(color))+1],size=5)

          rgl::segments3d(c(node[input[[1]][i,1],1]+90,node[input[[1]][i,2],1]+90),
                          c(node[input[[1]][i,1],2]+126,node[input[[1]][i,2],2]+126),
                          c(node[input[[1]][i,1],3]+72,node[input[[1]][i,2],3]+72)
                          ,col=color[(i%%length(color))+1],lwd=3)
        }
      }

      else{
        cat("p value of kccc between ", input[[4]][i,1], " and ", input[[4]][i,2], " is: ", input[[2]][[i]], "\n")

        for(j in 1:length(view)){
          rgl::next3d(clear = FALSE)
          rgl::plot3d(c(node[input[[1]][i,1],1]+90,node[input[[1]][i,2],1]+90),
                      c(node[input[[1]][i,1],2]+126,node[input[[1]][i,2],2]+126),
                      c(node[input[[1]][i,1],3]+72,node[input[[1]][i,2],3]+72),
                      add=TRUE,col="blue",size=5)
        }
      }
    }

    interest<-as.vector(as.numeric(names(table(input[[1]]))))
    for(j in 1:length(view)){
      rgl::next3d(clear = FALSE)
      for(i in 1:length(interest))
        rgl::text3d(node[interest[i],1]+90,node[interest[i],2]+126,node[interest[i],3]+72,node[interest[i],6])
    }
    rgl::par3d("windowRect"= c(0,50,1200,1200))
    if(screenShot==TRUE) rgl::rgl.snapshot("plot.png", fmt = "png", top = TRUE )
  }
}