Le but de ce travail est de mobiliser les différentes méthodes statistiques enseignées durant les cours et travaux pratiques afin d’analyser un jeu de données, et d’interpréter les résultats.
Les énoncés et questions sont fournis sous forme d’un document R markdown, qui servira également de modèle pour répondre auw questions.
En-dessous de chaque question nous vous demandons de fournir d’une part le code qui permet d’y repondre, et une interprétation sommaire des résultats (5 à 15 lignes par question).
Le rapport de travail consistera à remettre un document Rmd permettant aux évaluateurs de reproduire l’ensemble des résultats, et le rapport pdf produit au moyen de ce markdown. Ces deux documents sont à mettre dans une archive (nom.tgz) et à déposer sur moodle avant le lundi 11 mars à 3h00. Cette date-limite permettra de s’assurer que vous vous soyez tous affranchis des contrôles et ayez approfondi votre pratique de R avant le début des stages. Elle est cependant négociable; revenir vers nous (Jacques ou Anne) pour en discuter.
Le jeu de données se base sur une analyse de 125 patients souffrant de leucémie myéloïde aiguë (AML), disponible sur le site du projet The Cancer Genome Atlas (TCGA). Pour ce travail-ci nous nous focaliserons sur les données transcriptomiques (RNA-seq). Notez que le site TCGA fournit aussi des données de méthylome pour les mêmes patients.
Les données de comptage et métadonnées ont été téléchargées de la base de données Recount2 via la librairie R recount, et pré-moulinées pour les fournir dans un format directement utilisable.
Le jeu de données fourni consiste en 2 tableaux
Expression (fichier AML_counts_filtered-genes.tsv.gz) : comptages bruts des comptages de lectures courtes (reads) par gène (lignes) pour chaque patient (colonnes). Les gènes ont été filtrés pour écarter les lignes contenant trop de zeros et celles dont toutes les valeurs étaient inférieures à 10. Vous ne devez donc plus effectuer cette étape de filtrage.
Classes de patients (fichier AML_sample-classes.tsv.gz) : le pronostic associé à chaque patient, avec trois valeurs: “Favorable”, “Poor”, “Intermediate/Normal”.
Nous fournissons également pour information un tableau complet des donnés phénotypiques (AML_pheno.tsv.gz) extraites des données TCGA via Recount, mais celui-ci n’est pas nécessaire à la réalisation du travail personnel.
Le code ci-dessous définit l’endroit où les données et résultats seront sauvegardés dans votre compte. Vous pouvez bien entendu le modifier.
workdir <- "~/DU-Bii/m3/TGCA_AML_analysis"
dir.create(path = workdir, showWarnings = FALSE, recursive = TRUE)
message("Working directory\t", workdir)Pour vous éviter de perdre trop de temps avec la localisation des fichiers de données, nous fournissons nous-même les instructions qui permettent de les charger.
## Set the URL of the data folder
data.folder.URL <- "https://du-bii.github.io/study-cases/Homo_sapiens/TCGA_study-case/TCGA_AML_data/"
## Download a local copy of the expression file if not already done
## This avoids wasting time to re-download all the data at each run of the script.
expr.file <- "AML_counts_filtered-genes.tsv.gz" # Expression file name
expr.URL <- file.path(data.folder.URL, expr.file) # URL of expression file on the server
expr.path <- file.path(workdir, expr.file) # Full path of the local copy of the expression file
if (!file.exists(expr.path)) {
message("Downloading expression file from ", expr.URL)
download.file(url = expr.URL, destfile = expr.path)
}
## Load expression table
expr.table <- read.delim(
file = expr.path,
header = TRUE,
row.names = NULL)
message("Loaded expression table with " ,
nrow(expr.table), " rows (genes) x ",
ncol(expr.table), " columns")Le tableau d’expression contient
## Download a local copy of the sample classes if not already done
## This avoids wasting time to re-download all the data at each run of the script.
sample.class.file <- "AML_sample-classes.tsv.gz" # sample.class file name
sample.class.URL <- file.path(data.folder.URL, sample.class.file) # URL of sample.class file on the server
sample.class.path <- file.path(workdir, sample.class.file) # Full path of the local copy of the sample.class file
if (!file.exists(sample.class.path)) {
message("Downloading sample.class file from ", sample.class.URL)
download.file(url = sample.class.URL, destfile = sample.class.path)
}
## Load sample.class table
sample.class.table <- read.delim(
file = sample.class.path,
header = TRUE,
row.names = 1)
## Get the sample class in a vector
sample.classes <- as.vector(sample.class.table$sample.class)
## Set patient ID as names to each entry of this vector
names(sample.classes) <- row.names(sample.class.table)
## Report class sizes
message("Loaded sample.classes in a vector names sample.classes ")
kable(table(sample.classes), caption = "Samples per classe")| sample.classes | Freq |
|---|---|
| Favorable | 25 |
| Intermediate/Normal | 70 |
| Poor | 30 |
## Download a local copy of the pheno table if not already done
## This avoids wasting time to re-download all the data at each run of the script.
pheno.file <- "AML_pheno.tsv.gz" # pheno file name
pheno.URL <- file.path(data.folder.URL, pheno.file) # URL of pheno file on the server
pheno.path <- file.path(workdir, pheno.file) # Full path of the local copy of the pheno file
if (!file.exists(pheno.path)) {
message("Downloading pheno file from ", pheno.URL)
download.file(url = pheno.URL, destfile = pheno.path)
}
## Load pheno table
pheno.table <- read.delim(
file = pheno.path,
header = TRUE,
row.names = NULL)
message("Loaded pheno table with " ,
nrow(pheno.table), " rows (genes) x ",
ncol(pheno.table), " columns")## Histogram of the raw counts## histogram of log2-transformed countssample()), et dessinez un box plot des comptages bruts et log2-transformés. Commentez la figure en quelques lignes.## Box plot of the raw counts## Box plot of the log2-transformed counts## Draw PC plots for PC1 vs PC2, PC3 vs PC4 and PC5 vs PC6## Selection of differentially expressed genes with edgeR## P-value histogram## MA plot## Volcano plotAvec la fonction write.table(), exportez dans votre espace de travail
## Export DEG full table
## Export DEG table for genes declared positive only
## Export text file with DEG names (one gene name per line)Effectuez une analyse de classification supervisée en comparant soit deux méthodes différentes (par exemple Random Forest versus SVM) soit deux valeurs de parameètres pour une méthode donnée (par exemple les différents types de noyaux possibles pour SVM).
Pour cela, utilisez une stratégie de sous-échantilonnage aléatoire en utilisant 2/3 des échantillons pour l’apprentissage et 1/3 pour l’évaluation. Répétez l’opération 10 fois (par exemple au moyen d’une boucle for), et retentez à chaque itération l’indice de performance obtenu (par exemple “Accuracy”: somme des classifications correctes / total).
Pour chaque méthode ou valeur de paramètre, effectuez une analyse en multi-groupes (avec les 3 groupes d’AML), calculez la table de confusion, le taux de prédictions correctes (Accuracy) et d’erreurs de classification (Misclassification Error Rate, R).
Commentez les résultats.
## Multi-groups Supervised classification