Introduction: le jeu de Den Boer 2009

Comme cas d’étude, nous analyserons un jeu de données de Den Boer et collègues [@DenBoer:2009ik].

Les auteurs ont caractérisé le transcriptome de 190 patients souffrand de leucémie lymphoblastique aigue (Acute Lymphoblastic Leukaemia, ALL).

Cette leucémie est caractérisée par une prolifération clonale anormale, dans la moëlle osseuse, de progéniteurs de lymphocytes qui restent bloqués dans une phase précise de leur différenciation. La prolifération peut résulter de différentes mutations, dont certaines sont récurrentes chez différents patients. Le jeu de données de Den Boer inclut l’information concernant les problèmes génétiques particuliers associés à chaque échantillon d’un patient atteint de leucémie.

Dans leur article original, Den Boer et collègues ont défini une signature transcriptomique, c’est-à-dire une sélection de gènes dont les profils de transcription sont spécifiquement affectés dans une ou plusieurs sous-classes des leucémies. Les profils d’expression de ces gènes-signature peuvent ensuite être utilisés pour classer de nouveaux échantillons de patients atteints de leucémie, et prédire le type de perturbation génétique associée à leur cas.

Un avantage de ce jeu de données est qu’il inclut des groupes d’effectifs relativement élevés (4 groupes avec >30 patients), par rapport à d’autres publications dans le domaine de l’analyse transcriptomique.

Le tableau de donnée sur lequel nous travaillerons comporte 190 colonnes (une par patient) et ~22.000 lignes (une par gène) et les valeurs indiquent le niveau d’expression du gène chez le patient considéré.

On dispose également d’une table “pheno” qui fournit des informations sur chaque échantillon, et en particulier la classe de leucémie dont souffre ce patien (en fait cette classe est déterminée par son génotype plutôt que son phénotype, mais pour des raisons historiques les tableaux annotations des échantillons sont communément dénommés “pheno” dans le domaine des biopuces).

But du TP

Le but de ce TP sera de mettre en place une procédure qui combinera les tâches suivants:

  1. Télécharger les fichiers de données à partir d’un site Web.
  2. Mesurer une série de paramètres statistiques sur les gènes (lignes) et échantillons (colonnes).
  3. Détecter les gènes différentiellement exprimés entre deux groupes d’échantillons (correspondant à deux classes de leucémie).

Nous essaierons pour chaque étape - de modulariser le code en écrivant des fonctions ré-utilisables - d’utiliser des techniques efficaces en évitant les boucles, et en utilisant plutôt la fonction apply()

Exercices

Créaction d’un dossier local pour les données et résultats

  1. Ecrivez un bloc de code R qui définit (dans une variable dénommé dir.output) le chemin d’un dossier où vous sauvegarderez les données et les résultats de ce TP. Ce dossier sera défini comme un sous-dossier STAT/TP/DenBoer2009 à partir de la racine de votre compte. Créez ce dossier (sauf s’il existe déjà) et vérifiez son contenu en fichiers.

    Fonctions utiles:

    • Sys.getenv("HOME") indiquera la racine de votre compte (Unix et Mac OS X)
    • file.path() permet de construire un chemin en concaténant des chaînes de caractères
    • dir.create() pour créer un dossier. Lisez l’aide pour choisir les options appropriées.
    • list.files() pour obtenir la liste des fichiers dans un dossier donné
## Define the directory for the practical
dir.home <- Sys.getenv("HOME")
dir.output <- file.path(dir.home, "STAT", "TP", "DenBOer2009")
message(dir.output)
## Create the directory
dir.create(dir.output, recursive = TRUE, showWarnings = FALSE)
## Change directory to dir.output
setwd(dir.output)
## Print a message with output directory
message("Output directory for this practical: ", dir.output)
## Check it this directory already contains some files.
list.files(dir.output)

Si c’est la première fois que vous réalisez ce TP, la réponse devrait être

