Cuando necesitamos, para probar alguna herramienta bioinformática, datos de expresión de genes/proteínas que presenten diferencias significativas entre dos grupos de muestras, es decir, que tengan expresión diferencial (DE), podemos simularlos usando la función simData() del paquete {optBiomarker} de BioConductor.

Los grupos se denotan por H y D, las muestras se generan de manera única o se pueden replicar, se puede o no añadir expresión diferencial, se puede elegir la magnitud de esa expresión diferencial (fold change) y se pueden añadir variabilidades biológica y experimental.

La expresión diferencial se introduce añadiendo \(z·\delta\) a los datos del grupo 2, donde los valores \(\delta\) se generan a partir de una distribución normal truncada y \(z\) se selecciona aleatoriamente entre (-1,1) para caracterizar la sobre- o baja-expresión.

Para más detalle ver la ayuda: help('simData').

Ejemplo: Simulamos unos datos de expresión de 500 genes/proteinas/features en 40 muestras (10 del grupo H) replicadas 3 veces, con expresión diferencial de magnitud >= 4 (fold change mínimo) y con variabilidad experimental y biológica.

# Packages ----
if(!require('BiocManager')) {install.packages("BiocManager")}
if(!require('optBiomarker')) {BiocManager::install("optBiomarker")}
if(!require('pheatmap')) {install.packages("pheatmap")}
library(optBiomarker)
library(pheatmap)
# Simulate microarray data for two-group ----

# Parameters
nTrain <- 40      # total number of samples
nGr1 <- 10        # number of samples of group H
nBiom <- 500      # total number of genes/probes/proteins
nRep <- 3         # number of replicates of each sample
diffExpr <- TRUE  # introduce differential expression between the two groups
foldMin <- 4      # minimum value of fold changes between groups
sdW <- 2          # experimental (technical) variation (σ_e), log2 scale
sdB <- 1          # biological variation (σ_b), log2 scale

# Generate with simData function
data <- simData(nTrain = nTrain, nGr1 = nGr1, nBiom = nBiom, nRep = nRep,
                diffExpr = diffExpr, foldMin = foldMin, sdW = sdW, sdB = sdB,
                orderBiom = FALSE, baseExpr = NULL)

# Samples names
rownames(data$data) <- paste0('Sample',1:nTrain)

# View(data)

Podemos extraer información del ‘dataset’ data simulado, como por ejemplo los grupos de muestras groups y la expresión ex de los ‘features’ en ellas:

# Objects in data ----

# Groups
groups <- data$data[,1]
groups
 [1] H H H H H H H H H H D D D D D D D D D D D D D D D D D D D D D D D D D D D D
[39] D D
Levels: D H
table(groups)
groups
 D  H 
30 10 
prop.table(table(groups))
groups
   D    H 
0.75 0.25 
# Expression
ex <- data$data[,-1]
dim(ex)
[1]  40 500
summary(ex[,1:12])     # 12 features
   Biomarker1        Biomarker2        Biomarker3        Biomarker4     
 Min.   :-0.4063   Min.   :-0.6558   Min.   :-0.1167   Min.   :0.06227  
 1st Qu.: 1.3565   1st Qu.: 2.4971   1st Qu.: 2.5877   1st Qu.:2.08282  
 Median : 2.5264   Median : 3.5593   Median : 3.3573   Median :2.77575  
 Mean   : 2.8545   Mean   : 3.4804   Mean   : 3.6201   Mean   :3.15466  
 3rd Qu.: 3.9968   3rd Qu.: 4.4921   3rd Qu.: 4.7687   3rd Qu.:3.89671  
 Max.   : 7.7329   Max.   : 7.1315   Max.   : 8.5078   Max.   :7.27270  
   Biomarker5        Biomarker6      Biomarker7       Biomarker8    
 Min.   :-0.5439   Min.   :3.224   Min.   : 1.759   Min.   : 3.937  
 1st Qu.: 2.5462   1st Qu.:5.141   1st Qu.: 5.235   1st Qu.: 5.637  
 Median : 3.8457   Median :6.613   Median : 6.703   Median : 6.831  
 Mean   : 3.7688   Mean   :6.477   Mean   : 6.757   Mean   : 6.861  
 3rd Qu.: 5.2849   3rd Qu.:7.733   3rd Qu.: 8.249   3rd Qu.: 7.723  
 Max.   : 7.0455   Max.   :9.338   Max.   :10.596   Max.   :10.216  
   Biomarker9     Biomarker10      Biomarker11      Biomarker12   
 Min.   :2.821   Min.   : 2.935   Min.   : 2.386   Min.   :3.358  
 1st Qu.:5.612   1st Qu.: 5.337   1st Qu.: 6.068   1st Qu.:5.673  
 Median :6.327   Median : 6.331   Median : 6.751   Median :7.048  
 Mean   :6.386   Mean   : 6.569   Mean   : 6.795   Mean   :6.732  
 3rd Qu.:7.536   3rd Qu.: 8.037   3rd Qu.: 7.711   3rd Qu.:7.873  
 Max.   :9.225   Max.   :11.495   Max.   :11.140   Max.   :9.823  
summary(t(ex)[,1:8])   # 8 samples
    Sample1          Sample2          Sample3          Sample4     
 Min.   :0.3263   Min.   :0.2036   Min.   :0.8932   Min.   :1.117  
 1st Qu.:3.8148   1st Qu.:4.0307   1st Qu.:3.9230   1st Qu.:4.057  
 Median :4.8426   Median :5.0985   Median :5.0199   Median :4.985  
 Mean   :4.9059   Mean   :5.0462   Mean   :4.9859   Mean   :5.026  
 3rd Qu.:6.0093   3rd Qu.:6.1215   3rd Qu.:6.0911   3rd Qu.:6.028  
 Max.   :9.7553   Max.   :9.5783   Max.   :8.5971   Max.   :9.905  
    Sample5          Sample6          Sample7          Sample8      
 Min.   :0.7128   Min.   :0.2211   Min.   :0.8531   Min.   :0.6197  
 1st Qu.:3.8857   1st Qu.:4.0645   1st Qu.:3.9738   1st Qu.:4.0888  
 Median :5.0891   Median :5.0676   Median :4.9946   Median :5.0549  
 Mean   :5.0787   Mean   :5.0414   Mean   :5.0171   Mean   :5.0908  
 3rd Qu.:6.2239   3rd Qu.:6.0878   3rd Qu.:6.0450   3rd Qu.:6.1440  
 Max.   :9.2845   Max.   :9.4020   Max.   :9.2298   Max.   :9.5665

Con un boxplot visualizamos mejor la expresión y la variabilidad de todos los ‘features’ en cada muestra y con un heatmap visualizamos la expresión cruda de cada ‘feature’ en cada muestra, cómo se separan perfectamente los grupos y qué ‘features’ están bajo-, sobre- o no expresados.

# Samples boxplot
boxplot(t(ex), col = c('#00A9FF','#F8766D')[groups], las = 2, cex.axis = 0.75)

# Expression heatmap
anota <- data.frame(groups = groups)   # metadata dataframe for samples
rownames(anota) <- rownames(data$data)
ann_colors <- list(groups = c(D='#00A9FF',H='#F8766D'))
pheatmap(as.matrix(ex), annotation_row = anota, fontsize_row = 5, show_colnames = FALSE,
         annotation_colors = ann_colors, annotation_names_row = FALSE)