Commit 57abad16 authored by smorabit's avatar smorabit
Browse files

added support to automatically select soft power

parent 7e7167a9

DESCRIPTION

+0 −5
Original line number Diff line number Diff line
@@ -13,8 +13,3 @@ LazyData: true 
Roxygen: list(markdown = TRUE)

RoxygenNote: 7.1.2

URL: https://smorabit.github.io/scWGCNA/
 
Imports:

  Seurat (>= 4.0),

  WGCNA,

  igraph,

  dplyr

NAMESPACE

+3 −0
Original line number Diff line number Diff line
@@ -34,6 +34,8 @@ export(OverlapModulesDEGs) 
export(OverlapModulesMotifs)

export(PlotDendrogram)

export(PlotModulePreservation)

export(PlotModuleTraitCorrelation)

export(PlotSoftPowers)

export(ProjectModules)

export(ROCCurves)

export(ResetModuleColors)
@@ -49,6 +51,7 @@ export(SetDatExpr) 
export(SetMultiExpr)

export(SetupForWGCNA)

export(TestSoftPowers)

export(TestSoftPowersConsensus)

export(TransferModuleGenome)

export(construct_metacells)

export(metacells_by_groups)

R/WGCNA_functions.R

+30 −99
Original line number Diff line number Diff line
@@ -141,23 +141,22 @@ SetupForWGCNA <- function( 


#' TestSoftPowers

#'

#' This function gets the expression matrix from the metacell object.

#' Compute the scale-free topology model fit for different soft power thresholds

#'

#' @param seurat_obj A Seurat object

#' @param powers

#' @param outfile filepath for output pdf to be generated. Default is

#' @param figsize numeric determining the height and width of the output figure. Default is c(7,7)

#' @param powers numeric vector specifying soft powers to test

#' @param setDatExpr logical flag indicating whether to run setDatExpr

#' @param ... additional parameters passed to SetDatExpr

#' @keywords scRNA-seq

#' @export

#' @examples

#' TestSoftPowers(pbmc)

TestSoftPowers <- function(

  seurat_obj,

  use_metacells = TRUE,

  group.by=NULL, group_name=NULL,

  setDatExpr = TRUE,

  powers=c(seq(1,10,by=1), seq(12,30, by=2)),

  make_plot=TRUE, outfile="softpower.pdf", figsize=c(7,7)

  setDatExpr = TRUE,

  use_metacells = TRUE,

  group.by=NULL, group_name=NULL

){



  # add datExpr if not already added:
@@ -179,44 +178,6 @@ TestSoftPowers <- function( 
    )[[2]]

  );



  # Plot the results:

  if(make_plot){

    pdf(outfile, height=figsize[1], width=figsize[2], useDingbats=FALSE)



        colors = c("blue", "red","black")

        # Will plot these columns of the returned scale free analysis tables

        plotCols = c(2,5,6,7)

        colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "Mean connectivity",

        "Max connectivity");



        # Get the minima and maxima of the plotted points

        ylim = matrix(NA, nrow = 2, ncol = 4);

        for (col in 1:length(plotCols)){

          ylim[1, col] = min(ylim[1, col], powerTable$data[, plotCols[col]], na.rm = TRUE);

          ylim[2, col] = max(ylim[2, col], powerTable$data[, plotCols[col]], na.rm = TRUE);

        }



        # Plot the quantities in the chosen columns vs. the soft thresholding power

        par(mfcol = c(2,2));

        par(mar = c(4.2, 4.2 , 2.2, 0.5))

        cex1 = 0.7;



        for (col in 1:length(plotCols)){

          plot(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],

          xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],

          main = colNames[col]);

          addGrid();



          if (col==1){

            text(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],

            labels=powers,cex=cex1,col=colors[1]);

          } else

          text(powerTable$data[,1], powerTable$data[,plotCols[col]],

          labels=powers,cex=cex1,col=colors[1]);

        }

    dev.off()

  }



  # set the power table in Seurat object:

  seurat_obj <- SetPowerTable(seurat_obj, powerTable$data)
