Nous allons utiliser le package R appelé DESeq2 afin d’identifier les gènes qui sont statistiquement différentiellement exprimés (DEGs) entre les différentes conditions expérimentales. Ce script a été mis au point dans le laboratoire du Prof. Hanikenne dans le cadre de l’analyse du mutant frd3 et adapté par le Dr. Bouché pour ce cours.
Afin de réaliser l’analyse au cours de cette séance de TP/TD, nous vous invitons à télécharger la matrice brute de compte (“raw count matrix”), qui est disponible dans la section 2. II (read count and mapping) de ce site web. Vous pourrez ensuite réaliser les différentes opérations décrites ci-dessous.
Des informations détaillées à propos du package DESeq2 sont disponibles ici :
https://bioconductor.org/packages/release/bioc/html/DESeq2.html.
DESeq2 permet d’estimer la dépendance entre la variance et la moyenne (variance-mean dependence) dans les données de comptage provenant d’expériences de séquençage à haut débit et de tester l’expression différentielle. Pour les données de comptage, comme celles du séquençage de l’ARN, cette relation est particulièrement importante. Dans DESeq2, cette dépendance entre la variance et la moyenne est modélisée parce que les données de comptage (comme les niveaux d’expression génique) montrent souvent plus de variabilité à des niveaux d’expression plus élevés : à mesure que les niveaux d’expression augmentent, les comptages deviennent plus dispersés (c’est-à-dire que la variance augmente avec la moyenne).
Cette dépendance est souvent modélisée à l’aide de la distribution binomiale négative dans DESeq2, ce qui permet à la variance d’augmenter avec la moyenne, et donc de mieux représenter la surdispersion couramment observée dans les données biologiques. Cela permet d’effectuer des tests statistiques plus précis pour détecter l’expression différentielle. Le modèle généralement employé se base sur les modèles linéaires généralisés (GLM) de distribution binomiale négative.
Un modèle GLM très basique est le suivant :
# y=B0+B1x
# where
# B0 is the intercept (where the line crosses the y axis)
# B1 is the coefficient of the x variable
DESeq2 va d’abord modéliser les données pour estimer les coefficients à employer dans le modèle. Une fois les coefficients établis, il est possible de déterminer y pour chaque x.
Voici une brève description des différentes fonctions et arguments que nous allons utiliser. Pour une description complète, consultez le manuel.
DESeqDataSet is a subclass of RangedSummarizedExperiment, used to store the input values, intermediate calculations and results of an analysis of differential expression. The DESeqDataSet class enforces non-negative integer values in the “counts” matrix stored as the first element in the assay list. In addition, a formula which specifies the design of the experiment must be provided.
Usage
$ DESeqDataSetFromMatrix(countData, colData, design)
Arguments
design a formula which expresses how the counts for each gene depend on the variables in colData. Many R formula are valid, including designs with multiple variables, e.g., ~ group + condition, and designs with interactions, e.g., ~ genotype + treatment + genotype:treatment. See results for a variety of designs and how to extract results tables.
countData for matrix input: a matrix of non-negative integers.
colData for matrix input: a DataFrame or data.frame with at least a single column. Rows of colData correspond to columns of countData.
DESeq performs the differential expression analysis based on the Negative Binomial distribution. It performs a default analysis through the steps:
For complete details on each step, see the manual pages of the respective functions. After the DESeq function returns a DESeqDataSet object, results tables (log2 fold changes and p-values) can be generated using the results function. See the manual page for information on independent filtering and p-value adjustment for multiple test corrections.
Usage
$ DESeq(object)
Arguments
object a DESeqDataSet object, e.g. DESeqDataSetFromMatrix.
results extracts a result table from a DESeq analysis giving base means across samples, log2 fold changes, standard errors, test statistics, p-values and adjusted p-values.
resultsNames returns the names of the estimated effects (coefficents) of the model.
Usage
$results(object, contrast, lfcThreshold = 0, alpha = 0.1)
$ resultsNames(object)
Arguments
object a DESeqDataSet, on which one of the following functions has already been called: DESeq, nbinomWaldTest, or nbinomLRT.
contrast this argument specifies what comparison to extract from the object to build a results table. In our case, a character vector with exactly three elements: the name of a factor in the design formula, the name of the numerator level for the fold change, and the name of the denominator level for the fold change (simplest case).
lfcThreshold a non-negative value which specifies a log2 fold change threshold. The default value is 0, corresponding to a test that the log2 fold changes are equal to zero.
alpha the significance cutoff used for optimizing the independent filtering (by default 0.1). If the adjusted p-value cutoff (FDR) will be a value other than 0.1, alpha should be set to that value.
plotCounts allows creating plot of normalized counts for a single gene on log scale.
Usage
$ plotCounts(dds, gene, intgroup = "condition")
Arguments
dds a DESeqDataSet.
gene a character, specifying the name of the gene to plot.
intgroup interesting groups: a character vector of names in colData(x) to use for grouping.
Dans la localisation qui vous convient le mieux (e.g. bureau, document, folder du cours, etc), créez un dossier intitulé “rna-seq” (en minuscule).
Dans celui-ci, créez les dossiers vides suivants :
Téléchargez la matrice de comptage brute (tophat_root.matrix) dans la section “Downloads”. Attention, veillez à ce que la matrice de compte possède bien l’extension “.matrix”. Si l’extension est “.txt”, il faut retirer celle-ci.
Ajoutez la matrice de compte au dossier “data” de votre dossier rna-seq
Installez la dernière version de RStudio depuis le site RStudio ;
Ouvrez RStudio et créez un nouveau script R.
Enregistez le script dans le dossier “rna-seq” précédemment créé (à la racine de celui-ci), en utilisant le nom que vous voulez (e.g. rna-seq-frd3.R)
Lors de l’installation des packages, ceux-ci peuvent être disponibles soit sur le CRAN (ressources habituelles), soit sur BioConductor (qui regroupe un ensemble de softwares destinés à la bioinformatique) :
Installez les packages CRAN suivants :
# Install CRAN packages
install.packages("BiocManager") # Celui-ci nous permettra ensuite d'installer les packages Bioconductor
install.packages("VennDiagram")
install.packages("RColorBrewer")
Installez maintenant les packages disponibles sur Bioconductor :
# Install packages
BiocManager::install("openssl")
BiocManager::install("RCurl")
BiocManager::install("limma")
BiocManager::install("gplots")
BiocManager::install("EnhancedVolcano")
BiocManager::install("genefilter")
BiocManager::install("DESeq2")
Nous allons découvrir ci-dessous comment écrire un script R qui permet d’identifier des gènes différentiellement exprimés sur base d’une matrice de comptage. Vous devrez ajouter chaque nouvelle étape de manière séquentielle dans le même script R au sein de RStudio.
Nous allons d’abord chercher à identifier les gènes qui montrent une expression différentielle en fonction du génotype (WT vs mutant) indépendamment du traitement, puis ceux qui sont modulés par le traitement (Ctrl vs Treat) indépendamment du génotype, et enfin examiner l’effet du traitement sur chaque génotype. Chacune de ces analyses correspond à un design différent du test statistique, qui devra être pris en compte dans nos scripts.
Il est important de :
Step 1. Chargez les librairies et le jeu de donnée
# Create a new R Script
#
#Set-up the working directory
setwd("/Users/fredericbouche/Desktop/rna-seq") # Replace by your own file path
# This folder shold countain, as described above :
# - 1 folder called "data" with the table of raw count
# - 1 empty folder called "qc" to save the quality control data
# - 1 empty folder called "res_wt_cont_vs_treat" to save data from WT in cont vs treat
# - 1 empty folder called "res_wt_vs_mut" to save the coomparison "WT vs mutants"
# - 1 empty folder called "venn" to save gene list comparison data
#####
###### Library loading
library("DESeq2")
library("limma")
library("gplots")
library("RColorBrewer")
library("genefilter")
library("VennDiagram")
library("EnhancedVolcano")
#####
# The matrix is stored in a "data" folder located inside the wd setup above.
# This matrix countains the raw count of all the gene in all the samples.
countData=read.table("data/tophat_root.matrix",header=TRUE,row.names=1,sep="\t")
# Let's look at the file :
head(countData)
Notes:
- N'oubliez pas de mettre à jour le chemin d'accès des fichiers.
- Observez ce qu'il se passe.
Step 2. Décrivons le jeu de donnée : nous devons indiquer à quoi correspondent les titres des colonnes (i.e. Quel génotype ? Quel traitement ?)
# We need to described what are the different header (e.g. wt_Ctrl.1, etc)
#Describe the dataset for each variable (genot = genotype ; treat = treatment ; g_t = genotype x treatment)
genot=as.factor(rep(c("WT","mut"),each=6))
treat=as.factor((rep(rep(c("Ctrl","Treat"),each=3),2)))
g_t=as.factor(rep(c("WT_Ctrl", "WT_Treat", "mut_Ctrl", "mut_Treat"),each=3))
#Create a new table to store the dataset description in a data frame
colData=data.frame(g_t,genot,treat,row.names=names(countData))
# Let's look at this table :
colData
Step 3. Construisons un modèle afin d’analyser l’effet du génotype. Pour cela nous allons utiliser la fonction DESeqDataSetFromMatrix, qui utilise une matrice de comptage associée à une table de description (créée juste avant) afin de générer un objet DESEQ :
#Build model using our dataset and the sample information
gt=DESeqDataSetFromMatrix(countData = countData,colData = colData, design = ~ g_t) # ici, le design g_t indique que cette cette colonne là qu'il faut utiliser afin de déterminer quels sont les échantillons à comparer entre eux.
# The, let's compute models, dispersion analysis and make the fitting to the model
gt_DESeq<- DESeq(gt) # Faites attention aux différents informations que R vous fait parvenir quand aux manipulations qui sont réalisées par cette fonction.
Step 4. QC : Déterminons le nombre de gènes qui ne sont pas détectés du tout dans l’analyse :
# Commencons par regarder la dimension de la table de comptage et stockons cela dans un objet (dim_init)
dim(countData)
dim_init <-dim(countData)
# Filtrons pour retirer les gènes qui ne sont pas détectés, et comptons les :
keep = rowSums(countData) > 0
countData_filtered <- countData[keep,]
dim(countData_filtered)
# Faisons maintenant la différence entre les deux :
dim(countData)[1]-(dim(countData_filtered)[1])
Combien de gènes n’ont-ils pas été détectés dans l’analyse ? Cela vous semble-t-il plausible ?
Step 5. QC : Afin de vérifier que les réplicas techniques sont effectivement proches les uns des autres, générons une heatmap indiquant les distances entre les échantillons :
### To carry out the quality control, we firt need to log transform the data
rld<-rlog(gt_DESeq)
### Now, let's build a heatmap of sample distances
#Build sample distance
sampleDist <- dist(t(assay(rld)))
#Build heatmap
sampleDistMatrix<-as.matrix(sampleDist)
# Déterminons le nom des colomnes et des rangées de la heatmap
rownames(sampleDistMatrix)<-paste(rld$g_t)
colnames(sampleDistMatrix)<-NULL
rownames(sampleDistMatrix) <- colnames(sampleDistMatrix) <- rownames(colData)
# Créons une palette de couleur
colours=colorRampPalette(rev(brewer.pal(9, "Blues")))(300)
# Enfin, générons la heatmap
png(filename = "qc/heatmap_sampledist_Treat_root.png", width = 1500,
height = 1500, units = "px", res = 150)
heatmap.2(sampleDistMatrix, dendrogram = "both", trace = "none", col = colours,
main = "Treat Root Sample Distance", margin=c(6, 8))
dev.off()
Step 6. Construisons un graphique de type PCA afin de mettre en évidence la proximité entre les différents réplicas de chaque génotype et traitement :
png("qc/pca_1.png", width=800, height=500, pointsize = 60)
plotPCA(rld,intgroup=c("g_t"),ntop=1000) # "intgoup" indique quels groupes comparer / "ntop" indique que la PCA utilise les 1000 gènes présentant le plus de variance dans le jeu de données
dev.off()
Step 7. Commencons par identifier les DEGs, pour le WT, entre la conditions standard (control) et l’excès de Zn (treatment):
# Comparons le WT en condition traitées vs contrôle
WT_treat_vs_cont=results(gt_DESeq, contrast = c("g_t", "WT_Treat", "WT_Ctrl"), lfcThreshold = 1, alpha = 0.05)
# Regardons l'objet obtenu :
head(WT_treat_vs_cont,5)
# La fonction summary permet de vérifier le résultat de la comparaison :
summary(WT_treat_vs_cont)
# Mettons tout cela sous la forme d'une table que l'on sauvegardera dans nos résultats (ici, tous les résultats sont sauvgardés, que les gènes soient DEG ou non) :
write.table(WT_treat_vs_cont,"res_wt_cont_vs_treat/pairwise_root_WT_Treat_vs_Ctrl.txt",sep="\t")
# Nous pouvons aussi créer des listes séparées pour les gènes induits (up) et réprimés (down) avec un logFC>1 (= FC>2) et une p.value ajustée < 0.05 :
WT_treat_vs_cont_up<- WT_treat_vs_cont[which(WT_treat_vs_cont$log2FoldChange > 1 & WT_treat_vs_cont$padj<0.05),]
WT_treat_vs_cont_down<- WT_treat_vs_cont[which(WT_treat_vs_cont$log2FoldChange < -1 & WT_treat_vs_cont$padj<0.05),]
# Et les sauvegarder (attention a bien utiliser une nomenclature claire afin de ne pas entrainer de confusion dans le sens de régulation des gènes ):
write.table(WT_treat_vs_cont_up,"res_wt_cont_vs_treat/higher_root_WT_Treat_vs_Ctrl.txt",sep="\t")
write.table(WT_treat_vs_cont_down,"res_wt_cont_vs_treat/lower_root_WT_Treat_vs_Ctrl.txt",sep="\t")
Step 8. It is often useful to visualize the data using Volcano plots.
# Let's draw a volcano plot
png("res_wt_cont_vs_treat/Volcano_WT_treat_Vs_ctrl.png",width=750,height=750)
EnhancedVolcano(WT_treat_vs_cont,
xlim = c(-10, 10),
lab = rownames(WT_treat_vs_cont),
x = 'log2FoldChange',
y = 'pvalue')
dev.off()
Step 9. Même chose, mais pour la comparaison entre le WT et le mutant :
res_wt_vs_mut=results(gt_DESeq, contrast = c("g_t", "mut_Ctrl", "WT_Ctrl"), lfcThreshold = 1, alpha = 0.05)
summary(res_wt_vs_mut)
head(res_wt_vs_mut,2)
# Créons à nouveau une table pour tous les gènes (ici, tous les résultats sont sauvgardés, que les gènes soient DEG ou non) :
write.table(res_wt_vs_mut,"res_wt_vs_mut/pairwise_root_Ctrl_mut_vs_WT.txt",sep="\t")
# Filtrons cela pour ne conserver que les gènes statistiquement induits et réprimés :
res_mut_vs_wt_up<- res_wt_vs_mut[which(res_wt_vs_mut$log2FoldChange > 1 & res_wt_vs_mut$padj<0.05),]
res_mut_vs_wt_down<- res_wt_vs_mut[which(res_wt_vs_mut$log2FoldChange < -1 & res_wt_vs_mut$padj<0.05),]
# Et sauvegardons le sous la forme d'un fichier txt :
write.table(res_mut_vs_wt_up,"res_wt_vs_mut/higher_root_Ctrl_mut_vs_WT.txt",sep="\t")
write.table(res_mut_vs_wt_down,"res_wt_vs_mut/lower_root_Ctrl_mut_vs_WT.txt",sep="\t")
# Créons maintenant le volcano plot
png("res_wt_vs_mut/Volcano_WT_vs_mutant.png",width=750,height=750)
EnhancedVolcano(res_wt_vs_mut,
xlim = c(-13, 13),
lab = rownames(res_wt_vs_mut),
x = 'log2FoldChange',
y = 'pvalue')
dev.off()
Pour la suite des analyses, passez à l’onglet 3 : analyses downstream.
Une chaîne YouTube appelée StatQuest explique de manière simple de nombreux concepts statistiques, qu’ils soient simples ou complexes. Un index complet de leurs vidéos est disponible ici : Statquest Videos index. Voici quelques vidéos qui pourraient vous aider à comprendre les statistiques utilisées par le package DESEQ2 aux différentes étapes de l’analyse, ainsi que les différents types de normalisation des lectures (RPKM, FPKM, TPM) :
Library normalization :
Independent filetering :
RNA-seq - The Problem with Technical replicates :
RPKM, FPKM, and TPM :