Une introduction complète à votre génome avec la pile SciPy

La bioinformatique est un domaine interdisciplinaire qui développe des méthodes et des outils logiciels pour analyser et comprendre les données biologiques.

Plus simplement, vous pouvez simplement le considérer comme une science des données pour la biologie.

Parmi les nombreux types de données biologiques, les données génomiques sont l'une des plus largement analysées. Surtout avec les progrès rapides des technologies de séquençage d'ADN (NGS) de nouvelle génération, le volume de données génomiques a augmenté de façon exponentielle. Selon Stephens, Zachary D et al., l'acquisition de données génomiques se fait à l'échelle de l'exaoctet par an.

Dans cet article, je présente un exemple d'analyse d'un fichier GFF3 pour le génome humain avec SciPy Stack. Generic Feature Format Version 3 (GFF3) est le format de fichier texte standard actuel pour le stockage des caractéristiques génomiques. En particulier, dans cet article, vous apprendrez à utiliser la pile SciPy pour répondre aux questions suivantes sur le génome humain :

1. Quelle part du génome est incomplète ?
2. Combien y a-t-il de gènes dans le génome ?
3. Combien de temps dure un gène typique?
4. À quoi ressemble la distribution des gènes parmi les chromosomes ?

Le dernier fichier GFF3 pour le génome humain peut être téléchargé ici. Le fichier README fourni dans ce répertoire fournit une brève description de ce format de données, et une spécification plus approfondie se trouve ici.

Nous utiliserons Pandas, un composant majeur de la pile SciPy fournissant des structures de données rapides, flexibles et expressives, pour manipuler et comprendre le fichier GFF3.

Installer

Tout d'abord, configurons un environnement virtuel avec la pile SciPy installée. Ce processus peut prendre du temps s'il est construit manuellement à partir de la source, car la pile implique de nombreux packages, dont certains dépendent d'un code FORTRAN ou C externe. Ici, je recommande d'utiliser Miniconda, ce qui rend la configuration très facile.

 wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

L'indicateur -b sur la ligne bash lui indique de s'exécuter en mode batch. Une fois les commandes ci-dessus utilisées pour installer avec succès Miniconda, démarrez un nouvel environnement virtuel pour la génomique, puis installez la pile SciPy.

 mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

Notez que nous n'avons spécifié que les 3 packages que nous allons utiliser dans cet article. Si vous souhaitez que tous les packages soient répertoriés dans la pile SciPy, ajoutez-les simplement à la fin de la commande conda create .

Si vous n'êtes pas sûr du nom exact d'un paquet, essayez conda search . Activons l'environnement virtuel et démarrons IPython.

 source activate venv/ ipython

IPython est un remplacement beaucoup plus puissant de l'interface d'interpréteur Python par défaut, donc tout ce que vous faisiez dans l'interpréteur Python par défaut peut également être fait dans IPython. Je recommande fortement à tous les programmeurs Python, qui n'utilisent pas encore IPython, de l'essayer.

Télécharger le fichier d'annotations

Notre configuration étant maintenant terminée, téléchargeons le fichier d'annotation du génome humain au format GFF3.

Il s'agit d'environ 37 Mo, un très petit fichier comparé au contenu informationnel d'un génome humain, qui est d'environ 3 Go en texte brut. En effet, le fichier GFF3 ne contient que l'annotation des séquences, tandis que les données de séquence sont généralement stockées dans un autre format de fichier appelé FASTA. Si vous êtes intéressé, vous pouvez télécharger FASTA ici, mais nous n'utiliserons pas les données de séquence dans ce tutoriel.

 !wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

Le préfixe ! indique à IPython qu'il s'agit d'une commande shell au lieu d'une commande Python. Cependant, IPython peut également traiter certaines commandes shell fréquemment utilisées comme ls , pwd , rm , mkdir , rmdir même sans préfixe ! .

En jetant un coup d'œil à l'en-tête du fichier GFF, vous verrez de nombreuses lignes de métadonnées/pragmas/directives commençant par ## ou #! .

Selon le README, ## signifie que les métadonnées sont stables, tandis que #! signifie que c'est expérimental.

Plus tard, vous verrez également ### , qui est une autre directive avec une signification encore plus subtile basée sur la spécification.

Les commentaires lisibles par l'homme sont censés être après un seul # . Pour plus de simplicité, nous traiterons toutes les lignes commençant par # comme des commentaires, et les ignorerons simplement lors de notre analyse.

 ##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

La première ligne indique que la version du format GFF utilisée dans ce fichier est 3.