@@ -225,15 +186,26 @@ TestSoftPowers <- function( 






#' TestSoftPowersConsensus

#'

#' Compute the scale-free topology model fit for different soft power thresholds separately for each input dataset

#'

#' @param seurat_obj A Seurat object

#' @param powers numeric vector specifying soft powers to test

#' @param setDatExpr logical flag indicating whether to run setDatExpr

#' @param ... additional parameters passed to SetDatExpr

#' @keywords scRNA-seq

#' @export

#' @examples

#' TestSoftPowers(pbmc)

TestSoftPowersConsensus <- function(

  seurat_obj,

  powers=c(seq(1,10,by=1), seq(12,30, by=2)),

  setDatExpr = TRUE,

  use_metacells = TRUE,

  group.by=NULL, group_name=NULL,

  multi.group.by = NULL,

  multi_groups = NULL,

  setDatExpr = TRUE,

  powers=c(seq(1,10,by=1), seq(12,30, by=2)),

  make_plot=TRUE, outfile="softpower", figsize=c(7,7)

  multi_groups = NULL

){



  # add multiExpr if not already added:
@@ -271,43 +243,6 @@ TestSoftPowersConsensus <- function( 
    powerTable$data$group <- cur_group

    powerTables[[cur_group]] <- powerTable$data



    # Plot the results:

    if(make_plot){

      pdf(paste0(outfile, '_', cur_group, '.pdf'), height=figsize[1], width=figsize[2], useDingbats=FALSE)



          colors = c("blue", "red","black")

          # Will plot these columns of the returned scale free analysis tables

          plotCols = c(2,5,6,7)

          colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "Mean connectivity",

          "Max connectivity");



          # Get the minima and maxima of the plotted points

          ylim = matrix(NA, nrow = 2, ncol = 4);

          for (col in 1:length(plotCols)){

            ylim[1, col] = min(ylim[1, col], powerTable$data[, plotCols[col]], na.rm = TRUE);

            ylim[2, col] = max(ylim[2, col], powerTable$data[, plotCols[col]], na.rm = TRUE);

          }



          # Plot the quantities in the chosen columns vs. the soft thresholding power

          par(mfcol = c(2,2));

          par(mar = c(4.2, 4.2 , 2.2, 0.5))

          cex1 = 0.7;



          for (col in 1:length(plotCols)){

            plot(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],

            xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],

            main = colNames[col]);

            addGrid();



            if (col==1){

              text(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],

              labels=powers,cex=cex1,col=colors[1]);

            } else

            text(powerTable$data[,1], powerTable$data[,plotCols[col]],

            labels=powers,cex=cex1,col=colors[1]);

          }

      dev.off()

    }

  }



  # merge the power tables
