ACP

DESeq2 usage with R
Club Bioinfo
Author

Hélène Polvech

Published

October 3, 2019

Introduction

Méthode statistique permettant d’explorer des données dites multivariées (données avec plusieurs variables).

Chaque variable pourrait être considérée comme une dimension différente. L’ACP synthétise cette information en seulement quelques nouvelles variables appelées composantes principales. L’information contenue dans un jeu de données correspond à la variance ou l’inertie totale qu’il contient.

L’objectif de l’ACP est d’identifier les directions (i.e., axes principaux ou composantes principales) le long desquelles la variation des données est maximale.

En résumé, l’analyse en composantes principales (PCA = ACP ) permet:

  • d’identifier des “profils cachés” dans un jeu de données,
  • de réduire les dimensions des données en enlevant la redondance des données,
  • d’identifier les variables corrélées

Source : http://www.sthda.com/french/articles/38-methodes-des-composantes-principales-dans-r-guide-pratique/73-acp-analyse-en-composantes-principales-avec-r-l-essentiel/

Vocabulaire

Dans le cadre de l’étude de l’expression différentielle avec une ACP, plusieurs termes sont utilisés :

  • Individus = Échantillons
  • Variables = Gènes
  • Dimensions = Composantes Principales ( ~ Axe )
  • Valeurs propres = Variances des composantes principales
  • Coordonnées = Coordonnées des variables pour créer un nuage de points.
  • Cos² = Cosinus carré des variables. Représente la qualité de représentation des variables sur le graphique de l’ACP.
  • Contribution = « Poids » apporté par la variable ou l’individu en pourcentage

Librairies R

      library("FactoMineR") # PCA
        library("factoextra") # PCA
        library("ggplot2") # Graphiques
        library(RColorBrewer) # Couleurs

Les données IN

Tous types de tableau de comptages (NGS, qPCR, Intensités, etc.)

      cts <- read.csv2("../data/ClubBioInfo_PCA_readscounts_normalise_19311genes_v2.csv", header=T,
        stringsAsFactors=F, sep=",", dec=",", row.names = "symbol")

        str(cts)
    
      ## 'data.frame': 19311 obs. of 21 variables:
        ## $ hESC_n1: num 707.53 29.48 441.14 10.72 6.43 ...
        ## $ hESC_n2: num 799.27 1.3 461.57 3.91 3.91 ...
        ## $ hESC_n3: num 1254.22 36.82 569 9.78 13.23 ...
        ## $ T2_n1 : num 897.81 0 251.45 25.79 8.87 ...
        ## $ T2_n2 : num 1233.087 0 107.826 13.824 0.922 ...
        ## $ T2_n3 : num 1029.51 0 149.14 13.56 2.71 ...
        ## $ T3_n1 : num 1249.83 19.43 227.46 6.48 0 ...
        ## $ T3_n2 : num 1216.8 0 158.7 15.4 0 ...
        ## $ T3_n3 : num 1210.59 1.46 173.36 18.94 0 ...
        ## $ Orgd_n1: num 675.25 2.35 253.66 17.62 2.35 ...
        ## $ Orgd_n2: num 743.55 2.63 274.56 27.59 1.31 ...
        ## $ Orgd_n3: num 787.7 0 220.38 15.82 4.52 ...
        ## $ T1_n3 : num 982.78 0 311.15 23.08 7.96 ...
        ## $ T1_n1 : num 1114.3 9.15 172.91 11.89 7.32 ...
        ## $ T1_n2 : num 1118.08 0 218.27 5.94 0 ...
        ## $ T4_n1 : num 1014.23 0 229.06 6.61 3.3 ...
        ## $ T4_n2 : num 936.5 0.9 214.1 0.9 0 ...
        ## $ T4_n3 : num 892.04 0 288.28 8.7 7.61 ...
        ## $ T5_n1 : num 804.28 11.05 266.25 7.73 7.73 ...
        ## $ T5_n2 : num 414.41 0 365.89 6.56 2.62 ...
        ## $ T5_n3 : num 413.74 0 313.22 11.65 2.91 ...

Quantité de variance expliquée par chaque axe principal : les valeurs propres

Identification & choix du nombre d’axes à étudier en fonction du pourcentage de variance.

      res.pca <- PCA(t(cts), graph=F, scale.unit = T )
        eig.val <- get_eigenvalue(res.pca)
  eigenvalue variance.percent cumulative.variance.percent
