Libraries utilisées

ATTENTION : certaines libraries ne sont pas installées par défaut dans R, il faut les télécharger. Pour cela, il est nécessaire de trouver sur internet le package auquel elle correspond et télécharger avec la commande : install.packages(‘nom_package’) sans oublier de mettre les guillemets.

library(tidyverse) #Manipulation de données
library(limma) #Calculs des GDE 
library(NOISeq) #Normalisation des counts
library(ggplot2) #Représentation des données
library(ggrepel) #Mise en page des labels ggplot
library(ggpubr) #Mise en page des ggplot pour publication 
library(plotly) #pour faire des ggplots interactifs 
library(heatmaply) #Pour faire des heatmaps interactifs
library(rsconnect)

Importer les données de RNAseq

Les fichiers FASTq ont été préalablement traitées et les raw count rassemblés dans un fichier csv que l’on va charger.

RNAseq <- read.csv(file='RNAseq.csv', sep=';') #On charge le fichier brut des raw count 
colnames(RNAseq)[colnames(RNAseq) == 'ï..gene_symbol'] <- c('Gene.Symbol') #on renome la première colonne

#on construit un tableau avec les infos des échantillons 
Info <- data.frame(Nom_ech = c('DC JK 1','DC JK 2','DC JK 3','DC JK LPS 1','DC JK LPS 2','DC JK LPS 3','DC C91 1','DC C91 2','DC C91 3','DC C91 LPS 1','DC C91 LPS 2','DC C91 LPS 3','DC C81 1','DC C81 2','DC C81 3','DC C81 LPS 1', 'DC C81 LPS 2','DC C81 LPS 3'), 
                   HTLV = c(F,F,F,F,F,F,T,T,T,T,T,T,T,T,T,T,T,T), 
                   LPS = c(F,F,F,T,T,T,F,F,F,T,T,T,F,F,F,T,T,T), 
                   Condition = colnames(RNAseq[,-1]), 
                   Doneur = as.factor(c('A','B','C','A','B','C','A','B','C','A','B','C','A','B','C','A','B','C')))

Normalisation des données

norm.data <- tmm(RNAseq[,-1]) #on enlève le nom des gènes et on normalise par la méthode TMM
boxplot(norm.data, outline = F) #on visualise si la normalisation est ok

Analyse en composantes principales

On ne sait pas quels sont les échantillons qui se ressemblent le plus. Pour avoir une première idées des groupes qui sont présents on calcul des indices de corrélation que l’on représente en heatmaps. On confirmera avec une ACP cette heatmap.

correlation <- cor(norm.data) #on calcul les indices de corrélation
colnames(correlation) <- Info$Nom_ech #on ajoute aux colonne le noms des échantillons
rownames(correlation) <- Info$Nom_ech #on ajoute aux lignes le nomdes échantillons 
heatmaply(correlation, key.title = 'Indice de corrélation') #calcul des indices de corrélations entre échantillons et on le représente sous forme de heatmap interactive

Il semblerait que les échantillons se regroupent en deux populations mais restent finalement très proches les unes des autres.

ACP <- prcomp(t(norm.data)) #on transpose le tableau des counts normalisés et on calcule les composantes principales (réduction de dimensions) 
table.ACP <- as.data.frame(ACP$x) %>% mutate(Echantillon = colnames(RNAseq[-1])) %>% cbind(Info) #on récupère le tableau des résultats de la réduction de dimensions et on ajoute le nom des échantillons dans une colonne puis on ajoute les informations sur les échantillons du tableau Info
ggplot(table.ACP, aes(x=PC1, y=PC2, col=LPS, label=Doneur))+ geom_point(size=3) + geom_text_repel() + theme_pubr() #représentation graphique de l'ACP

Cette ACP classe les échantillons selon les donneurs et le traitement au LPS. On souhaite voir l’effet de l’infection par HTLV. Pour cela on réalise une ACP cette fois si en séparant les groupes traités ou non au LPS.

#ACP sur le groupe non traité au LPS
ACP_NL <- prcomp(t(norm.data[,Info$Condition[Info$LPS==F]])) #on calcule une ACP sur les échantillons qui ne sont pas traités au LPS
table.ACP_NL <- data.frame(ACP_NL$x) #on récupère le tableau des résultats de l'ACP
table.ACP_NL <- cbind(table.ACP_NL, Info[Info$LPS==F,]) #on ajoute les infos des échantillons qui ne sont pas traités au LPS
A <- ggplot(table.ACP_NL, aes(x=PC1, y=PC2, col=HTLV, label=Nom_ech)) + geom_point(size=3) + geom_text_repel() + theme_pubr() #on stock le graph dans une variable

#ACP sur le gourpe traité au LPS
ACP_L <- prcomp(t(norm.data[,Info$Condition[Info$LPS==T]]))
table.ACP_L <- data.frame(ACP_L$x)
table.ACP_L <- cbind(table.ACP_L, Info[Info$LPS==T,])
B <- ggplot(table.ACP_L, aes(x=PC1, y=PC2, col=HTLV, label=Nom_ech)) + geom_point(size=3) + geom_text_repel() + theme_pubr()

ggarrange(A,B, labels=c('A','B'), common.legend=T)