{"id":1709,"date":"2015-11-16T14:03:22","date_gmt":"2015-11-16T19:03:22","guid":{"rendered":"http:\/\/bioinfo.iric.ca\/fr\/?p=1709"},"modified":"2017-04-29T23:05:37","modified_gmt":"2017-04-30T03:05:37","slug":"effectuer-une-pca-sur-les-donnees-de-leucegene","status":"publish","type":"post","link":"https:\/\/bioinfo.iric.ca\/fr\/effectuer-une-pca-sur-les-donnees-de-leucegene\/","title":{"rendered":"Effectuer une PCA sur les donn\u00e9es de Leuc\u00e9g\u00e8ne"},"content":{"rendered":"

<\/b>
\nGEO est une source extr\u00eamement riche de donn\u00e9es de profils transcriptionnels, mais t\u00e9l\u00e9charger et pr\u00e9parer ces donn\u00e9es constituent bien souvent un obstacle pour les apprentis bioinformaticiens. La d\u00e9monstration qui suit devrait faciliter vos premiers pas, j’utiliserai le jeu de donn\u00e9es de Leuc\u00e9g\u00e8ne. Une fois ces donn\u00e9es charg\u00e9es et pr\u00eates \u00e0 \u00eatre utilis\u00e9es dans R, je pr\u00e9senterai une perspective simplifi\u00e9e mais\u00a0pratique de l’utilisation de la PCA (Principal Component Analysis) pour faire de l’analyse exploratoire d’un ensemble de profiles transcriptionnels.<\/p>\n

Chargement des\u00a0donn\u00e9es<\/b><\/h2>\n

Un jeu de donn\u00e9es de 285 profils d’expression de leuc\u00e9mie my\u00e9lo\u00efde aigu\u00eb (LMA) peut \u00eatre t\u00e9l\u00e9charg\u00e9 de GEO \u00e0 cette adresse<\/a>\u00a0(Accession: GSE67040). Cette cohorte a \u00e9t\u00e9 s\u00e9quenc\u00e9e (RNA-Seq) dans le cadre du projet Leuc\u00e9g\u00e8ne (site web<\/a>) et repr\u00e9sente 285 patients. Une fois t\u00e9l\u00e9charg\u00e9 de GEO, le fichier \u00ab\u00a0GSE67040_RAW.tar\u00a0\u00bb peut \u00eatre \u00ab\u00a0d\u00e9tar\u00e9\u00a0\u00bb sur Mac en double cliquant dessus (sur linux, vous pouvez utiliser la commande \u00ab\u00a0tar xvf GSE67040_RAW.tar\u00a0\u00bb). Cela g\u00e9n\u00e9rera\u00a0285 fichiers .txt, chacun contenant l’expression de 24 553 g\u00e8nes. Les fichiers ressembleront \u00e0 ceci:<\/p>\n

chr1 11874 14361 DDX11L1 0.04970 1400\r\nchr1 14410 29370 WASH7P 1.95957 59227\r\nchr1 69091 70008 OR4F5 0.00620 100\r\nchr1 134773 140566 LOC729737 2.24401 215729\r\nchr1 567705 567793 MIR6723 0.00000 0<\/em><\/span><\/pre>\n
Pour charger un seul fichier en R, vous pouvez\u00a0utiliser:<\/p>\n
\r\ndata <- read.delim (\"GSM1636929_14H038_RPKM.txt\", header=F, col.names=c(\"chr\", \"start\", \"end\", \"sym\", \"rpkm\", \"count\"))\r\nhead (data)\r\n<\/span><\/code><\/pre>\n

\nchr start end sym rpkm count<\/code><\/em>
\n 1 chr1 11874 14361 DDX11L1 0.04970 1400<\/code><\/em>
\n 2 chr1 14410 29370 WASH7P 1.95957 59227<\/code><\/em>
\n 3 chr1 69091 70008 OR4F5 0.00620 100<\/code><\/em>
\n 4 chr1 134773 140566 LOC729737 2.24401 215729<\/code><\/em>
\n<\/span><\/p>\n
Pour charger tous les fichiers, nous devons boucler \u00e0 travers les noms de fichier\u00a0contenant\u00a0un motif donn\u00e9 (\u00ab\u00a0GSM.*_RPKM\u00a0\u00bb) et concat\u00e9ner les r\u00e9sultats dans un data frame<\/em>. Il ne faut pas oublier d’ajouter une colonne pour identifier l’\u00e9chantillon pour lequel\u00a0les valeurs d’expression ont \u00e9t\u00e9 mesur\u00e9es. Les quelques lignes de code suivantes prendront 5 \u00e0 10 minutes \u00e0 s’ex\u00e9cuter (la fonction \u00ab\u00a0rbind\u00a0\u00bb est plut\u00f4t lente. Le code pourrait \u00eatre optimis\u00e9, mais je laisse cette t\u00e2che aux lecteurs plus curieux!) :<\/p>\n
all.fn <- list.files (pattern=\"GSM.*_RPKM.txt\")\r\nall.data <- data.frame()\r\ndata.len <- sapply (all.fn, function (fn) {\r\n    pat <- unlist (strsplit (fn, \"_\"))[2] ## Extraction de l'identificateur patient \r\n    data <- read.delim (fn, header=F, col.names=c(\"chr\", \"start\", \"end\", \"sym\", \"rpkm\", \"count\"))\r\n    data$pat <- pat                       ## Ajout d'une colonne identifiant le patient\r\n    all.data <<- rbind (all.data, data)   ## Ajout des diff\u00e9rents profils dans all.data\r\n    nrow (data)                           ## Pour v\u00e9rifier que les profils contiennent le m\u00eame nombre de g\u00e8nes, voir data.len\r\n}<\/span><\/code><\/pre>\n
Les donn\u00e9es dans all.data sont pr\u00e9sent\u00e9es dans un format que l’on appelle \u00ab\u00a0long\u00a0\u00bb: chaque rang\u00e9e repr\u00e9sente une seule valeur d’expression et les premi\u00e8res colonnes servent \u00e0 identifier le g\u00e8ne. Ici, nous combinerons les colonnes #1-4 et utiliserons le r\u00e9sultat comme nom de g\u00e8ne, nous conserverons la colonne #5 pour l’expression (mesur\u00e9e en RPKM). Pour les utilisateurs de l’excellent package de graphiques ggplot2, vous aurez reconnu que ce format long correspond \u00e0 la forme pr\u00e9f\u00e9r\u00e9e par le package. Ce format est aussi tr\u00e8s utile lorsque l’on travaille avec les fonctions \u00ab\u00a0lm\u00a0\u00bb et \u00ab\u00a0glm\u00a0\u00bb! Toutefois, nous pr\u00e9f\u00e9rons ici\u00a0travailler avec un format de matrice plus traditionnel o\u00f9 nous avons les g\u00e8nes en rang\u00e9es et les patients en colonnes. Je convertirai donc mes donn\u00e9es du format long au format \u00ab\u00a0large\u00a0\u00bb (wide<\/em>). Cela se fait facilement avec la fonction\u00a0\u00ab\u00a0acast\u00a0\u00bb du package<\/em> \u00ab\u00a0reshape2\u00a0\u00bb (le package reshape offre des fonctions pour effectuer des conversions des formats long \u00e0 large et de large \u00e0 long). Finalement, je transformerai les RPKM en utilisant une transformation \u00ab\u00a0started-log\u00a0\u00bb (introduite par J. Tuckey en 1957<\/a>!).<\/p>\n
library (reshape2)\r\nz <- acast (all.data, chr + start + sym ~ pat, value.var=\"rpkm\")\r\nzl <- log10 (z + 0.001)<\/span><\/code><\/pre>\n
Analyse en composantes principales (PCA)<\/b><\/p>\n
Les donn\u00e9es sont charg\u00e9es et sont pr\u00eates \u00e0 \u00eatre analys\u00e9es… mais par o\u00f9 commencer? C’est une bonne habitude de prendre un peu de temps pour explorer les donn\u00e9es. Un bon outil \u00e0 appliquer aux donn\u00e9es pour cela est la PCA. Il existe un grand nombre de r\u00e9f\u00e9rences expliquant les bases\u00a0math\u00e9matiques de cette m\u00e9thode (pour les lecteurs curieux, je vous r\u00e9f\u00e8re \u00e0 la pr\u00e9sentation exhaustive de Jolliffe 1986<\/a>). Pour ma part, j’aime penser \u00e0 la PCA comme une m\u00e9thode pour d\u00e9finir une s\u00e9rie de signatures de g\u00e8nes dans lesquelles \u00e0 chaque g\u00e8ne correspond un poids lui conf\u00e9rant une importance sp\u00e9cifique dans une signature donn\u00e9e.<\/p>\n
La premi\u00e8re signature (appel\u00e9e premi\u00e8re composante dans le jargon de la PCA ou PC1 dans R) identifie les g\u00e8nes qui peuvent le mieux diff\u00e9rencier les \u00e9chantillons du jeu de donn\u00e9es. Si l’on prend chaque profil transcriptionnel et qu’on multiplie l’expression des g\u00e8nes par les poids donn\u00e9s par la premi\u00e8re composante et qu’on\u00a0fait\u00a0la somme du r\u00e9sultat pour obtenir un score unique, ce score repr\u00e9sentera le maximum de\u00a0variance qu’il est possible de repr\u00e9senter par une seule valeur. Le seconde composante (ou signature) sera orthogonale \u00e0 la premi\u00e8re, c’est-\u00e0-dire qu’elle capturera des variations\u00a0qui sont\u00a0ind\u00e9pendantes de celles captur\u00e9es par la premi\u00e8re composante. Et ainsi de suite pour les autres composantes.<\/p>\n
En utilisant la PCA, nous remplacerons ainsi chaque profil transcriptionnel (24 553 expressions par \u00e9chantillon) par un ensemble de valeurs (les scores dont il \u00e9tait question plus haut) qui pr\u00e9servent les variations observ\u00e9es dans le jeu de donn\u00e9es complet.<\/p>\n
pr <- prcomp (t (zl), center=T)\r\nx <- predict (pr)<\/span><\/code><\/pre>\n
La variable \u00ab\u00a0pr$rotation\u00a0\u00bb contient les signatures (composantes), une par colonne, et la fonction \u00ab\u00a0predict\u00a0\u00bb est utilis\u00e9e pour calculer les \u00ab\u00a0scores\u00a0\u00bb de ces signatures pour chacun de nos \u00e9chantillons.<\/p>\n
plot (x[,c(1,2)], pch=20)\r\ntext (x[,c(1,2)], rownames (x))<\/span><\/code><\/pre>\n
Avec ce type de repr\u00e9sentation, il est tr\u00e8s facile d’identifier des valeurs extr\u00eames (outliers<\/em>) dans un jeu de donn\u00e9e de grande dimensionalit\u00e9. Si des \u00e9chantillons forment des groupes facilement identifiables, ils deviendraient aussi apparents. Le code suivant illustre comment la PCA peut \u00eatre utile comme premi\u00e8re \u00e9tape exploratoire. Il est possible d’identifier les \u00e9chantillons de leuc\u00e9mie appartenant aux classes FAB<\/a> M0-1 (M0 et M1 combin\u00e9es), M5 et M6. Ces types de leuc\u00e9mie sont reconnus pour leur profil d’expression caract\u00e9ristique. La sortie de la PCA montre clairement que la PC1 (la \u00ab\u00a0signature\u00a0\u00bb principale identifi\u00e9e par la PCA) permet d’identifier les M0-1 et que la PC2 permet de distinguer les M5 des M6.<\/p>\n
p.M6 <- scan (what=character())\r\n04H017 04H039 05H097 07H024 10H066 11H183 07H034 03H028 12H067 12H117\r\np.M5 <- scan (what=character ())\r\n03H041 07H089 07H134 08H089 10H026 12H096 14H031 01H001 02H026 02H032 02H033 03H016 03H067 04H080 04H111 04H115 04H118 04H121 05H025 05H065 05H128 05H143 06H010 06H016 06H073 06H074 07H041 07H045 07H063 08H011 08H021 08H063 08H108 08H139 09H008 09H010 09H024 09H032 09H098 09H102 10H022 10H055 10H058 10H127 11H006 11H083 11H095 11H192 13H019 13H141 13H185 03H052 04H041 05H013 05H094 05H181 07H055 07H155 08H004 08H049 08H104 09H042 09H106 09H111 10H029 10H031 10H072 10H161 11H126\r\np.M01 <- scan (what=character ())\r\n02H003 02H017 03H030 04H096 05H179 06H133 06H135 07H038 08H118 09H018 09H079 10H014 10H038 11H002 11H014 11H046 11H138 11H170 12H030 12H149 13H039 13H048 13H056 13H083 14H012 14H019 14H027 02H009 02H025 02H043 02H053 02H060 02H066 03H024 03H033 03H036 03H081 03H090 03H094 03H116 03H119 04H001 04H006 04H024 04H025 04H048 04H055 04H103 04H108 04H112 04H117 04H120 04H127 04H133 04H135 04H138 05H008 05H022 05H030 05H034 05H078 05H111 05H163 05H184 05H195 06H028 06H029 06H088 06H089 06H103 06H144 06H151 07H060 07H062 07H107 07H112 07H125 07H133 07H135 07H137 07H151 07H152 07H158 07H160 08H012 08H034 08H048 08H053 08H056 08H065 08H082 08H085 08H113 08H129 09H002 09H022 09H026 09H031 09H040 09H043 09H045 09H048 09H070 09H078 09H083 09H088 09H090 09H113 09H115 10H040 10H048 10H056 10H068 10H070 10H092 10H095 10H107 10H109 10H113 10H115 11H010 11H019 11H035 11H058 11H097 11H103 11H107 11H129 11H135 11H142 11H151 11H186 11H232 11H240 12H010 12H021 12H039 12H056 12H079 12H091 12H116 12H171 12H172 12H173 12H175 12H176 13H073 13H080 13H139 13H150 13H158 13H169 14H001 14H017 14H020\r\nplot (x[,c(1,2)], pch=20)\r\npoints (x[(rownames (x) %in% p.M6), c(1,2)], pch=20, col=\"yellow\")\r\npoints (x[(rownames (x) %in% p.M5), c(1,2)], pch=20, col=\"orange\")\r\npoints (x[(rownames (x) %in% p.M01), c(1,2)], pch=20, col=\"blue\")<\/span><\/code><\/pre>\n
\"final_pca\"<\/a>
Figure 1<\/strong>: Deux composantes principales de la PCA appliqu\u00e9es \u00e0 285 \u00e9chantillons de patients de Leuc\u00e9g\u00e8ne. \u00a0Les couleurs correspondent \u00e0 la classe FAB; les M0 et M1 sont repr\u00e9sent\u00e9s en bleu, les M5 en orange et les M6 en jaune.<\/p><\/div>\n","protected":false},"excerpt":{"rendered":"
GEO est une source extr\u00eamement riche de donn\u00e9es de profils transcriptionnels, mais t\u00e9l\u00e9charger et pr\u00e9parer ces donn\u00e9es constituent bien souvent un obstacle pour les apprentis bioinformaticiens. La d\u00e9monstration qui suit devrait faciliter vos premiers pas, j’utiliserai le jeu de donn\u00e9es de Leuc\u00e9g\u00e8ne. Une fois ces donn\u00e9es charg\u00e9es et pr\u00eates \u00e0 \u00eatre utilis\u00e9es dans R, je pr\u00e9senterai une perspective simplifi\u00e9e mais\u00a0pratique de l’utilisation de la PCA (Principal Component Analysis) pour faire de l’analyse exploratoire d’un ensemble de profiles transcriptionnels. Chargement des\u00a0donn\u00e9es  […]<\/a><\/p>\n","protected":false},"author":5,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[69,24],"tags":[127,150],"class_list":["post-1709","post","type-post","status-publish","format-standard","hentry","category-analyse-de-donnees","category-langage-r","tag-expression-de-genes","tag-pca"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/posts\/1709","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/users\/5"}],"replies":[{"embeddable":true,"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/comments?post=1709"}],"version-history":[{"count":15,"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/posts\/1709\/revisions"}],"predecessor-version":[{"id":3299,"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/posts\/1709\/revisions\/3299"}],"wp:attachment":[{"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/media?parent=1709"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/categories?post=1709"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/bioinfo.iric.ca\/fr\/wp-json\/wp\/v2\/tags?post=1709"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}