Exercices en R - Analyse des données de biopuces de Den Boer et al. (2009)

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

2. 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.

3. 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)

4. 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]

5. 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)

6. 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"] 

7. 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. ")

8. 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 
102.012  22.517 125.399 

9. 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 
  8.258   0.058   8.462 

kable(head(welch.apply), caption = "First rows of the table resulting from the application of Welch test to each gene. ")

## 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)

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. ")

print(apply.time)

   user  system elapsed 
  8.258   0.058   8.462 

plot(welch.apply[, "p"], student.apply[, "p"], log="xy")

10. Dessinez un volcano plot avec les résultats.

Le volcano plot indique en abscisse la "taille d'effet" (dans notre cas, la différence de moyennes), et l'ordonnée  $-log_{10}(p)$, où $p$ est la p-valeur.

## Define the threshold on p-value
alpha <- 0.05 ## Adopté à l'unanimité par les étudiants

## Draw a XY plot with the difference between means (X) versus the p-value
plot(welch.apply$d, welch.apply$p, xlab="m2 - m1", ylab="p", col="grey")
abline(v=0)
abline(h=alpha, col="red")

## Draw a volcano plot with Welch results
plot(welch.apply$d, -log10(welch.apply$p), xlab="m2 - m1", ylab="-log10(p)", col="grey", main="Volcano plot of the nominal p-values")
abline(v=0)
abline(h=-log10(alpha), col="red")

11. Dessinez un histogramme des p-valeurs.

## Histogram of p-values for Den Boer 2009 expression matrix
hist(welch.apply$p, breaks=seq(from=0, to=1, by=0.05), main="P-value histogram", xlab="P-value", ylab="Nb genes", las=1, col="#BBDDFF")

Exercices supplémentaires:

a. Générez une matrice de données aléatoires de la même taille que la matrice d'expression, où chaque ligne comporte des données normales sous hypothèse nulle ($H_0: \mu_1 = \mu_2$). Effectuez le test de Student et dessinez le volcano plot et l'histogramme des p-valeurs. A titre d'exemple, choisissons $\mu_1 = \mu_2 = 0$, $\sigma_1 = \sigma_2 = 1$. 
b. Générez de la même façon une autre matrice de données sous hypothèse alternative, avec $\mu_1 = 0$, $\mu_2 = 2$, $\sigma_1 = \sigma_2 = 1$, effectuez un test de Student et dessinez le volcano plot et l'histogramme des p-valeurs. 
c. Sélectionnez 15000 lignes de la première matrice (sous $H_0$), 5000 lignes de la seconde (sous $H_A$), effectuez les tests de Student et dessinez le volcano plot et l'histogramme des p-valeurs.

## Get the dimensions of the expresion matrix
g <- nrow(expr.matrix) ## number of genes
s <- ncol(expr.matrix) ## Number of samples

## Generate a matrix with random data under H0
rand.h0 <- matrix(nrow=g, ncol=s, data = rnorm(n = g * s))

student.under.h0 <- as.data.frame(t(apply(X = rand.h0[, selected.samples], 1, TTestVector, group.labels = sample.labels[selected.samples], group1 = group1, alternative = "two.sided", var.equal = TRUE)))

## Plot the histogram of p-values for random normal data under H0
hist(student.under.h0$p, breaks=seq(from=0, to=1, by=0.05), main="Random numbers under H0", xlab="P-value", ylab="Nb genes", las=1, col="#BBDDFF")

## Draw a volcano plot with random normal data under H0
plot(student.under.h0$d, -log10(student.under.h0$p), xlab="m2 - m1", ylab="-log10(p)", col="grey", main="Random normal under H0")
abline(v=0)
abline(h=-log10(alpha), col="red")

12. Dessinez un histogramme (en densités) avec les valeurs de la statistique du test (\(t\)).

hist(welch.apply$t, breaks=100)

13. Corrections de tests multiples. Ajoutez au tableau de résultats les colonnes suivantes.

    • e-valeur (\(E = p \cdot N\), où \(N\) est le nombre de test)
    • p-valeur corrigée par la méthode de Bonferonni, calculé avec la fonction p.adjust() (tond vous lirez préalablement l’aide).
    • taux de de fausse découverte (False Discovery Rate, FDR), calculé avec la fonction p.adjust().

## E-value
g <- nrow(expr.matrix)
welch.apply$E <- welch.apply$p * g
welch.apply$Bonferroni <- p.adjust(welch.apply$p, method = "bonferroni")
welch.apply$BH <- p.adjust(welch.apply$p, method = "BH")
welch.apply$fdr <- p.adjust(welch.apply$p, method = "fdr")
plot(welch.apply$p, welch.apply$Bonferroni, log="xy", col="red")
points(welch.apply$p, welch.apply$BH, col="darkgreen")
points(welch.apply$p, welch.apply$fdr, col="blue", pch=20, cex=0.8)
grid()

## Negative control: all data under H0
student.under.h0$E <- student.under.h0$p * g
student.under.h0$Bonferroni <- p.adjust(student.under.h0$p, method = "bonferroni")
student.under.h0$BH <- p.adjust(student.under.h0$p, method = "BH")
student.under.h0$fdr <- p.adjust(student.under.h0$p, method = "fdr")



