Analyse RNAseq

Packages utilisés

library(DESeq2) #Calcul des GDE
library(NOISeq) #normalisation de counts
library(tidyverse) #manipulation des données (en toute légalité)
library(gprofiler2) #analyse des GO enrichis
library(org.Hs.eg.db) #base de donnée GO Homo Sapiens
library(clusterProfiler) #enrichment analysis
library(enrichplot) 

library(ggplot2) #représentation des données
library(ggpubr) #mise en page des représentations graphiques
library(plotly) #graphiques interactifs
library(ggsci) #palettes de couleur des journaux de publication
library(ggrepel) #mise en page des texte dans les graphs (pour pas qu'ils se montent dessus)
library(heatmaply) #heatmaply interactives

Importer fichier

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 
rownames(RNAseq) <- RNAseq$Gene

#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]),
                   Type_cell = c('JK','JK','JK','JK','JK','JK','C91PL','C91PL','C91PL','C91PL','C91PL','C91PL','C81PL','C81PL','C81PL','C81PL','C81PL','C81PL'),
                   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 <- as.matrix(tmm(RNAseq[,-1])) #on normalise suivant la méthode TMM (M.D. Robinson et al. 2010)
boxplot(norm.data, outline=F) #on vérifie si la normalisation est ok sur un boxplot 

Analyses en composantes principales

ACP de tous les échantillons

all_acp <- prcomp(t(norm.data)) # on réalise l'ACP
all_acp <- cbind(all_acp$x,Info) # on extrait les résultats que l'on ajoute aux informations des échantillons 
ggplot(all_acp, aes(x=PC1, y=PC2, col=`LPS`, label=`Nom_ech`)) + geom_point(size=3) + geom_text_repel() + theme_pubr() + scale_color_hue(labels = c("Non traité au LPS", "Traité au LPS")) + theme(legend.title = element_blank())#représentation graphique

ACP des échantillons traités ou non au LPS

#ACP sur les échantillons traités au LPS
lps_acp <- prcomp(t(norm.data[,Info$Condition[Info$LPS]])) # on réalise l'ACP sur les échantillons traités au LPS
lps_acp <- cbind(lps_acp$x,Info[Info$LPS,]) # on extrait les résultats que l'on ajoute aux informations des échantillons 
A <- ggplot(lps_acp, aes(x=PC1, y=PC2, col=`HTLV`, label=Nom_ech)) + geom_point(size=3) + geom_text_repel() + theme_pubr() + ggtitle('ACP des échantillons traités au LPS') + theme(plot.title = element_text(size = 12, hjust = 0.5)) +  scale_color_hue(labels = c("Non infecté par HTLV-1", "Infecté par HTLV-1")) + theme(legend.title = element_blank())  #représentation graphique

#ACP sur les échantillons non traités au LPS
nlps_acp <- prcomp(t(norm.data[,Info$Condition[Info$LPS==F]])) # on réalise l'ACP sur les échantillons traités au LPS
nlps_acp <- cbind(nlps_acp$x,Info[Info$LPS==F,]) # on extrait les résultats que l'on ajoute aux informations des échantillons 
B <- ggplot(nlps_acp, aes(x=PC1, y=PC2, col=`HTLV`, label=Nom_ech)) + geom_point(size=3) + geom_text_repel() + theme_pubr() + ggtitle('ACP des échantillons non traités au LPS') + theme(plot.title = element_text(size = 12, hjust = 0.5)) +  scale_color_hue(labels = c("Non infecté par HTLV-1", "Infecté par HTLV-1")) + theme(legend.title = element_blank())  #représentation graphique

#mise en page des graphs 
ggarrange(A,B, ncol=2, nrow=1, labels= c('A','B'), common.legend = T, font.label=list(size=18, color='black', face='bold', family=NULL), legend=c('bottom'))

Préparation des counts pour l’analyse DESeq2 des échantillons traités au LPS

colnames(RNAseq) <- c('Gene',Info$Nom_ech) #on renomme les colonnes
cts <- RNAseq[,-1] #on transforme en matrice les colonnes contenant les counts (on enlève les noms des gènes)
cts <- cts[,Info$Nom_ech[Info$LPS]] #on ne conserve que les échantillons traités au LPS
coldata <- Info[c(-1,-2,-3,-7,-8,-9,-13,-14,-15), c('HTLV', 'LPS')] #On récupère les informations sur l'état d'infection et des types cellulaires des échantillons conservés
rownames(coldata) <- colnames(cts) #on donne le même noms aux échantillons dans les deux tableaux 

Analyse DESeq2 : transformer les counts en un fichier utilisable par le package

