La bioinformatique est un domaine interdisciplinaire qui développe des méthodes et des outils logiciels d'analyse et de compréhension des données biologiques.

Plus simplement, vous pouvez simplement le considérer comme 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 l'avancement rapide des technologies de séquençage d'ADN (NGS) de nouvelle génération, le volume de données génomiques a augmenté de manière exponentielle. Selon Stephens, Zachary D et al. , l'acquisition de données génomiques se fait à l'échelle de l'exaoctet par an.

SciPy rassemble de nombreux modules Python pour le calcul scientifique, ce qui est idéal pour de nombreux besoins bioinformatiques. Tweet

Dans cet article, je démontre un exemple d'analyse d'un fichier GFF3 pour le génome humain avec le SciPy Stack. Format de fonction générique 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 partie du génome est incomplète?
2. Combien de gènes y a-t-il dans le génome?
3. Combien de temps dure un gène typique?
4. À quoi ressemble la distribution des gènes entre les chromosomes?

Le dernier fichier GFF3 pour le génome humain peut être téléchargé à partir de Ici . La LISEZ-MOI Le fichier qui vient dans ce répertoire fournit une brève description de ce format de données, et une spécification plus approfondie est trouvée 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 du code FORTRAN ou C externe. Ici, je recommande d'utiliser Miniconda , ce qui rend la configuration très simple.

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

Le -b L'indicateur sur la ligne bash lui indique de s'exécuter en mode batch. Une fois que les commandes ci-dessus ont été 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 voulez que tous les packages soient répertoriés dans la pile SciPy, ajoutez-les simplement à la fin de conda create commander.

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

source activate venv/ ipython

IPython est un remplacement nettement plus puissant de la valeur par défaut Interface d'interprétation Python , donc tout ce que vous faisiez dans l'interpréteur python par défaut peut également être fait dans IPython. Je recommande vivement à chaque programmeur Python, qui n’a pas encore utilisé IPython, de l’essayer.

Téléchargez le fichier d'annotation

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

Il fait environ 37 Mo, un très petit fichier comparé au contenu informationnel d'un génome humain, soit environ 3 Go en texte brut. En effet, le fichier GFF3 ne contient que l'annotation des séquences, alors 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 didacticiel.

!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 regardant en tête du fichier GFF, vous verrez de nombreuses lignes de métadonnées / pragmas / directives commençant par ## ou #!.

Selon le LISEZ-MOI , ## 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 le spécification .

Les commentaires lisibles par l'homme sont censés être après un seul #. Pour simplifier, nous traiterons toutes les lignes commençant par # comme commentaires, et ignorez-les 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 loin, de telles informations peuvent également être trouvées dans la partie corps du fichier.

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

Consortium de référence génomique (GCR) est un consortium international qui supervise les mises à jour et les améliorations de plusieurs assemblages de génomes de référence, y compris ceux pour l'homme, la souris, le poisson zèbre et le 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 , la source , type , début , fin , But , Plage , phase , les attributs . 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 commence de la première base à la 24 895 622e base.

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

Notre analyse n’aura pas besoin d’informations provenant des trois colonnes dont la valeur est . (c'est-à-dire le score, le volet et la phase), nous pouvons donc simplement les ignorer pour le moment.

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 l’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 a un très bon support pour la lecture de fichiers de type CSV.

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

Selon le spécification , ##FASTA indique la fin d'une partie d'annotation, qui sera suivie d'une ou plusieurs séquences au format FASTA (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=' ', comment='#', low_memory=False, header=None, names=col_names)

La dernière ligne ci-dessus charge le fichier GFF3 entier avec pandas.read_csv méthode.

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

Tout d'abord, nous informons les 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 le names L'argument n'est pas spécifié, Pandas utilisera des nombres incrémentiels comme noms pour chaque colonne.

sep=' ' indique aux Pandas que les colonnes sont séparées par des tabulations au lieu de 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, cela arrive). comment='#' signifie des lignes commençant par # sont considérés comme des commentaires et seront ignorés.

compression='gzip' indique aux Pandas que le fichier d'entrée est compressé au format gzip.

En plus, pandas.read_csv dispose d'un riche ensemble de paramètres permettant 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 des données 2D.

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

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

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 joliment 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 longues chaînes avec pd.set_option('display.max_colwidth', -1). De plus, Pandas a de nombreux options qui peut être personnalisé.

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

In [20]: df.info() 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 montre é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 se composent d'un mélange d'entiers, de flottants et de chaînes.

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

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