character(0)

Ceci indique que le dossier que vous venez de créer ne contient aucun fichier. Le résultat de list.files() est un vecteur contenant 0 chaînes de caractères.

Téléchargement des données

Créez un sous-dossier data dans votre dossier TP, et téléchargez-y les deux fichiers de données ci-dessous, sauf si ces fichiers existent déjà à cet endroit.

URL du dossier des données de Den Boer:

Dans ce dossier nous utiliserons les deux tables suivantes.

Vous pouvez éventuellement utiliser le fichier de Description des sous-types de cancers (abréviation, couleurs d’affichage):

Le dossier contient également un fichier de qualité des mesures (absent / marginal / présent), que nous ignorerons pour ce TP-ci:

Vous trouverez ci-dessous le code pour télécharger une copie des fichiers utiles dans votre dossier personnel.

Le bloc de code ci-dessous définit l’URL du dossier de données, à partir duquel nous allons télécharger la table de transcriptome + la description des échantillons (pheno).

## Define the URL of the example data sets
url.course <- "http://pedagogix-tagc.univ-mrs.fr/courses/ASG1"
microarray.data.url <- file.path(url.course, "data", "marrays")
## Check the results
print(paste("Course URL: ", url.course))
[1] "Course URL:  http://pedagogix-tagc.univ-mrs.fr/courses/ASG1"
print(paste("URL for the microarray data folder: ", microarray.data.url))
[1] "URL for the microarray data folder:  http://pedagogix-tagc.univ-mrs.fr/courses/ASG1/data/marrays"

We will now load the expression matrix, which is stored in the file named GSE13425_Norm_Whole.txt in the microarray data folder of the server. This file contains a table with one row per probeset (~gene) and one column per sample. Each sample was taken from a patient suffering from a particular type of Acute Lymphoblastic Leukemia (ALL).

We will once use the following functions:

  • file.path() to build the full URL to the expression table
  • download.file() to download the file from the URL and obtain a local copy. This will avoid us to download the whole file each time we restart the analysis.
  • read.table() to read the table from this URL. Note that this functions allows you to load a table either from a local folder (on your computer), or from any Web site.
## Get some help about the read.table fonction
# help(read.table)
## URL of the expression table on the course server
expr.url <- file.path(microarray.data.url, "GSE13425_Norm_Whole.txt") ## Build the full URL of the matrix
## Load expression values
dir.data <- file.path(dir.output, "data")
dir.create(dir.data, recursive = TRUE, showWarnings = FALSE)
expr.file <- file.path(dir.data, "DenBoer_2009_GSE13425_Norm_Whole.txt")
if (file.exists(expr.file)) {
message("Expression file already downloaded: ", expr.file)
} else {
message("Downloading matrix table from ", expr.url)
download.file(url = expr.url, destfile = expr.file)
message("Local expression file: ", expr.file)
}
## Load the expression table in memory
message("Loading expression table from ", expr.file)
expr.matrix <- read.table(expr.file,sep="\t", header = T, row.names = 1)

Exercices

  1. Adaptez le code ci-dessus pour télécharger une copie locale du fichier pheno.
## Get some help about the read.table fonction
# help(read.table)
## URL of the pheno table on the course server
pheno.url <- file.path(microarray.data.url, "phenoData_GSE13425.tab") ## Build the full URL of the matrix
## Load pheno values
pheno.file <- file.path(dir.data, "phenoData_GSE13425.tab")
if (file.exists(pheno.file)) {
message("pheno file already downloaded: ", pheno.file)
} else {
message("Downloading matrix table from ", pheno.url)
download.file(url = pheno.url, destfile = pheno.file)
message("Local pheno file: ", pheno.file)
}
## Load the pheno table in memory
message("Loading pheno table from ", pheno.file)
pheno.matrix <- read.table(pheno.file,sep="\t", header = T, row.names = 1)
# View(pheno.matrix)
## Check that the order ofthe pheno table rows is identifcal to the order of the columns in the expression table. For this, we count the number of differences
sum(colnames(expr.matrix) != rownames(pheno.matrix))
[1] 0
  1. Chargez les fichiers d’expression et de description des échantillons dans des objets de type data.frame(), et récupérez dans un vecteur les sous-types de cancer associés à chaque échantillon.

  2. Comptez le nombre d’échantillons par sous-classe de cancer.