Viennent ensuite des résumés de toutes les régions de séquence. Comme nous le verrons plus tard, ces informations peuvent également être trouvées dans la partie corps du fichier.

Les lignes commençant par #! affiche des informations sur la version particulière du génome, GRCh38.p7, à laquelle s'applique ce fichier d'annotation.

Genome Reference Consortium (GCR) est un consortium international qui supervise les mises à jour et les améliorations de plusieurs assemblages de génomes de référence, notamment ceux de l'homme, de la souris, du poisson zèbre et du poulet.

En parcourant ce fichier, voici les premières lignes d'annotation.

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

Les colonnes sont seqid , source , type , start , end , score , strand , phase , attributes . Certains d'entre eux sont très faciles à comprendre. Prenons la première ligne comme exemple :

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

Il s'agit de l'annotation du premier chromosome avec un seqid de 1, qui part de la première base jusqu'à la 24 895 622e base.

En d'autres termes, le premier chromosome est long d'environ 25 millions de bases.

Notre analyse n'aura pas besoin des informations des trois colonnes avec une valeur de . (c'est-à-dire score, volet et phase), nous pouvons donc simplement les ignorer pour l'instant.

La dernière colonne d'attributs indique que le chromosome 1 a également trois noms d'alias, à savoir CM000663.2, chr1 et NC_000001.11. C'est essentiellement à quoi ressemble un fichier GFF3, mais nous ne les inspecterons pas ligne par ligne, il est donc temps de charger le fichier entier dans Pandas.

Pandas est bien adapté pour traiter le format GFF3 car il s'agit d'un fichier délimité par des tabulations, et Pandas prend très bien en charge la lecture de fichiers de type CSV.

Notez une exception au format délimité par des tabulations lorsque le GFF3 contient ##FASTA .

Selon la spécification, ##FASTA indique la fin d'une partie d'annotation, qui sera suivie d'une ou plusieurs séquences au format FASTA (un format non délimité par des tabulations). Mais ce n'est pas le cas du fichier GFF3 que nous allons analyser.

 In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep='\t', comment='#', low_memory=False, header=None, names=col_names)

La dernière ligne ci-dessus charge l'intégralité du fichier pandas.read_csv .

Comme il ne s'agit pas d'un fichier CSV standard, nous devons personnaliser un peu l'appel.

Tout d'abord, nous informons Pandas de l'indisponibilité des informations d'en-tête dans le GFF3 avec header=None , puis nous spécifions le nom exact de chaque colonne avec names=col_names .

Si l'argument names n'est pas spécifié, Pandas utilisera des nombres incrémentiels comme noms pour chaque colonne.

sep='\t' indique à Pandas que les colonnes sont séparées par des tabulations au lieu de séparées par des virgules. La valeur de sep= peut en fait être une expression régulière (regex). Cela devient pratique si le fichier à portée de main utilise des séparateurs différents pour chaque colonne (oh oui, ça arrive). comment='#' signifie que les lignes commençant par # sont considérées comme des commentaires et seront ignorées.

compression='gzip' indique à Pandas que le fichier d'entrée est compressé avec gzip.

De plus, pandas.read_csv dispose d'un riche ensemble de paramètres qui permettent de lire différents types de formats de fichiers de type CSV.

Le type de la valeur renvoyée est un DataFrame , qui est la structure de données la plus importante dans Pandas, utilisée pour représenter les données 2D.

Pandas a également une structure de données en Series et en Panel pour les données 1D et 3D, respectivement. Veuillez vous référer à la documentation pour une introduction aux structures de données de Pandas.

Jetons un coup d'œil aux premières entrées avec la méthode .head .

 In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

La sortie est bien formatée dans un format tabulaire avec des chaînes plus longues dans la colonne des attributs partiellement remplacées par ... .

Vous pouvez configurer Pandas pour ne pas omettre les chaînes longues avec pd.set_option('display.max_colwidth', -1) . De plus, Pandas dispose de nombreuses options personnalisables.

Ensuite, obtenons quelques informations de base sur cette trame de données avec la méthode .info .

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

Cela montre que le GFF3 a 2 601 848 lignes annotées, et chaque ligne a neuf colonnes.

Pour chaque colonne, il affiche également leurs types de données.

Ce début et cette fin sont de type int64, des entiers représentant des positions dans le génome.

Les autres colonnes sont toutes de type object , ce qui signifie probablement que leurs valeurs consistent en un mélange d'entiers, de flottants et de chaînes.

La taille de toutes les informations est d'environ 178,7+ Mo stockés en mémoire. Celui-ci s'avère plus compact que le fichier non compressé, qui fera environ 402 Mo. Une vérification rapide est illustrée ci-dessous.

 gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

