• Acquérir des compétences d’analyse de données génomiques
  • Méthodes (normalisation, tests à l’effet du génome, contrôle des faux positifs)
  • Démarches (comment choisir une méthode d’analyse ?)
• Se familiariser avec R et Rstudio
  • Mettre en oeuvre (limma, DESeq2)
  • Editer des rapports (Rmarkdown)

Modalités pédagogiques

• Apprentissage par mise en situation (tutoriel)
• Programme modulable selon questions

: environnement pour l’analyse de données et la visualisation

• R est un logiciel libre
• Site web : https://www.r-project.org
• Nov. 2021 : > 18 000 packages https://cran.r-project.org/
• Communautés R :
  • https://www.r-bloggers.com/,
  • https://rladies.org/,
  • https://r-toulouse.netlify.app/
  • https://stackoverflow.com/questions/
  • https://r-graph-gallery.com/

: interface pour R (IDE : Integrated Development Environment)

• Site web : https://www.rstudio.com/
• Environnement ergonomique pour l’analyse de données
• Outils pour l’édition de rapport d’études
  • Rmarkdown : http://rmarkdown.rstudio.com/
• Préparation de la session de travail
  • Créez un projet pour l’analyse des données RNA-seq

• 1 - Données RNA-seq
• 2 - Tests à l’échelle du génome
  • 2.1 - Normalisation des données
  • 2.2 - Modèles pour données de comptage
  • 2.3 - Transformation des données
• 3 - Sélection de gènes

Nombres de reads par gène et par échantillon:

rawdta <- read.table("DataTest_foie_15ind_4cat4ind_10742genes_countData.txt",
                      sep="\t",header=TRUE,stringsAsFactors=TRUE)
head(rawdta[,1:6],n=3) %>% kbl(caption="Extrait des données") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Gènes :

nbgenes <- nrow(rawdta) ; nbgenes

[1] 10742

Echantillons :

nbsamples <- ncol(rawdta) ; nbsamples

[1] 15

sample_labs <- colnames(rawdta) ; head(sample_labs)

[1] "F10_G_B" "F11_M_H" "F13_G_H" "F14_M_H" "F15_G_B" "F16_M_B"

Deux lignées x deux régimes :

genotype <- factor(substring(sample_labs,first=5,last=5))
diet <- factor(substring(sample_labs,first=7,last=7))
covariates <- data.frame(Genotype=genotype,Diet=diet)
rownames(covariates) <- colnames(rawdta)
summary(covariates)

 Genotype Diet 
 G:8      B:7  
 M:7      H:8  

Objectif : identifier les gènes dont l’expression moyenne dépend des conditions expérimentales (lignée, régime, …)

Méthode : tests d’un même effet sur chaque gène

df <- data.frame(LPCAT3=as.matrix(rawdta)["LPCAT3",],Diet=covariates$Diet)
df %>% ggplot() + geom_boxplot(aes(x=Diet,y=LPCAT3),fill="orange") +
  ggtitle("Gène LPCAT3 : Répartition des données de comptage par régime") +
  xlab("Régime") + ylab("Nombre de reads")

Modèle : \(Y_{ij} = \mu + \alpha_{i} + \epsilon_{ij}\)

où

• \(Y_{ij}\) : nombre de reads pour un gène dans le \(j\)ème échantillon du groupe \(i\) (régime B ou H)
• \(\mu\) : expression moyenne dans le groupe 1 (régime B)
• \(\alpha_2\) (fold-change) : différence d’expressions moyennes entre groupes 1 et 2
• \(\epsilon_{ij}\) suit une loi normale d’écart-type \(\sigma\)

Test pour le gène LPCAT3

dta <- data.frame(y=as.matrix(rawdta)["LPCAT3",],Diet=covariates$Diet)
mod0 <- lm(y~1,data=dta)
mod1 <- lm(y~Diet,data=dta)
anova(mod0,mod1) %>% kbl(caption="Table d'analyse de la variance") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

• Signal biologique : Sum of Sq
• Bruit : RSS

Test pour tous les gènes

my_ftest <- function(y,groupe) {
   mod0 <- lm(y~1,data=data.frame(y=y,Groupe=groupe))
   mod1 <- lm(y~Groupe,data=data.frame(y=y,Groupe=groupe))
   anova(mod0,mod1)[2,"Pr(>F)"]
}
pval_Fisher <- apply(X=rawdta,MARGIN=1,FUN=my_ftest,groupe=covariates$Diet)
pval_Fisher["LPCAT3"]

 LPCAT3 