À partir d'une 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'un dataframe. 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 espace, la première méthode (df.seqid) ne fonctionnera pas.

La sortie montre qu'il y a 194 seqids uniques, qui incluent les chromosomes 1 à 22, X, Y et l'ADN de mitochondrie (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ées avec succès dans les chromosomes.

Pour ceux qui ne connaissent pas 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 à régions répétitives complexes dans le génome .

Jetons ensuite un œil à la colonne source.

La LISEZ-MOI dit que la source est un qualificatif 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 value_counts méthode, 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 ce post .

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

Quelle partie du génome est incomplète?

Les informations sur chaque chromosome entier se trouvent dans les entrées de la source GRCh38, alors filtrons d'abord le reste et affectons 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 True et False valeurs pour chaque entrée avec le même index que df. Le passer à df[] ne renverra que les entrées dont les valeurs correspondantes sont True.

Il y a 194 touches dans df[] pour lesquels df.source == 'GRCh38'.

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

Ensuite, nous sélectionnons au hasard 10 entrées avec le sample méthode 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 chromosome. Pour calculer la fraction du génome incomplète, nous devons d'abord connaître la longueur de l'ensemble du génome, 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, l'original gdf est juste une tranche de df, et la modifier directement entraînerait SettingWithCopyWarning (voir Ici pour plus de détails).

Nous calculons ensuite la longueur de chaque entrée et la rajoutons à gdf comme une 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 le découpage fonctionne dans les deux dernières commandes.

Tout d'abord, nous créons une liste de chaînes qui couvre toutes les séquences de séquences bien assemblées, qui sont tous des chromosomes et des mitochondries. On utilise alors le isin pour filtrer toutes les entrées dont le seqid est dans le chrs liste.

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

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

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

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

Ici, nous nous concentrons sur les entrées de source ensembl, havana et ensembl_havana car elles sont à la place de 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

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

Ensuite, un comptage rapide des valeurs montre que la majorité des entrées sont l'exon, la séquence codante (CDS) et la région non traduite (UTR).

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

Plus précisément, quels sont leurs noms et que font-ils? Pour répondre à ces questions, nous devons examiner attentivement les informations de 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 sous forme de liste de paires valeur / étiquette séparées par des points-virgules. Les informations qui nous intéressent le plus sont le nom du gène, son ID et sa description, et nous les extrairons avec une expression régulière (regex).

import re RE_GENE_NAME = re.compile(r'Name=(?P.+?);') 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.+?); , +? est utilisé à la place de + parce que nous voulons qu'il ne soit pas gourmand et que la recherche s'arrête au premier point-virgule; sinon, le résultat correspondra au dernier point-virgule.

De plus, l'expression régulière est d'abord compilée avec re.compile au lieu d'être utilisé 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 rajouter 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=(?PENSG.+?);') 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.+?);') 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)

L'expression régulière pour RE_GENE_ID est un peu plus spécifique puisque nous savons que chaque gene_id doit commencer par ENSG, où ENS veux dire ensembl et G signifie gène.

Pour les entrées sans 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, attributes indique 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 rapide .head look montre que la colonne des attributs a bel et bien disparu, et trois nouvelles colonnes: gene_name, gene_id et desc ont été ajoutés.

Par curiosité, voyons si tout 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,)

De manière surprenante, le nombre de noms de gènes est plus petit que celui des ID de gène, ce qui indique que certains nom de gène doivent correspondre à plusieurs ID de gène. 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, alors 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 attentif de tous les gènes SCARNA20 montre qu’ils sont tous différents.

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

Cependant, leurs descriptions 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 ID de gènes, et environ 0,15% des noms qui 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 un length colonne vers 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 5200 bases

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

Parce 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 des DataFrames ou des séries.

Dans ce cas, cela indique que nous voulons un histogramme (kind='hist') avec 50 bins, et que l'axe y soit sur une échelle log (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 dépasser 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 est nommé 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 gros gènes du génome humain.

En effet! Nous venons de vérifier cela 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 longueurs des deux cas extrêmes sont distantes de cinq ordres de grandeur (2,3 millions contre 8), ce qui est énorme et peut être une indication du niveau de diversité de la vie.

Un seul gène peut être traduit en de nombreuses 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 je voudrais parler est la distribution des gènes entre les chromosomes, qui sert également d’exemple pour l’introduction de .merge méthode pour combiner deux DataFrames. Intuitivement, les chromosomes plus longs hébergent probablement plus de gènes. Voyons si cela 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é le chrs variable de la section précédente, et l'a utilisé pour filtrer les séquences non assemblées. Sur la base de la sortie, le plus gros 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 source insdc (qui étaient filtrés auparavant lors de la génération de edf où nous ne considérions 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...