D'un point de vue de haut niveau, nous avons chargé l'intégralité du fichier GFF3 dans un objet DataFrame en Python, et toutes nos analyses suivantes seront basées sur cet objet unique.

Voyons maintenant en quoi consiste le seqid de la première colonne.

 In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid est un moyen d'accéder aux données de colonne à partir d'une trame de données. Une autre façon est df['seqid'] , qui est une syntaxe plus générale, car si le nom de la colonne est un mot-clé réservé Python (par exemple class ) ou contient un . ou un espace, la première méthode ( df.seqid ) ne fonctionnera pas.

La sortie montre qu'il existe 194 seqids uniques, qui incluent les chromosomes 1 à 22, X, Y et l'ADN mitochondrial (MT) ainsi que 169 autres seqids.

Les seqids commençant par KI et GL sont des séquences d'ADN - ou des échafaudages - dans le génome qui n'ont pas été assemblés avec succès dans les chromosomes.

Pour ceux qui ne sont pas familiers avec la génomique, c'est important.

Bien que le premier projet de génome humain soit sorti il ​​y a plus de 15 ans, le génome humain actuel est encore incomplet. La difficulté d'assemblage de ces séquences est en grande partie due aux régions répétitives complexes du génome.

Ensuite, regardons la colonne source.

Le README indique que la source est un qualificateur de texte libre destiné à décrire l'algorithme ou la procédure d'exploitation qui a généré cette fonctionnalité.

 In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

Ceci est un exemple d'utilisation de la méthode value_counts , qui est extrêmement utile pour un comptage rapide des variables catégorielles.

D'après le résultat, nous pouvons voir qu'il y a sept valeurs possibles pour cette colonne, et la majorité des entrées dans le fichier GFF3 proviennent de havana, ensembl et ensembl_havana.

Vous pouvez en savoir plus sur la signification de ces sources et les relations entre elles dans cet article.

Pour garder les choses simples, nous nous concentrerons sur les entrées des sources GRCh38, havana, ensembl et ensembl_havan.a.

Quelle part du génome est incomplète ?

Les informations sur chaque chromosome entier se trouvent dans les entrées de la source GRCh38, donc commençons par filtrer le reste et attribuons le résultat filtré à une nouvelle variable gdf .

 In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

Le filtrage est facile dans Pandas.

Si vous inspectez la valeur évaluée à partir de l'expression df.source == 'GRCh38' , il s'agit d'une série de valeurs True et False pour chaque entrée avec le même index que df . Le passer à df[] ne renverra que les entrées où leurs valeurs correspondantes sont True.

Il y a 194 clés dans df[] pour lesquelles df.source == 'GRCh38' .

Comme nous l'avons vu précédemment, il existe également 194 valeurs uniques dans la colonne seqid , ce qui signifie que chaque entrée dans gdf correspond à un seqid particulier.

Ensuite, nous sélectionnons au hasard 10 entrées avec la méthode de l' sample pour regarder de plus près.

Vous pouvez voir que les séquences non assemblées sont de type supercontig tandis que les autres sont de type chromosome. Pour calculer la fraction de génome incomplète, nous devons d'abord connaître la longueur du génome entier, qui est la somme des longueurs de tous les seqids.

 In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

Dans l'extrait ci-dessus, nous avons d'abord fait une copie de gdf avec .copy() . Sinon, le gdf origine n'est qu'une tranche de df , et le modifier directement entraînerait SettingWithCopyWarning (voir ici pour plus de détails).

Nous calculons ensuite la longueur de chaque entrée et l'ajoutons à gdf en tant que nouvelle colonne nommée "longueur". La longueur totale s'avère être d'environ 3,1 milliards et la fraction de séquences non assemblées est d'environ 0,37 %.

Voici comment fonctionne le découpage dans les deux dernières commandes.

Tout d'abord, nous créons une liste de chaînes qui couvre tous les seqids de séquences bien assemblées, qui sont tous des chromosomes et des mitochondries. Nous utilisons ensuite la méthode isin pour filtrer toutes les entrées dont les seqid sont dans la liste chrs .