# View(pheno.matrix)
sample.labels <- as.vector(pheno.matrix$sample.labels)
class.sizes <- sort(table(sample.labels), decreasing = TRUE)
  1. Sélectionnez deux groupes comportant au moins 30 échantillons chacun, et créez une table d’expression et une table pheno limitées à ces deux groupes. A partir de ce moment vous travaillerez systématiquement avec ces deux tables.
## Select the labels of the two most populated classes
group1 <- names(class.sizes[1])
group2 <- names(class.sizes[2])
## Select the columns of the expression matrix for each group
group1.expr <- expr.matrix[, sample.labels == group1]
dim(group1.expr)
[1] 22283    44
group2.expr <- expr.matrix[, sample.labels == group2]
dim(group2.expr)
[1] 22283    44
selected.samples <- (sample.labels == group1) | (sample.labels == group2)
group.labels <- sample.labels[selected.samples]
  1. En utilisant la fonction apply(), calculez pour chaque gène les paramètres suivants, et stockez-les dans une data.frame().

    • \(m_1, m_2\) moyennes d’échantillons des deux groupes;
    • \(s_2.est, s_2.est\) estimations des écarts-types des populations dont sont extraits les échantillons du 1er et 2d groupe, respectivement;
    • \(d\) différence entre les moyennes des deux groupes.
gene.stats <- data.frame(
m1 = apply(group1.expr, 1, mean),
m2 = apply(group2.expr, 1, mean),
s1.pop = apply(group1.expr, 1, sd),
s2.pop = apply(group2.expr, 1, sd)
)
gene.stats$d <- gene.stats$m2 - gene.stats$m1
# View(gene.stats)
hist(gene.stats$d, breaks=100)

  1. Sélectionnez un gène au hasard, récupérez dans un vecteur la ligne correspondante dans la table d’expression, et effectuez un test de comparaison de moyennes entre les deux groupes.