## Draw a volcano plot with random normal data under H0
plot(student.under.h0$d, -log10(student.under.h0$E), xlab="m2 - m1", ylab="-log10(E)", col="grey", main="Random normal under H0")
abline(v=0)
abline(h=-log10(alpha), col="red")

## Histogram of E-values for Den Boer
hist(welch.apply$E, breaks=g *seq(from=0, to=1, by=0.05), main="E-value histogram", xlab="E-value", ylab="Nb genes", las=1, col="#BBDDFF")

## Plot the histogram of E-values for random normal data under H0
hist(student.under.h0$E, breaks=g*seq(from=0, to=1, by=0.05), main="Random numbers under H0", xlab="E-value", ylab="Nb genes", las=1, col="#BBDDFF")

par(mfrow=c(1,1))

14. Dessinez un graphe qui compare la p-valeur nominale (en abscisse) et les différentes corrections (en ordonnée). Utilisez les axes logarithmiques pour faire apparaître les valeurs pertinentes (très petites).

## Plot the volcano and p-value histograms of E-values
par(mfrow=c(2,2))
## Draw a volcano plot with Welch results
plot(welch.apply$d, -log10(welch.apply$E), xlab="m2 - m1", ylab="-log10(E)", col="grey", main="Den Boer, E-value volcano")
abline(v=0)
abline(h=-log10(alpha), col="red")## Histogram of p-values for Den Boer 2009 expression matrix

15. Sélectionez les gènes significatifs avec un seuil \(\alpha=0.05\) sur le taux de fausse découverte.

selected.genes <- welch.apply$fdr < alpha
sum(selected.genes)

[1] 2762

16. Sélectionnez une matrice restreinte aux profils d’expression des gènes différentiellement exprimés pour les deux groupes d’intérêt (les deux groupes dans la même matrice).

selected.expr.matrix <- expr.matrix[selected.genes, selected.samples]
dim(selected.expr.matrix)

[1] 2762   88

Ici s’arrête la matière à connaître pour l’examen: le reste du tuto est fourni pour le plaisir.

17. Calculez – avec la fonction cor() – une matrice de corrélation entre échantillons de cette matrice de gènes différentiellement exprimés. Générez un diagramme de cette matrice avec heatmap().

## Compute Pearson correlation between samples
dim(selected.expr.matrix)

[1] 2762   88

## Compute inter-sample correlation
sample.cor <- cor(selected.expr.matrix)
dim(sample.cor)

[1] 88 88

heatmap(sample.cor, scale = "none", col=bluered(100), symm=TRUE)

18. Utilisez la fonction hclust() pour effectuer un clustering des échantillons, sur base de cette matrice de corrélation.

## Compute inter-gene correlation
gene.cor <- 1 - cor(t(selected.expr.matrix))
dim(gene.cor)

[1] 2762 2762

## Build a gene tree by running hierarchical clustering between genes
gene.tree <- hclust(d = as.dist(gene.cor), method = "complete", members = NULL)
plot(gene.tree)

19. Avec la fonction heatmap(), dessinez une carte thermique (heat map) des gènes différentiellement exprimés, en activant le clustering sur les colonnes (samples) et lignes (gènes). Raffinez ensuite le graphique avec la fonction heatmap2() de la librairie gplots (installez cette librairie si nécessaire).

## Gene-wise centering for the heatmap. The gene-wise (row-wise) mean is substracted fromeach expression value
selected.expr.matrix.gene.centered <- selected.expr.matrix - apply(selected.expr.matrix,1, mean)


# library(corrplot)
# corrplot(sample.cor[1:10, 1:10])
x <- 
selected.expr.matrix.gene.centered

## Add class label to each gene name
colnames(x) <- paste(colnames(x), sample.labels[selected.samples])

## Compute maximal absolute value of x for the range of the heatmap
absmax <- ceiling(max(abs(range(x))))

## Use 99th percentile as limit for the color scale
zlim <- quantile(x = abs(unlist(x)), probs = 0.99)

## Draw a heatmap with one row per gene and one column per sample
heatmap(as.matrix(x), cexCol = 0.5, col=redgreen(100), scale="none", labRow = FALSE, symm=FALSE, margins=c(5,1), zlim=c(-zlim, zlim), Rowv = as.dendrogram(gene.tree))

20. Au moyen de la fonction strsplit(), extrayez les noms des gènes significatifs à partir des étiquettes de lignes de la matrice d’expression. Ces étiquettes contiennent 2 identifiants séparés par un caractère \(|\) (tube): nom de gène et identifiant de la sonde sur la biopuce.

gene.rownames <- rownames(expr.matrix)

gene.desc <- t(data.frame(strsplit(x = gene.rownames, split = '|', fixed=TRUE)))
rownames(gene.desc) <- rownames(expr.matrix)
colnames(gene.desc) <- c("gene.name", "affy.ID")
# dim(gene.desc)
kable(head(gene.desc))

## Selected gene names
selected.gene.names <- as.vector(gene.desc[selected.genes, "gene.name"])

selected.affy.IDs <- as.vector(gene.desc[selected.genes, "affy.ID"])

21. Analysez l’enrichissement fonctionnel de ces gènes avec les deux outils suivants:

    • gProfiler: http://biit.cs.ut.ee/gprofiler/
    • DAVID: https://david.ncifcrf.gov/

22. ACP (PCA), sparse PCA