Dim.1 3771.8238 21.505353 21.50535
Dim.2 2520.6070 14.371441 35.87679
Dim.3 1638.2416 9.340564 45.21736
Dim.4 1418.1466 8.085675 53.30303
Dim.5 1010.0338 5.758788 59.06182
Dim.6 798.0465 4.550125 63.61195
Dim.7 740.2814 4.220773 67.83272
Dim.8 683.0535 3.894484 71.72720
Dim.9 622.9754 3.551944 75.27915
Dim.10 537.0363 3.061955 78.34110
Dim.11 476.9013 2.719091 81.06019
Dim.12 442.4067 2.522417 83.58261
Dim.13 411.9748 2.348907 85.93152
Dim.14 407.0720 2.320953 88.25247
Dim.15 397.6561 2.267268 90.51974
Dim.16 384.8314 2.194147 92.71388
Dim.17 345.0296 1.967213 94.68110
Dim.18 317.9959 1.813079 96.49418
Dim.19 313.1212 1.785285 98.27946
Dim.20 301.7651 1.720538 100.00000

Lorsque eigenvalue > 1 , indique que la composante principale représente plus de variance par rapport à une seule variable d’origine.

      fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 30))

Résultat des variables ( gènes )

      var <- get_pca_var(res.pca)
        var
    
      ## Principal Component Analysis Results for variables
        ## ===================================================
        ## Name Description
        ## 1 "$coord" "Coordinates for the variables"
        ## 2 "$cor" "Correlations between variables and dimensions"
        ## 3 "$cos2" "Cos2 for the variables"
        ## 4 "$contrib" "contributions of the variables"
  • var$coord: coordonnées des variables pour créer un nuage de points.
  • var$cos2: cosinus carré des variables. Représente la qualité de représentation des variables sur le graphique de l’ACP. Il est calculé comme étant les coordonnées au carré: var.cos2 = var.coord * var.coord .
  • var$contrib: contient les contributions (en pourcentage), des variables, aux composantes principales. La contribution d’une variable (var) à une composante principale donnée: (var.cos2 * 100) / (total cos2 du composant) .
      summary(var$contrib[,"Dim.1"])
    
      ## Min. 1st Qu. Median Mean 3rd Qu. Max.
        ## 0.0000000 0.0004117 0.0025756 0.0051784 0.0080659 0.0258162
    
      fviz_contrib(res.pca, choice = "var", axes = 1, top = 20)
    

      fviz_pca(res.pca, select.var = list(name = c("NANOG", "VSX2","BMP4")),
        col.ind = "cos2", axes = c(1, 2), gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")) +
        labs(title ="genes")
    

      fviz_pca_var(res.pca, col.var = "contrib", title="Cercle de corrélation",
        gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
    

      fviz_pca_var(res.pca, col.var = "cos2", title="Cercle de corrélation",
        gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

  • Contribution : Les variables corrélées avec PC1 (i.e., Dim.1) et PC2 (i.e., Dim.2) sont les plus importantes pour expliquer la variabilité dans le jeu de données.

  • Qualité : Un cos2 élevé indique une bonne représentation de la variable sur les axes principaux en considération. Dans ce cas, la variable est positionnée à proximité de la circonférence du cercle de corrélation.

Description des dimensions

      res.desc <- dimdesc(res.pca, axes = c(1:3), proba = 0.05)

        str(res.desc)
    
      ## List of 3
        ## $ Dim.1:List of 1
        ## ..$ quanti: num [1:8863, 1:2] NA NA NA NA NA NA NA NA NA NA ...
        ## .. ..- attr(*, "dimnames")=List of 2
        ## .. .. ..$ : chr [1:8863] NA NA NA NA ...
        ## .. .. ..$ : chr [1:2] "correlation" "p.value"
        ## $ Dim.2:List of 1
        ## ..$ quanti: num [1:6860, 1:2] NA NA NA NA NA NA NA NA NA NA ...
        ## .. ..- attr(*, "dimnames")=List of 2
        ## .. .. ..$ : chr [1:6860] NA NA NA NA ...
        ## .. .. ..$ : chr [1:2] "correlation" "p.value"
        ## $ Dim.3:List of 1
        ## ..$ quanti: num [1:4867, 1:2] NA NA NA NA NA NA NA NA NA NA ...
        ## .. ..- attr(*, "dimnames")=List of 2
        ## .. .. ..$ : chr [1:4867] NA NA NA NA ...
        ## .. .. ..$ : chr [1:2] "correlation" "p.value"
  correlation p.value
CTPS1 0.9867843 0
CACYBP 0.9838903 0
POLR3K 0.9824619 0
DDX18 0.9804360 0
XRCC5 0.9782874 0
CSE1L 0.9766512 0
  correlation p.value
ZBTB10 -0.9032808 0
LIX1L -0.9146378 0
RABGAP1 -0.9150203 0
DAZAP2 -0.9185193 0
TOM1L2 -0.9339979 0
YAF2 -0.9442530 0

Résultats des individus (échantillons)

      ind <- get_pca_ind(res.pca)
        ind
    
      ## Principal Component Analysis Results for individuals
        ## ===================================================
        ## Name Description
        ## 1 "$coord" "Coordinates for the individuals"
        ## 2 "$cos2" "Cos2 for the individuals"
        ## 3 "$contrib" "contributions of the individuals"
    
      ind.coord <- ind$coord

        fviz_pca_ind (res.pca, col.ind = "cos2",
        gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
        repel = TRUE # Évite le chevauchement de texte
        )
    

      fviz_pca_ind (res.pca, pointsize = "cos2", axes = c(2, 3),
        pointshape = 21, fill = "#E7B800",
        repel = TRUE # Évite le chevauchement de texte
        )
    

      # Contribution totale sur PC1 et PC2
        fviz_contrib(res.pca, choice = "ind", axes = 1:2)

Graphique détaillé avec ggplot2

      coldata <- as.data.frame(matrix(nc=1, nr=21))
        colnames(coldata) <- "condition"
        rownames(coldata) <- colnames(cts)

        sampleCondition <- as.factor(c("hESC", "hESC", "hESC", "T2", "T2", "T2", "T3", "T3", "T3", "Orgd", "Orgd",
        "Orgd", "T1", "T1", "T1", "T4", "T4", "T4", "T5", "T5", "T5"))

        coldata$condition <- sampleCondition
  condition
hESC_n1 hESC
hESC_n2 hESC
hESC_n3 hESC
T2_n1 T2
T2_n2 T2
T2_n3 T2
T3_n1 T3
T3_n2 T3
T3_n3 T3
Orgd_n1 Orgd
Orgd_n2 Orgd
Orgd_n3 Orgd
T1_n3 T1
T1_n1 T1
T1_n2 T1
T4_n1 T4
T4_n2 T4
T4_n3 T4
T5_n1 T5
T5_n2 T5
T5_n3 T5
  Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
hESC_n1 133.95915 13.04094 8.8076541 15.92419 14.48629
hESC_n2 149.20792 -1.87774 0.1076187 14.43236 27.57976
hESC_n3 146.31142 -27.34490 -7.0822876 38.75226 -19.06069
T2_n1 -39.32805 -29.15552 -19.0136773 29.36034 -36.37820
T2_n2 -40.26080 -44.10873 -29.9457092 17.75943 30.56527
T2_n3 -51.37032 -53.16016 -38.6935852 15.01031 45.09895 
      ind.coldata <- data.frame(coldata,ind.coord[,c("Dim.1","Dim.2")])
  condition Dim.1 Dim.2
hESC_n1 hESC 133.959145 13.040938
hESC_n2 hESC 149.207917 -1.877740
hESC_n3 hESC 146.311424 -27.344901
T2_n1 T2 -39.328055 -29.155519
T2_n2 T2 -40.260802 -44.108735
T2_n3 T2 -51.370323 -53.160157
T3_n1 T3 -33.298809 -34.339007
T3_n2 T3 -31.548471 -37.695157
T3_n3 T3 -44.152702 -59.312252
Orgd_n1 Orgd -12.418703 102.434365
Orgd_n2 Orgd -16.594854 95.537522
Orgd_n3 Orgd -20.886353 106.599277
T1_n3 T1 27.065045 19.434908
T1_n1 T1 -13.938422 -19.835694
T1_n2 T1 -2.880837 -17.921199
T4_n1 T4 -9.336367 -24.826319
T4_n2 T4 -11.213221 -40.806605
T4_n3 T4 -3.293285 -37.578867
T5_n1 T5 -25.449104 -3.628655
T5_n2 T5 -49.011041 50.407089
T5_n3 T5 -51.562182 44.136706 
      p <- ggplot(data=ind.coldata, aes(x=Dim.1, y=Dim.2, colour=condition))
        p <- p + geom_point(size=4)
        print(p)