0.00013 

data.frame(Pval=pval_Fisher) %>% ggplot() + aes(x = Pval) +
  geom_histogram(fill="orange",bins=15,col="white",breaks=seq(from=0,to=1,by=0.025)) +
  ggtitle("Répartition des p-values de l'effet régime - Test de Fisher") + 
  xlab("p-value") + ylab("Effectifs")

Combien de gènes positifs ?

positives_Fisher <- pval_Fisher<=0.05
select <- which(positives_Fisher)
length(select)

[1] 376

Combien de gènes positifs si l’expression moyenne de tous les gènes était la même dans les deux régimes ?

Sources d’hétérogénéité

• Variations technologiques :
• profondeur de séquençage (comparaison entre groupes d’échantillons),
• longueur des gènes (comparaison des expressions de différents gènes)

• Nature des données de comptage : \(\sigma^{2} = f(\mu_{i})\).

Nombre de reads pour un pseudo-échantillon de référence

prs <- exp(rowMeans(log(rawdta)))
head(cbind(rawdta[,1:6],PRS=prs),n=3) %>% 
  kbl(caption="Extrait des données avec PRS") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Ratios au pseudo-échantillon de référence

ratios_to_prs <- rawdta/tcrossprod(prs,rep(1,nbsamples))
head(ratios_to_prs[,1:6],n=3) %>% 
  kbl(caption="Extrait des ratios au PRS") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Facteurs de normalisation par échantillon

norm_factors <- apply(X=ratios_to_prs,MARGIN=2,FUN=median,na.rm=TRUE)
head(norm_factors) %>% t() %>%
  kbl(caption="Facteurs de normalisation") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Données de comptage normalisées

norm_dta <- rawdta/tcrossprod(rep(1,nbgenes),norm_factors)
head(norm_dta[,1:6],n=3) %>% 
  kbl(caption="Données de comptage normalisées") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Normalisation par le package DESeq2

dds <- DESeqDataSetFromMatrix(countData = rawdta,colData = covariates,design = ~ Diet)
dds <- estimateSizeFactors(dds)
norm_dta_deseq2 <- counts(dds, normalized=TRUE)
norm_dta_deseq2 <- round(norm_dta_deseq2)
head(norm_dta_deseq2[,1:6],n=3) %>% 
  kbl(caption="Données de comptage normalisées (DESeq2)") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Test pour tous les gènes

pval_Fisher <- apply(X=norm_dta_deseq2,MARGIN=1,FUN=my_ftest,groupe=covariates$Diet)
pval_Fisher["LPCAT3"]

  LPCAT3 
1.63e-07 

data.frame(Pval=pval_Fisher) %>% ggplot() + aes(x = Pval) +
  geom_histogram(fill="orange",bins=15,col="white",breaks=seq(from=0,to=1,by=0.025)) +
  ggtitle("Répartition des p-values de l'effet régime - Test de Fisher") + 
  xlab("p-value") + ylab("Effectifs")

Combien de gènes positifs ?

positives_Fisher <- pval_Fisher<=0.05
select <- which(positives_Fisher)
length(select)

[1] 1756

Rapport signal-bruit

my_wb_sd <- function(y,groupe) {
   mod0 <- lm(y~1,data=data.frame(y=y,Groupe=groupe))
   mod1 <- lm(y~Groupe,data=data.frame(y=y,Groupe=groupe))
   anova_tab <- anova(mod0,mod1)
   within_sd <- sqrt(anova_tab[2,"RSS"]/anova_tab[2,"Res.Df"])
   between_sd <- sqrt(anova_tab[2,"Sum of Sq"]/anova_tab[2,"Df"])
   return(list(Wsd=within_sd,Bsd=between_sd))
}
wb_sd <- apply(X=norm_dta_deseq2,MARGIN=1,FUN=my_wb_sd,groupe=covariates$Diet)
unlist(wb_sd["LPCAT3"])

LPCAT3.Wsd LPCAT3.Bsd 
       481       4854 

Rapport signal-bruit