@@ -325,26 +260,22 @@ TestSoftPowersConsensus <- function( 
#' This function constructs a co-expression network from a Seurat object

#'

#' @param seurat_obj A Seurat object

#' @param soft_power

#' @param

#' @param

#' @param

#' @param

#' @param

#' @param

#' @param

#' @param

#' @param soft_power the soft power used for network construction. Automatically selected by default.

#' @param tom_outdir path to the directory where the TOM will be written

#' @param ... additional parameters passed to SetDatExpr and blockwiseConsensusModules

#' @keywords scRNA-seq

#' @export

#' @examples

#' ConstructNetwork(pbmc)

ConstructNetwork <- function(

  seurat_obj, soft_power=NULL, use_metacells=TRUE,

  setDatExpr=TRUE, group.by=NULL, group_name=NULL,

  seurat_obj, soft_power=NULL,

  tom_outdir="TOM",

  use_metacells=TRUE,

  setDatExpr=TRUE, group.by=NULL,

  group_name=NULL,

  consensus = FALSE,

  multi.group.by = NULL,

  multi_groups = NULL,

  tom_outdir="TOM",

  blocks=NULL, maxBlockSize=30000, randomSeed=12345, corType="pearson",

  consensusQuantile=0.3, networkType = "signed", TOMType = "unsigned",

  TOMDenom = "min", scaleTOMs = TRUE, scaleQuantile = 0.8,

R/plotting.R

+183 −0
Original line number Diff line number Diff line 
#' PlotSoftPowers

#'

#' Plot Soft Power Threshold results

#'

#' @param seurat_obj A Seurat object

#' @param selected_power power to highlight in the plots

#' @param point_size the size of the points in the plot

#' @param text_size the size of the text in the plot

#' @param plot_connectivity logical indicating whether to plot the connectivity in addition to the scale free topplogy fit.

#' @param wgcna_name The name of the WGCNA experiment in seurat_obj

#' @keywords scRNA-seq

#' @export

#' @examples

#' PlotSoftPowers

PlotSoftPowers <- function(

  seurat_obj,

  selected_power = NULL,

  point_size = 5,

  text_size=3,

  plot_connectivity = TRUE,

  wgcna_name = NULL

){



  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  pt <- GetPowerTable(seurat_obj, wgcna_name)



  if("group" %in% colnames(pt)){

    print("here")

    # select soft power for each group:

    power_tables <- pt %>% dplyr::group_split(group)

    soft_powers <- sapply(power_tables, function(power_table){

      power_table %>% subset(SFT.R.sq >= 0.8) %>% .$Power %>% min

    })



  } else{



    # get soft power:

    if(is.null(selected_power)){

      soft_power <- pt %>% subset(SFT.R.sq >= 0.8) %>% .$Power %>% min

    } else{

      soft_power <- selected_power

    }

    soft_powers <- NULL

    power_tables <- list("power" = pt)



  }



  plot_list <- list()



  for(i in 1:length(power_tables)){



    pt <- power_tables[[i]]



    if(!is.null(soft_powers)){

      soft_power <- soft_powers[i]

      print(i)

      print(soft_power)



      pt <- pt %>% dplyr::select(-group)

    }



    print(head(pt))



    # get other params:

    sft_r <- as.numeric(pt[pt$Power == soft_power,'SFT.R.sq'])

    mean_k <- as.numeric(pt[pt$Power == soft_power,'mean.k.'])

    median_k <- as.numeric(pt[pt$Power == soft_power,'median.k.'])

    max_k <- as.numeric(pt[pt$Power == soft_power,'max.k.'])



    # set color of text

    pt$text_color <- ifelse(

      pt$Power == soft_power, 'white', 'black'

    )



    # plot for soft power thresh:

    p1 <- pt %>% ggplot(aes(x=Power, y=SFT.R.sq)) +

      geom_rect(

        data = pt[1,],

        aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0.8), fill='grey80', alpha=0.8, color=NA

      ) +

      geom_hline(yintercept = sft_r, linetype='dashed') +

      geom_vline(xintercept = soft_power, linetype= 'dashed') +

      geom_point(

        data = pt[pt$Power == soft_power,c('Power', 'SFT.R.sq')],

        aes(x=Power, y=SFT.R.sq),

        inherit.aes=FALSE,

        color = 'black',

        size=point_size

      ) +

      geom_text(label=pt$Power, color = pt$text_color, size=text_size) +

      scale_y_continuous(limits = c(0,1), breaks=c(0, 0.2, 0.4, 0.6, 0.8, 1)) +

      ylab('Scale-free Topology Model Fit') +

      xlab('Soft Power Threshold') +

      theme(

        axis.line.x = element_blank(),

        axis.line.y = element_blank(),

        panel.border = element_rect(colour = "black", fill=NA, size=1)

      )



      if(plot_connectivity){



        # plot for mean connectivity:

        p2 <- pt %>% ggplot(aes(x=Power, y=mean.k.)) +

          geom_hline(yintercept = mean_k, linetype='dashed') +

          geom_vline(xintercept = soft_power, linetype= 'dashed') +

          geom_point(

            data = pt[pt$Power == soft_power,c('Power', 'mean.k.')],

            aes(x=Power, y=mean.k.),

            inherit.aes=FALSE,

            color = 'black',

            size=point_size

          ) +

          geom_text(label=pt$Power, color = pt$text_color, size=text_size) +

          scale_y_continuous(labels=scales::comma) +

          ylab('Mean Connectivity') +

          xlab('Soft Power Threshold') +

          theme(

            axis.line.x = element_blank(),

            axis.line.y = element_blank(),

            panel.border = element_rect(colour = "black", fill=NA, size=1)

          )



          # plot for medianan connectivity:

          p3 <- pt %>% ggplot(aes(x=Power, y=median.k.)) +

            geom_hline(yintercept = median_k, linetype='dashed') +

            geom_vline(xintercept = soft_power, linetype= 'dashed') +

            geom_point(

              data = pt[pt$Power == soft_power,c('Power', 'median.k.')],

              aes(x=Power, y=median.k.),

              inherit.aes=FALSE,

              color = 'black',

              size=point_size

            ) +

            geom_text(label=pt$Power, color = pt$text_color, size=text_size) +

            scale_y_continuous(labels=scales::comma) +

            ylab('Median Connectivity') +

            xlab('Soft Power Threshold') +

            theme(

              axis.line.x = element_blank(),

              axis.line.y = element_blank(),

              panel.border = element_rect(colour = "black", fill=NA, size=1)

            )



          # plot for mean connectivity:

          p4 <- pt %>% ggplot(aes(x=Power, y=max.k.)) +

            geom_hline(yintercept = max_k, linetype='dashed') +

            geom_vline(xintercept = soft_power, linetype= 'dashed') +

            geom_point(

              data = pt[pt$Power == soft_power,c('Power', 'max.k.')],

              aes(x=Power, y=max.k.),

              inherit.aes=FALSE,

              color = 'black',

              size=point_size

            ) +

            geom_text(label=pt$Power, color = pt$text_color, size=text_size) +

            scale_y_continuous(labels=scales::comma) +

            ylab('Max Connectivity') +

            xlab('Soft Power Threshold') +

            theme(

              axis.line.x = element_blank(),

              axis.line.y = element_blank(),

              panel.border = element_rect(colour = "black", fill=NA, size=1)

            )





          plot_list[[i]] <- list(p1,p2,p3,p4)



      } else{

        plot_list[[i]] <- p1

      }



  }



  if(length(plot_list) == 1){

    print('return')

    return(plot_list[[1]])

  }

  plot_list



}







#' PlotDendrogram

#'

_pkgdown.yml

+13 −1
Original line number Diff line number Diff line
@@ -71,8 +71,15 @@ reference: 
  contents:

  - '`ConstructMetacells`'

  - '`MetacellsByGroups`'

- title: Network Analysis

  desc: Functions for constructing the co-expression network

  contents:

    - '`TestSoftPowers`'

    - '`TestSoftPowersConsensus`'

    - '`PlotSoftPowers`'

    - '`ConstructNetwork`'

- title: Network Visualization

  desc: Functions for visualizing the scWGCNA co-expression network

  desc: Functions for visualizing the co-expression network

  contents:

    - '`ModuleNetworkPlot`'

    - '`HubGeneNetworkPlot`'
@@ -96,6 +103,11 @@ reference: 
  contents:

  - '`ModulePreservation`'

  - '`PlotModulePreservation`'

- title: Module Trait Correlation

  desc: Functions for performing module trait correlation analysis

  contents:

  - '`ModuleTraitCorrelation`'

  - '`PlotModuleTraitCorrelation`'

- title: Seurat wrappers

  desc: Wrapper functions to run Seurat commands on the metacell data

  contents:

CRA Git | Maintained and supported by SUSTech CRA and CCSE