Cet exemple montre également comment combiner deux conditions lors du filtrage avec &; l'opérateur logique pour «ou» serait |.

Notez que les parenthèses autour de chaque condition sont obligatoires, et cette partie de la syntaxe dans Pandas est différente de Python, qui aurait été littéralement and et or.

Ensuite, empruntons le gdf DataFrame de la section précédente comme source pour la longueur de chaque 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

Les colonnes qui ne sont pas pertinentes pour l'analyse sont supprimées pour plus de clarté.

Oui, .drop peut également prendre une liste de colonnes et les supprimer en une seule opération.

Notez que la ligne avec un seqid de MT est toujours là; nous y reviendrons plus tard. La prochaine opération que nous effectuerons est de fusionner les deux ensembles de données en fonction des valeurs de 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

Depuis chr_gene_counts est toujours un objet Series, qui ne prend pas en charge une opération de fusion, nous devons d'abord le convertir en objet DataFrame avec .to_frame.

.reset_index() convertit l'index d'origine (c'est-à-dire seqid) en une nouvelle colonne et réinitialise l'index actuel en tant que nombres incrémentiels basés sur 0.

La sortie de cdf.head(2) montre à quoi il ressemble. Ensuite, nous avons utilisé le .merge méthode pour combiner les deux DataFrame sur la colonne seqid (on='seqid').

Après la fusion gdf et cdf, le MT l'entrée est toujours manquante. En effet, par défaut, .merge opère une jointure interne, tandis que la jointure gauche, droite ou externe sont disponibles en réglant le how paramètre.

Veuillez vous référer au Documentation pour plus de détails.

Plus tard, vous constaterez peut-être qu'il existe également un .join méthode. .merge et .join sont similaires mais ont des API différentes.

Selon le fonctionnaire Documentation dit

La méthode DataFrame.join associée utilise la fusion en interne pour les jointures index-on-index et index-on-column (s), mais se joint aux index par défaut plutôt que d'essayer de se joindre sur des colonnes communes (le comportement par défaut pour la fusion). Si vous rejoignez sur index, vous souhaiterez peut-être utiliser DataFrame.join pour vous épargner de la saisie.

Fondamentalement, .merge est plus polyvalent et utilisé par .join.

Enfin, nous sommes prêts à calculer la corrélation entre le chromosome length et 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

Par défaut .corr calcule le Corrélation de Pearson entre toutes les paires de colonnes dans un dataframe.

Mais nous n'avons qu'une seule paire de colonnes dans ce cas, et la corrélation s'avère positive - 0,73.

En d'autres termes, plus le chromosome est gros, plus il est susceptible d'avoir plus de gènes. Tracons également les deux colonnes après avoir trié les paires de valeurs par length:

ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count', style='o-') # 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))

Comme le montre l'image ci-dessus, même s'il s'agit d'une corrélation positive dans l'ensemble, cela ne vaut pas pour tous les chromosomes. En particulier, pour les chromosomes 17, 16, 15, 14, 13, la corrélation est en fait négative, ce qui signifie que le nombre de gènes sur le chromosome diminue à mesure que la taille du chromosome augmente.

Conclusions et recherches futures

Cela termine notre tutoriel sur la manipulation d'un fichier d'annotation pour le génome humain au format GFF3 avec la pile SciPy. Les outils que nous avons principalement utilisés incluent IPython, Pandas et matplotlib. Au cours du didacticiel, nous avons non seulement appris certaines des opérations les plus courantes et les plus utiles dans Pandas, mais nous avons également répondu à des questions très intéressantes sur notre génome. En résumé:

1. Environ 0,37% du génome humain est encore incomplet, même si le premier projet est sorti il ​​y a plus de 15 ans.
2. Il y a environ 42000 gènes dans le génome humain basés sur ce particulier Fichier GFF3 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 gros, 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 venons de gratter la surface. Si vous souhaitez approfondir votre exploration, 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 a plus d'un transcrit?
2. Combien d'isoformes une transcription a-t-elle généralement?
3. Combien d'exons, CDS et UTR ont généralement une transcription? De quelles tailles s'agit-il?
4. Est-il possible de catégoriser les gènes en fonction de leur fonction comme décrit dans la colonne de description?