w_sd <- unlist(lapply(wb_sd,function(x) x$Wsd))
b_sd <- unlist(lapply(wb_sd,function(x) x$Bsd))
data.frame(w_sd,b_sd) %>% ggplot() + xlab("Différences d'expressions moyennes (log)") + 
  ylab("Ecarts-types intra-régimes (log)") + 
  geom_point(aes(x=log(b_sd),y=log(w_sd)),col="orange") +
  ggtitle("Relation entre signal et bruit") + geom_abline(intercept=0,slope=1)

Modèle de Poisson d’analyse de la variance à un facteur : \(Y_{ij} \sim {\cal P} ( \lambda_{i} )\)

où, pour un gène,

• \(Y_{ij}\) est le nombre de reads dans le \(j\)ème échantillon du groupe \(i\)
• \(\log ( \lambda_{i} ) = \mu + \alpha_{i}\)
• \(\mu\) : log-expression moyenne dans le groupe 1 (régime B)
• \(\alpha_2\) (log-fold-change) : log-ratio d’expressions moyennes entre groupes 1 et 2

Propriété : \(Var(Y_{ij}) = \lambda_{i}\)

Relation variance-moyenne : \(Var(Y_{ij}) \geq \lambda_{i}\)

mu_B <- rowMeans(norm_dta_deseq2[,covariates$Diet=="B"])
v_B <- apply(X=norm_dta_deseq2[,covariates$Diet=="B"],MARGIN = 1,FUN=var)
data.frame(Mean=mu_B,Var=v_B) %>% ggplot() + 
  xlab("Expressions moyennes (log) - Régime B") + ylab("Variance intra-régimes (log)") + 
  geom_point(aes(x=log(Mean),y=log(Var)),col="orange") +
  ggtitle("Relation entre variance et moyenne") + geom_abline(intercept=0,slope=1)

Test de \(\chi^{2}\) pour le gène LPCAT3

dta <- data.frame(y=norm_dta_deseq2["LPCAT3",],Diet=covariates$Diet)
mod0 <- glm(y~1,data=dta,family="poisson")
mod1 <- glm(y~Diet,data=dta,family="poisson")
anova(mod0,mod1,test="Chisq") %>% kbl(caption="Table d'analyse de la déviance") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

• Signal biologique : Deviance
• Bruit : Resid. Dev

Test pour tous les gènes

my_chi2_test <- function(y,groupe) {
   dta <- data.frame(y=y,Groupe=groupe)
   mod0 <- glm(y~1,data=dta,family="poisson")
   mod1 <- glm(y~Groupe,data=dta,family="poisson")
   anova(mod0,mod1,test="Chisq")[2,"Pr(>Chi)"]
}
pval_Poisson <- apply(X=norm_dta_deseq2,MARGIN=1,FUN=my_chi2_test,groupe=covariates$Diet)
pval_Poisson["LPCAT3"]

LPCAT3 
     0 

data.frame(Pval=pval_Poisson) %>% ggplot() + aes(x = Pval) +
  geom_histogram(fill="orange",bins=15,col="white",breaks=seq(from=0,to=1,by=0.025)) +
  ggtitle("Répartition des p-values de l'effet régime - Modèle de Poisson") + 
  xlab("p-value") + ylab("Effectifs")

Combien de gènes positifs ?

positives_Poisson <- pval_Poisson<=0.05
select <- which(pval_Poisson<=0.05)
length(select)

[1] 7219

Concordance Fisher-Poisson

table(positives_Fisher,positives_Poisson)

                positives_Poisson
positives_Fisher FALSE TRUE
           FALSE  3519 5466
           TRUE      3 1753

Modèle Binomial Négatif d’analyse de la variance : \(Y_{ij} \sim {\cal NB} ( \mu_{i} , \theta )\)

où, pour un gène,

• \(Y_{ij}\) est le nombre de reads dans le \(j\)ème échantillon du groupe \(i\)
• \(\log ( \mu_{i} ) = \mu + \alpha_{i}\)
• \(\mu\) : log-expression moyenne dans le groupe 1 (régime B)
• \(\alpha_2\) (log-fold-change) : log-ratio d’expressions moyennes entre groupes 1 et 2
• \(\theta\) : paramètre de dispersion

Propriété : \(Var (Y_{ij}) = \mu_{i} + \frac{\mu_{i}^2}{\theta}\)

Test de \(\chi^{2}\) pour le gène LPCAT3