dds <- DESeqDataSetFromMatrix(countData = round(cts), 
                              colData = coldata,
                              design = ~ HTLV) #transforme les data normalisées par TMM d'une matrice en objet DSE
dds$HTLV <- factor(dds$HTLV, levels = c("FALSE","TRUE")) #on indique quelles sont les conditions que l'on compare 

On peut éventuellement rajouter des filtres avant l’analyse des GDE

Analyse des GDE

res <- DESeq(dds) %>% results() #calcul des GDE et on stock dans une variable les résultats du calculs

Tableau des GDE en condition traité au LPS

p_med <- 0.05 #seuil de la p.value ajustée à paritr duquel on considère que le gène est DE
fc_med <- 2 #seuil du log2FC à paritr duquel on considère que le gène est DE
LPS_GDE <- as.data.frame(res) %>% #on fait un tableau des résultats obtenus 
                  drop_na(padj) #on enlève les lignes dont padj contient un NA
LPS_GDE <- filter(LPS_GDE, padj<p_med&abs(log2FoldChange)>fc_med) #on ne récupère les lignes des gènes qui sont dans les seuils définis plus hauts
LPS_GDE$Gene <- rownames(LPS_GDE)
print(paste(c('Nombre de GDE ='), dim(LPS_GDE)[1]))

## [1] "Nombre de GDE = 541"

print(paste(c('Nombre de gènes surexprimés ='), dim(LPS_GDE[LPS_GDE$log2FoldChange>0,])[1]))

## [1] "Nombre de gènes surexprimés = 440"

print(paste(c('Nombre de gènes sousexprimés ='), dim(LPS_GDE[LPS_GDE$log2FoldChange<0,])[1]))

## [1] "Nombre de gènes sousexprimés = 101"

Représentation en Volcano plot

volcano <- as.data.frame(res) %>% #on prends les résultats du calcul de DEG et on en fait un tableau 
              mutate(significatif = padj<p_med&abs(log2FoldChange)>fc_med) %>% # on ajoute une colonne pour savoir quels sont les gènes qui sont dans les seuils définis 
                mutate(Label=RNAseq$Gene) %>% #on rajoute le nom des gènes dans une colonne (plus simple que d'utiliser les rownames)
                    drop_na(padj) #on enlève les lignes dont la colonne padj contient des NA

volcano$GDE <- ifelse(volcano$significatif,c('Gènes différentiellement exprimés \n p.val.adjd < 0.05 - Log2FC > 2'),c('Gènes non différentiellement exprimés')) #on associe une légende à chaque gène différentiellement exprimé

volcaplot <- ggplot(volcano, aes(x=log2FoldChange, y=-log10(pvalue), col=GDE, label = Label)) + 
  geom_point() +
  theme_pubr() + 
  theme(legend.title = element_blank()) +
  scale_color_manual(values = c("#CF4E9CFF", "#2E2A2BFF"))#on en fait un graph

ggplotly(volcaplot)

GO Enrichment analysis des gènes surexprimés_LPS infecté vs. non infecté

Analyse des GO enrichis

up <- subset(LPS_GDE, log2FoldChange > 0) #on sélectionne les gènes différentiellement surexprimés
#calculs des GO enrichis 
ego <- enrichGO(gene = rownames(up), #liste de gènes dont on souhaite connaitre quels sont les go enrichis 
                universe = rownames(res), #liste de gènes de référence
                OrgDb = org.Hs.eg.db, #base de donnée de GO de l'espèce humaine
                ont           = "BP", #on sélectionne une partie des GO : ceux de biological process 
                pAdjustMethod = "BH", #méthode de calculs de pvalue ajustée 
                pvalueCutoff= 0.01,
                qvalueCutoff= 0.05,
                keyType = 'SYMBOL', #on lui donne en entrée des symbols de gènes 
        readable= TRUE)
dotplot(ego, showCategory = 5) + theme_pubr() + theme(legend.position = 'right') + geom_text_repel(data=ego, aes(label=ID), max.overlaps=Inf) #on réalise un plot des 10 GO les plus enrichis

heatplot(ego, showCategory = 5)

cnetplot(ego, node_label='category', cex_label_gene=1.2)

## Warning: Removed 5 rows containing missing values (geom_text_repel).

Enrichment analysis des gènes surexprimés

gp_up <- gost(row.names(up), organism = "hsapiens", custom_bg =rownames(res))

## Detected custom background input, domain scope is set to 'custom'

gp_up_link <- gost(row.names(up), organism = "hsapiens", as_short_link =  T)
resultgpup <- gp_up$result %>% arrange(desc(term_size))
gostplot(gp_up, interactive = T, capped=F)

```