bcm-uga / penda

PenDA is a rank-based method for personalized differential analysis.

penda

An R package that performs personalized differential analysis of omics data. The penda package provide different methods to perform differential analysis.

PenDA, a rank-based method for personalized differential analysis: Application to lung cancer Magali Richard, Clémentine Decamps, Florent Chuffart, Elisabeth Brambilla, Sophie Rousseaux, Saadi Khochbin, Daniel Jost. PLoS Comput Biol, 2020 May 11;16(5):e1007869. https://doi.org/10.1371/journal.pcbi.1007869.

Installation

To get the current development version from github:

git clone https://github.com/bcm-uga/penda.git
cd penda 
R

Build package

install.packages("devtools")
# mixtools htmltools scales yaml lazyeval plyr rlang ggplot2 gtools caTools KernSmooth
devtools::load_all(); devtools::document(); devtools::install()

Build vignettes

setwd("vignettes")
rmarkdown::render("vignette_simulation.Rmd")
# source(knitr::purl(("vignette_simulation.Rmd")))

#Depends of your number of controls:
rmarkdown::render("vignette_penda.Rmd")
rmarkdown::render("vignette_penda_1ctrl.Rmd")

penda's Issues

When using penda to draw ROC curve, an error occurred

Question:
I draw the ROC curve according to the method in the PENDA paper, but the maximum value of FPR in the ROC curve is only 0.3, and the ROC curve cannot be obtained (1,1).
Why is this? The following is the data and code I used.

Result

Paper
PenDA, a rank-based method for personalized differential analysis: Application to lung cancer

Data
tcga_data_ctrl.rds
tcga_data_case.rds

Code

library(penda)
library(data.table)

tcga_data_ctrl <- readRDS("C:/Users/liuyachen02/Desktop/rworkspace/data/penda_data/tcga_data_ctrl.rds")
tcga_data_case <- readRDS("C:/Users/liuyachen02/Desktop/rworkspace/data/penda_data/tcga_data_case.rds")
SIMU_NUM = 10


if(SIMU_NUM < dim(tcga_data_ctrl)[2]){
  sample_list = sample(dim(tcga_data_ctrl)[2], SIMU_NUM)
}else{
  print("ERROR：仿真样本的数量（SIMU_NUM）应该小于正常样本的数量")
}

data_simu = tcga_data_ctrl[,sample_list]
data_ctrl = tcga_data_ctrl[,-sample_list]

threshold_dataset = 0.99
Penda_dataset = penda::make_dataset(data_ctrl, 
                                    tcga_data_case, 
                                    detectlowvalue = TRUE, 
                                    detectNA = TRUE, 
                                    threshold = threshold_dataset)
data_ctrl <- Penda_dataset$data_ctrl
data_case <- Penda_dataset$data_case
data_simu <- data_simu[rownames(data_ctrl),]

S_MAX = 30
threshold_LH = 0.99
L_H_list <- penda::compute_lower_and_higher_lists(data_ctrl, 
                                                  threshold = threshold_LH, 
                                                  s_max = S_MAX)

L <- L_H_list
H <- L_H_list

DEREGULARATED_PROP = 0.3

simulation = penda::simplified_simulation(data = data_simu, 
                                          proportion = DEREGULARATED_PROP, 
                                          threshold = 100,
                                          modifier = 50,
                                          factor = 4)

quantile_values = c(0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.15, 0.2, 0.3, 0.4)
factor_values = c(1, 1.2, 1.4)
which_quantile = penda::choose_quantile(data_ctrl, 
                                        simulation, 
                                        factor_values = factor_values, 
                                        quantile_values = quantile_values)

best_quantile = penda::select_quantile_param(which_quantile, FDR_max = 0.05)
best_quant = best_quantile$quantile
best_fact = best_quantile$factor

threshold_values = seq(0.05, 0.95, 0.05)
which_threshold = penda::choose_threshold(data_ctrl, 
                                          L_H_list, 
                                          iterations = 30, 
                                          simulation, 
                                          threshold_values, 
                                          best_quant, 
                                          best_fact)

plot(which_threshold[,'FPR'], which_threshold[,'TPR'])