dta <- data.frame(y=norm_dta_deseq2["LPCAT3",],Diet=covariates$Diet)
mod0 <- glm.nb(y~1,data=dta)
mod1 <- glm.nb(y~Diet,data=dta)
anova(mod0,mod1,test="Chisq") %>% kbl(caption="Test du rapport de vraisemblance") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

• Signal biologique : LR stat.
• Bruit : -2 x log-lik.

Test pour tous les gènes

my_chi2_test <- function(y,groupe) {
   dta <- data.frame(y=y,Groupe=groupe)
   mod0 <- glm.nb(y~1,data=dta)
   mod1 <- glm.nb(y~Groupe,data=dta)
   anova(mod0,mod1,test="Chisq")[2,"Pr(Chi)"]
}
pval_nb <- apply(X=norm_dta_deseq2,MARGIN=1,FUN=my_chi2_test,groupe=covariates$Diet)
pval_nb["LPCAT3"]

  LPCAT3 
1.06e-08 

data.frame(Pval=pval_nb) %>% ggplot() + aes(x = Pval) +
  geom_histogram(fill="orange",bins=15,col="white",breaks=seq(from=0,to=1,by=0.025)) +
  ggtitle("Répartition des p-values de l'effet régime - Modèle Binomial Négatif") + 
  xlab("p-value") + ylab("Effectifs")

Combien de gènes positifs ?

positives_nb <- pval_nb <= 0.05
select <- which(positives_nb)
length(select)

[1] 2321

Concordance Fisher - Binomial Négatif

table(positives_Fisher,positives_nb)

                positives_nb
positives_Fisher FALSE TRUE
           FALSE  8415  571
           TRUE      6 1750

dds <- DESeq(dds)
res <- results(dds)
mcols(res)$description

[1] "mean of normalized counts for all samples" "log2 fold change (MLE): Diet H vs B"      
[3] "standard error: Diet H vs B"               "Wald statistic: Diet H vs B"              
[5] "Wald test p-value: Diet H vs B"            "BH adjusted p-values"                     

data.frame(Pval=res$pvalue) %>% ggplot() + aes(x = Pval) +
  geom_histogram(fill="orange",bins=15,col="white",breaks=seq(from=0,to=1,by=0.025)) +
  ggtitle("Répartition des p-values de l'effet régime - DESeq2") + xlab("p-value") + 
  ylab("Effectifs")

Combien de gènes positifs ?

pval_deseq2 <- res$pvalue
positives_deseq2 <- pval_deseq2 <= 0.05
select <- which(positives_deseq2)
length(select)

[1] 1995

Objectif : transformer les données pour pouvoir utiliser le modèle linéaire classique

Transformation logarithmique : \(Z_{ij} = \mbox{log}_{2} ( Y_{0} + Y_{ij} )\)

vsd <- vst(dds, blind=TRUE)
pseudo_counts <- assay(vsd)
head(pseudo_counts[,1:6],n=3) %>% kbl(caption="Table d'analyse de la déviance") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Modèle d’analyse de la variance à l’échelle du génome

design = model.matrix(~Diet,data=covariates)   
fit = lmFit(pseudo_counts,design)
head(fit$coefficients,n=3) %>% 
  kbl(caption="Coefficients estimés") %>%
  kable_paper(bootstrap_options = "striped", full_width = F)

Test de Fisher pour le gène k

\(F_{k}=\frac{\frac{\mbox{RSS}_{0k}-\mbox{RSS}_{k}}{\mbox{df}_{0k} - \mbox{df}_{k}}}{\frac{\mbox{RSS}_{k}}{\mbox{df}_{k}}}\)

Test de Fisher modéré pour le gène k

\(F^{*}_{k}=\frac{\frac{\mbox{RSS}_{0k}-\mbox{RSS}_{k}}{\mbox{df}_{0k} - \mbox{df}_{k}}}{(1-p_{k})s^{2}_{0} + p_{k} \frac{\mbox{RSS}_{k}}{\mbox{df}_{k}}}\)

Méthode d’estimation de \(p_{k}\) et \(s^{2}_{0}\) : empirical Bayes (fonction eBayes)

fit <- eBayes(fit)
pval_limma <- fit$p.value[,"DietH"]
pval_limma["LPCAT3"]

  LPCAT3 
2.29e-08 

Gènes positifs