Un signe moins ( - ) est ajouté au début de l'index pour inverser la sélection, car nous voulons en fait tout ce qui n'est pas dans la liste (c'est-à-dire que nous voulons les non assemblés commençant par KI et GL)…

Remarque : Puisque les séquences assemblées et non assemblées se distinguent par la colonne de type, la dernière ligne peut alternativement être réécrite comme suit pour obtenir les mêmes résultats.

 gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

Combien y a-t-il de gènes ?

Ici, nous nous concentrons sur les entrées de source ensembl, havana et ensembl_havana puisqu'elles sont celles auxquelles appartiennent la majorité des entrées d'annotation.

 In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

La méthode isin est à nouveau utilisée pour le filtrage.

Ensuite, un comptage rapide des valeurs montre que la majorité des entrées sont des exons, des séquences codantes (CDS) et des régions non traduites (UTR).

Ce sont des éléments sous-géniques, mais nous recherchons principalement le nombre de gènes. Comme indiqué, il y en a 42 470, mais nous voulons en savoir plus.

Concrètement, comment s'appellent-ils et que font-ils ? Pour répondre à ces questions, nous devons examiner attentivement les informations contenues dans la colonne des attributs.

 In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

Ils sont formatés comme une liste de paires balise-valeur séparées par des points-virgules. Les informations qui nous intéressent le plus sont le nom du gène, l'ID et la description du gène, et nous les extrairons avec une expression régulière (regex).

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

Tout d'abord, nous extrayons les noms des gènes.

Dans la regex Name=(?P<gene_name>.+?); , +? est utilisé à la place de + car nous voulons qu'il soit non gourmand et que la recherche s'arrête au premier point-virgule ; sinon, le résultat correspondra jusqu'au dernier point-virgule.

De plus, la regex est d'abord compilée avec re.compile au lieu d'être utilisée directement comme dans re.search pour de meilleures performances car nous l'appliquerons à des milliers de chaînes d'attributs.

extract_gene_name sert de fonction d'assistance à utiliser dans pd.apply , qui est la méthode à utiliser lorsqu'une fonction doit être appliquée à chaque entrée d'une trame de données ou d'une série.

Dans ce cas particulier, nous voulons extraire le nom du gène pour chaque entrée dans ndf.attributes et ajouter les noms à ndf dans une nouvelle colonne appelée gene_name .

Les identifiants et la description des gènes sont extraits de la même manière.

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P<desc>.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

La regex pour RE_GENE_ID est un peu plus spécifique puisque nous savons que chaque gene_id doit commencer par ENSG , où ENS signifie ensemble et G signifie gène.

Pour les entrées qui n'ont pas de description, nous renverrons une chaîne vide. Une fois que tout est extrait, nous n'utiliserons plus la colonne des attributs, alors supprimons-la pour que les choses restent propres avec la méthode .drop :

 In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

Dans l'appel ci-dessus, les attributes indiquent la colonne spécifique que nous voulons supprimer.

axis=1 signifie que nous supprimons une colonne au lieu d'une ligne ( axis=0 par défaut).

inplace=True signifie que la suppression est effectuée sur le DataFrame lui-même au lieu de renvoyer une nouvelle copie avec la colonne spécifiée supprimée.

Un coup d'œil rapide .head montre que la colonne des attributs a bien disparu et que trois nouvelles colonnes : gene_name , gene_id et desc ont été ajoutées.

Par curiosité, voyons si tous les gene_id et gene_name sont uniques :

 In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

Étonnamment, le nombre de noms de gènes est inférieur à celui des identifiants de gènes, ce qui indique que certains noms de gènes doivent correspondre à plusieurs identifiants de gènes. Découvrons ce qu'ils sont.

 In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

Nous regroupons toutes les entrées par la valeur de gene_name , puis comptons le nombre d'éléments dans chaque groupe avec .count() .

Si vous inspectez la sortie de ndf.groupby('gene_name').count() , toutes les colonnes sont comptées pour chaque groupe, mais la plupart d'entre elles ont les mêmes valeurs.

Notez que les valeurs NA ne seront pas prises en compte lors du comptage, donc ne prenez que le nombre de la première colonne, seqid (nous utilisons .ix[:, 0] pour nous assurer qu'il n'y a pas de valeurs NA).

Triez ensuite les valeurs de comptage avec .sort_values ​​et inversez l'ordre avec .ix[::-1] .

En conséquence, un nom de gène peut être partagé avec jusqu'à sept ID de gène.

 In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

Un examen plus approfondi de tous les gènes SCARNA20 montre qu'ils sont en effet tous différents.

Bien qu'ils partagent le même nom, ils sont situés à différentes positions du génome.

Leurs descriptions, cependant, ne semblent pas très utiles pour les distinguer.

Le point ici est de savoir que les noms de gènes ne sont pas uniques pour tous les identifiants de gènes, et environ 0,15 % des noms sont partagés par plusieurs gènes.

Combien de temps dure un gène typique ?

Semblable à ce que nous avons fait lorsque nous essayions de comprendre l'incomplétude du génome, nous pouvons facilement ajouter une colonne de length à ndf :

 In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() calcule quelques statistiques simples basées sur les valeurs de longueur :

• La longueur moyenne d'un gène est d'environ 36 000 bases

• La longueur médiane d'un gène est d'environ 5 200 bases

• Les longueurs minimale et maximale des gènes sont respectivement d'environ huit et 2,3 millions de bases.

Étant donné que la moyenne est beaucoup plus grande que la médiane, cela implique que la distribution des longueurs est biaisée vers la droite. Pour avoir un aspect plus concret, traçons la distribution.

 import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

Pandas fournit une interface simple à matplotlib pour rendre le traçage très pratique avec les DataFrames ou les séries.

Dans ce cas, il est indiqué que nous voulons un tracé d'histogramme ( kind='hist' ) avec 50 bacs et que l'axe y soit sur une échelle logarithmique ( logy=True ).

À partir de l'histogramme, nous pouvons voir que la majorité des gènes se trouvent dans le premier bac. Cependant, certaines longueurs de gènes peuvent être supérieures à deux millions de bases. Découvrons ce qu'ils sont:

 In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

Comme vous pouvez le voir, le gène le plus long s'appelle CNTNAP2, qui est l'abréviation de contactin Associated Protein-like 2. Selon sa page wikipedia,

Ce gène englobe près de 1,6 % du chromosome 7 et est l'un des plus grands gènes du génome humain.

En effet! Nous venons de le vérifier nous-mêmes. En revanche, qu'en est-il des plus petits gènes ? Il s'avère qu'ils peuvent être aussi courts que huit bases.

 In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

Les durées des deux cas extrêmes sont distantes de cinq ordres de grandeur (2,3 millions contre 8), ce qui est énorme et qui peut être une indication du niveau de diversité de la vie.

Un seul gène peut être traduit en plusieurs protéines différentes via un processus appelé épissage alternatif, quelque chose que nous n'avons pas exploré. Ces informations se trouvent également dans le fichier GFF3, mais en dehors de la portée de cet article.

Distribution des gènes parmi les chromosomes

La dernière chose dont j'aimerais discuter est la distribution des gènes entre les chromosomes, qui sert également d'exemple pour l'introduction de la méthode .merge pour combiner deux DataFrames. Intuitivement, les chromosomes plus longs hébergent probablement plus de gènes. Voyons si c'est vrai.

 In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

Nous avons emprunté la variable chrs de la section précédente et l'avons utilisée pour filtrer les séquences non assemblées. Sur la base de la sortie, le plus grand chromosome 1 a en effet le plus de gènes. Bien que le chromosome Y ait le plus petit nombre de gènes, ce n'est pas le plus petit chromosome.

Notez qu'il semble n'y avoir aucun gène dans la mitochondrie (MT), ce qui n'est pas vrai.

Un peu plus de filtrage sur le premier DataFrame df renvoyé par pd.read_csv montre que tous les gènes MT proviennent de la source insdc (qui ont été filtrés auparavant lors de la génération d' edf où nous n'avons considéré que les sources de havana, ensembl ou ensembl_havana).

 In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of seqid.

 In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and gene_count .

 In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s)) 

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. In summary:

1. Environ 0,37 % du génome humain est encore incomplet même si la première ébauche est sortie il y a plus de 15 ans.
2. Il y a environ 42 000 gènes dans le génome humain d'après ce fichier GFF3 particulier que nous avons utilisé.
3. La longueur d'un gène peut aller de quelques dizaines à plus de deux millions de bases.
4. Les gènes ne sont pas répartis uniformément entre les chromosomes. Dans l'ensemble, plus le chromosome est grand, plus il héberge de gènes, mais pour un sous-ensemble de chromosomes, la corrélation peut être négative.

Le fichier GFF3 est très riche en informations d'annotation, et nous n'avons fait qu'effleurer la surface. Si vous êtes intéressé par une exploration plus approfondie, voici quelques questions avec lesquelles vous pouvez jouer :

1. Combien de transcriptions un gène possède-t-il généralement ? Quel pourcentage de gènes ont plus d'un transcrit ?
2. Combien d'isoformes une transcription a-t-elle généralement ?
3. Combien d'exons, de CDS et d'UTR une transcription a-t-elle généralement ? De quelles tailles sont-ils ?
4. Est-il possible de catégoriser les gènes en fonction de leur fonction comme décrit dans la colonne de description ?