## Number of genes
my.gene <- sample(x = row.names(expr.matrix), size = 1)
welch.result <- t.test(x = group1.expr[my.gene, ], y = group2.expr[my.gene,], alternative = "two.sided", var.equal = FALSE)
## Get all the attributes of the t.test() result
attributes(welch.result)
$names
[1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"    "null.value"  "alternative" "method"      "data.name"  

$class
[1] "htest"
result <- c(welch.result$statistic,
p = welch.result$p.value,
welch.result$parameter,
welch.result$estimate)
names(result) <- c("t", "p", "df", "m1", "m2")
result["d"] <- result["m2"] - result["m1"]
  1. Ecrivez une fonction R dénommée t.test.wrapper qui prend en entrée 2 vecteurs (expression et classes de cancer, respectivement), effectue un test de Student, et retourene un vecteur avec les statistiques pertinentes.
TTestVector <- function(values, group.labels, group1 = group.labels[1], ...) {
## Select samples of the first and second group, respectively
group1.values <- values[group.labels == group1]
group2.values <- values[group.labels != group1]
t.result <- t.test(x = group1.values, y = group2.values, ...)
## Collect the relevant statistics in the t.test result
result <- c(t.result$statistic,
p = t.result$p.value,
t.result$parameter,
t.result$estimate)
names(result) <- c("t", "p", "df", "m1", "m2")
result["d"] <- result["m2"] - result["m1"]
result["s1.est"] <- sd(group1.values)
result["s2.est"] <- sd(group2.values)
result["s.diff"] <- result["s2.est"] - result["s1.est"]
result["s.diff.percent"] <- 100*result["s.diff"] /max(result["s2.est"], result["s1.est"])
return(result)
}
## Run Welch test with my nice function
my.gene <- sample(x = row.names(expr.matrix), size = 1)
welch.result <- TTestVector(values = expr.matrix[my.gene,selected.samples],
group.labels = sample.labels[selected.samples],
alternative = "two.sided", var.equal = FALSE)
## Just for fun, run a Student test with the same gene
student.result <- TTestVector(values = expr.matrix[my.gene,selected.samples],
group.labels = sample.labels[selected.samples],
alternative = "two.sided", var.equal = TRUE)
## Compare the results obtained with Student and Welch t.test()
kable(data.frame(
Welch= welch.result,
Student = student.result), caption="Comparison between the results of a Welch and Student test on a randomly picked up gene. ")
Comparison between the results of a Welch and Student test on a randomly picked up gene.
Welch Student
t 2.6051051 2.6051051
p 0.0109993 0.0108210
df 77.8207622 86.0000000
m1 6.7281818 6.7281818
m2 6.2061364 6.2061364
d -0.5220455 -0.5220455
s1.est 0.7726887 0.7726887
s2.est 1.0816099 1.0816099
s.diff 0.3089213 0.3089213
s.diff.percent 28.5612459 28.5612459
  1. Ecrivez une boucle forqui applique ce même test à chaque ligne de la table d’expression, et récupère dans un tableau les résultats du test, avec une ligne par gène, une colonne par statistique pertinente produite par t.test(). Mesurez le temps d’exécution des 22.000 tests de Student avec cette boucle.
## Prepare a data frame to host all the results, with one row per gene
## and one column per field in the result of TTestVector
## (we use the previous result to adapt the number of columns.
welch.loop <- data.frame(matrix(nrow = nrow(expr.matrix), ncol=length(welch.result)))
names(welch.loop) <- names(welch.result) ## Set column names from previous result
row.names(welch.loop) <- row.names(expr.matrix) ## Set row names as gene names
head(welch.loop) ## Check that the data frame has the right shape and contains NA values
                  t  p df m1 m2  d s1.est s2.est s.diff s.diff.percent
DDR1|1007_s_at   NA NA NA NA NA NA     NA     NA     NA             NA
RFC2|1053_at     NA NA NA NA NA NA     NA     NA     NA             NA
HSPA6|117_at     NA NA NA NA NA NA     NA     NA     NA             NA
PAX8|121_at      NA NA NA NA NA NA     NA     NA     NA             NA
GUCA1A|1255_g_at NA NA NA NA NA NA     NA     NA     NA             NA
UBA7|1294_at     NA NA NA NA NA NA     NA     NA     NA             NA
## Run a
loop.time <- system.time(
for (my.gene in rownames(group1.expr)) {
welch.loop[my.gene,] <- TTestVector(
values=expr.matrix[my.gene, selected.samples],
group.labels = sample.labels[selected.samples],
alternative = "two.sided", var.equal = FALSE)
}
)
print(loop.time)
      user     system    elapsed 
    99.156     19.794 -59068.290 
  1. Utilisez cette fonction et la fonction apply pour appliquer le test de Student à chaque gène, et mesurz le temps d’exécution.
## Run all the t tests with the apply function
apply.time <- system.time(
welch.apply <- as.data.frame(t(apply(X = expr.matrix[, selected.samples], 1, TTestVector, group.labels = sample.labels[selected.samples], group1 = group1, alternative = "two.sided", var.equal = FALSE)))
)
print(apply.time)
   user  system elapsed 
  6.259   0.049   6.318 
kable(head(welch.apply), caption = "First rows of the table resulting from the application of Welch test to each gene. ")
First rows of the table resulting from the application of Welch test to each gene.
t p df m1 m2 d s1.est s2.est s.diff s.diff.percent
DDR1|1007_s_at -0.2044225 0.8386294 68.50664 4.679091 4.697045 0.0179545 0.2897458 0.5054428 0.2156970 42.674858
RFC2|1053_at -0.2171182 0.8286358 85.12354 3.132954 3.142727 0.0097727 0.2001230 0.2215733 0.0214503 9.680902
HSPA6|117_at -1.8028725 0.0758283 68.16601 2.563864 2.662955 0.0990909 0.1801836 0.3169446 0.1367610 43.149814
PAX8|121_at -0.5359874 0.5933759 84.46846 6.160682 6.179318 0.0186364 0.1517096 0.1737199 0.0220103 12.669978
GUCA1A|1255_g_at -1.2019848 0.2327599 83.68006 2.165000 2.188636 0.0236364 0.0842063 0.0996177 0.0154114 15.470531
UBA7|1294_at -6.0525510 0.0000000 77.01918 5.365682 5.854546 0.4888636 0.3074304 0.4387851 0.1313546 29.935985
## Nice formatting for the knit report (a bit tricky, because proc_time)
time.frame <- data.frame(
loop = as.factor(loop.time),
apply = as.factor(apply.time)
)
time.frame$loop <- as.numeric(as.vector(as.factor(loop.time)))
time.frame$apply <- as.numeric(as.vector(as.factor(apply.time)))
kable(time.frame, caption = paste("Time spent for", nrow(group1.expr), "Welch tests with apply:", apply.time), digits=3)
Time spent for 22283 Welch tests with apply: 6.259 Table: Time spent for 22283 Welch tests with apply: 0.049000000000003 Table: Time spent for 22283 Welch tests with apply: 6.31800000000658 Table: Time spent for 22283 Welch tests with apply: 0 Table: Time spent for 22283 Welch tests with apply: 0
loop apply
user.self 99.156 6.259
sys.self 19.794 0.049
elapsed -59068.290 6.318
user.child 0.000 0.000
sys.child 0.000 0.000

Exercice supplémentaire optionnel: à titre de curiosité, nous pouvons également appliquer le test de Student au même jeu de données (en combinant apply() et TTestVector()) et comparer (sur un graphe aux axes logarithmiques) les p-valeurs obtenues par Welch et Student, respectivement.

## Run all the t tests with the apply function
student.apply.time <- system.time(
student.apply <- as.data.frame(t(apply(X = expr.matrix[, selected.samples], 1, TTestVector, group.labels = sample.labels[selected.samples], group1 = group1, alternative = "two.sided", var.equal = TRUE)))
)
kable(head(student.apply), caption = "First rows of the table resulting from the application of Student test to each gene. ")
First rows of the table resulting from the application of Student test to each gene.
t p df m1 m2 d s1.est s2.est s.diff s.diff.percent
DDR1|1007_s_at -0.2044225 0.8385063 86 4.679091 4.697045 0.0179545 0.2897458 0.5054428 0.2156970 42.674858
RFC2|1053_at -0.2171182 0.8286305 86 3.132954 3.142727 0.0097727 0.2001230 0.2215733 0.0214503 9.680902
HSPA6|117_at -1.8028725 0.0749109 86 2.563864 2.662955 0.0990909 0.1801836 0.3169446 0.1367610 43.149814
PAX8|121_at -0.5359874 0.5933509 86 6.160682 6.179318 0.0186364 0.1517096 0.1737199 0.0220103 12.669978
GUCA1A|1255_g_at -1.2019848 0.2326688 86 2.165000 2.188636 0.0236364 0.0842063 0.0996177 0.0154114 15.470531
UBA7|1294_at -6.0525510 0.0000000 86 5.365682 5.854546 0.4888636 0.3074304 0.4387851 0.1313546 29.935985
print(apply.time)
   user  system elapsed 
  6.259   0.049   6.318 
plot(welch.apply[, "p"], student.apply[, "p"], log="xy")