positives_limma <- pval_limma <= 0.05
sum(positives_limma)

[1] 1786

data.frame(Limma=-log10(pval_limma),Deseq2=-log10(pval_deseq2)) %>% ggplot() + 
  xlab("P-values limma (-log)") + ylab("P-values DESeq2 (-log)") + 
  geom_point(aes(x=Limma,y=Deseq2),col="orange") +
  ggtitle("Relation entre p-values Limma et Deseq2") + geom_abline(intercept=0,slope=1)

Densité de répartition des p-values : \(f (p ) = \pi_{0} \mathbb{1}_{[0,1]} ( p ) + ( 1 - \pi_{0} ) f_{1} ( p )\) où \(\pi_{0}\) est la proportion de gènes non-impactés par l’effet testé.

data.frame(Pval=pval_limma) %>% ggplot() + aes(x = Pval,y=..density..) +
  geom_histogram(fill="orange",bins=15,col="white",breaks=seq(from=0,to=1,by=0.025)) +
  ggtitle("Répartition des p-values de l'effet régime - Test de Fisher") + 
  xlab("p-value") + ylab("Effectifs") + geom_hline(yintercept=0.7)

Estimation du nombre de gènes non-impactés par l’effet testé

m <- length(pval_limma)
pi0 <- pval.estimate.eta0(pval_limma,diagnostic.plot=FALSE)
m0 <- m*pi0 ; m0

[1] 7524

Procédure

• Etape 1 - Pour le \(k\)ème gène, calcul de la p-value \(p_{k}\)
• Etape 2 - Choix d’un seuil sur la p-value : si \(p_{k} \leq t\) alors le gène \(k\) est positif

Nombre de gènes positifs si \(t=0.05\)

positives_limma <- pval_limma <= 0.05
sum(positives_limma)

[1] 1786

Définition : un gène est faux positif si la procédure de test conclut que l’effet testé est significatif, alors que cet effet n’existe pas.

• Nombre de gènes positifs : \(\mbox{P}_{t}\)
• Nombre de gènes faux positifs : \(\mbox{V}_{t}\) (non-observable)
• Proportion de faux positifs : \(\mbox{FDP}_{t} = \frac{V_{t}}{P_{t}}\) (non-observable)

Taux de faux positifs : \(\mbox{FDR}_{t}=\mathbb{E} ( \mbox{FDP}_{t})\)

Comment choisir le seuil \(t\) pour garantir \(\mbox{FDR}_{t} \leq 0.05\) ?

Supposons que l’effet testé n’existe pas pour \(m_{0}\) gènes, parmi \(m\) gènes.

Pour chacun de ces \(m_{0}\) gènes, \(\mathbb{P} ( p_{k} \leq t) = t\)

Comme \(V_{t} = \# \left\{ k , p_{k} \leq t \right\}\), \(\mathbb{E} ( V_{t} ) = m_{0} t\).

Estimation du FDR : \(\widehat{FDR}_{t} = \frac{m_{0}t}{P_{t}}\)

Exemple avec \(t=0.05\)

P <- sum(pval_limma <= 0.05)
FDR <- m0*0.05/P ; FDR

[1] 0.211

Principe : le seuil \(t\) est le plus grand possible pour lequel \(\widehat{\mbox{FDR}}_{t} \leq 0.05\)

Sélection pour un contrôle du FDR au seuil de 0.05

fdr_limma <- pi0*p.adjust(pval_limma,method="BH")
sum(fdr_limma <= 0.05)

[1] 220

Seuil de décision

p_sort <- sort(pval_limma)
fdr <- m0*p_sort/(1:m)
seuil <- max(p_sort[fdr<=0.05]) ; seuil

[1] 0.00145

volcanoplot(fit,coef=2)
abline(h = -log10(seuil),col = "orange",lty =2,lwd = 2)

Gènes postifs et ayant un ratio d’expressions moyennes suffisamment grand

log_fc <- fit$coefficients[,"DietH"]
select <- (fdr_limma <= 0.05) & (abs(log_fc) > 1)
sum(select)

[1] 21

x = pseudo_counts[select,]
heatmap(x, Rowv = FALSE, col = rev(heat.colors(16)))

heatmap(cor(t(x)), symm = TRUE, distfun = function(c) as.dist(1 - abs(c)), 
        col = rev(heat.colors(16)),trace="none",dendrogram="col")