Un'introduzione completa al tuo genoma con lo stack SciPy

La bioinformatica è un campo interdisciplinare che sviluppa metodi e strumenti software per l'analisi e la comprensione dei dati biologici.

Detto più semplicemente, puoi semplicemente pensarlo come scienza dei dati per la biologia.

Tra i molti tipi di dati biologici, quello di genomica è uno dei più analizzati. Soprattutto con il rapido progresso delle tecnologie di sequenziamento del DNA (NGS) di prossima generazione, il volume dei dati genomici è cresciuto in modo esponenziale. Secondo Stephens, Zachary D et al., l'acquisizione dei dati genomici è su scala exabyte all'anno.

In questo post, mostro un esempio di analisi di un file GFF3 per il genoma umano con SciPy Stack. Generic Feature Format Version 3 (GFF3) è l'attuale formato di file di testo standard per la memorizzazione di caratteristiche genomiche. In particolare, in questo post imparerai come utilizzare lo stack SciPy per rispondere alle seguenti domande sul genoma umano:

1. Quanto del genoma è incompleto?
2. Quanti geni ci sono nel genoma?
3. Quanto è lungo un gene tipico?
4. Che aspetto ha la distribuzione dei geni tra i cromosomi?

L'ultimo file GFF3 per il genoma umano può essere scaricato da qui. Il file README che arriva in questa directory fornisce una breve descrizione di questo formato di dati e una specifica più completa si trova qui.

Useremo Pandas, un componente importante dello stack SciPy che fornisce strutture di dati veloci, flessibili ed espressive, per manipolare e comprendere il file GFF3.

Impostare

Per prima cosa, configuriamo un ambiente virtuale con lo stack SciPy installato. Questo processo può richiedere molto tempo se compilato manualmente dal sorgente, poiché lo stack coinvolge molti pacchetti, alcuni dei quali dipendono da FORTRAN o codice C esterno. Qui, consiglio di utilizzare Miniconda, che rende l'installazione molto semplice.

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

Il flag -b sulla riga bash gli dice di essere eseguito in modalità batch. Dopo aver utilizzato i comandi precedenti per installare correttamente Miniconda, avviare un nuovo ambiente virtuale per la genomica e quindi installare lo stack SciPy.

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

Nota che abbiamo specificato solo i 3 pacchetti che useremo in questo post. Se vuoi che tutti i pacchetti siano elencati nello stack SciPy, aggiungili semplicemente alla fine del comando conda create .

Se non sei sicuro del nome esatto di un pacchetto, prova conda search . Attiviamo l'ambiente virtuale e avviamo IPython.

 source activate venv/ ipython

IPython è un sostituto significativamente più potente dell'interfaccia dell'interprete Python predefinita, quindi qualsiasi cosa tu abbia usato per fare nell'interprete Python predefinito può essere eseguita anche in IPython. Consiglio vivamente a tutti i programmatori Python, che non hanno ancora utilizzato IPython, di provarlo.

Scarica il file di annotazione

Con la nostra configurazione ora completata, scarichiamo il file di annotazione del genoma umano in formato GFF3.

Si tratta di circa 37 MB, un file molto piccolo rispetto al contenuto informativo di un genoma umano, che è di circa 3 GB in testo normale. Questo perché il file GFF3 contiene solo l'annotazione delle sequenze, mentre i dati della sequenza sono solitamente archiviati in un altro formato di file chiamato FASTA. Se sei interessato, puoi scaricare FASTA qui, ma non useremo i dati di sequenza in questo tutorial.

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

Il prefisso ! dice a IPython che questo è un comando shell invece di un comando Python. Tuttavia, IPython può anche elaborare alcuni comandi shell usati di frequente come ls , pwd , rm , mkdir , rmdir anche senza un prefisso ! .

Dando un'occhiata alla testa del file GFF, vedrai molte righe di metadati/pragma/direttive che iniziano con ## o #! .

Secondo il README, ## significa che i metadati sono stabili, mentre #! significa che è sperimentale.

Più avanti vedrai anche ### , che è un'altra direttiva con un significato ancora più sottile basato sulla specifica.

I commenti leggibili dall'uomo dovrebbero essere dopo un singolo # . Per semplicità, tratteremo tutte le righe che iniziano con # come commenti e le ignoreremo semplicemente durante la nostra analisi.

 ##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 prima riga indica che la versione del formato GFF utilizzata in questo file è 3.

Di seguito sono riportati i riepiloghi di tutte le regioni della sequenza. Come vedremo in seguito, tali informazioni si trovano anche nella parte del corpo del file.

Le righe che iniziano con #! mostra informazioni sulla build particolare del genoma, GRCh38.p7, a cui si applica questo file di annotazione.

Il Genome Reference Consortium (GCR) è un consorzio internazionale che supervisiona gli aggiornamenti e i miglioramenti a diversi gruppi di genomi di riferimento, inclusi quelli per l'uomo, il topo, il pesce zebra e il pollo.

Esaminando questo file, ecco le prime righe di annotazione.

 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 ...

Le colonne sono seqid , sorgente , tipo , inizio , fine , punteggio , strand , fase , attributi . Alcuni di loro sono molto facili da capire. Prendi la prima riga come esempio:

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

Questa è l'annotazione del primo cromosoma con seqid 1, che parte dalla prima base alla 24.895.622a base.

In altre parole, il primo cromosoma è lungo circa 25 milioni di basi.

La nostra analisi non avrà bisogno di informazioni dalle tre colonne con un valore di . (vale a dire punteggio, filamento e fase), quindi per ora possiamo semplicemente ignorarli.

L'ultima colonna degli attributi dice che anche il cromosoma 1 ha tre nomi alias, ovvero CM000663.2, chr1 e NC_000001.11. Questo è fondamentalmente l'aspetto di un file GFF3, ma non li ispezioniamo riga per riga, quindi è ora di caricare l'intero file in Pandas.

Pandas è adatto per gestire il formato GFF3 perché è un file delimitato da tabulazioni e Pandas ha un ottimo supporto per la lettura di file simili a CSV.

Nota un'eccezione al formato delimitato da tabulazioni è quando il GFF3 contiene ##FASTA .

Secondo la specifica, ##FASTA indica la fine di una porzione di annotazione, che sarà seguita da una o più sequenze in formato FASTA (non delimitato da tabulazioni). Ma questo non è il caso del file GFF3 che andremo ad analizzare.

 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)

L'ultima riga sopra carica l'intero file pandas.read_csv .

Poiché non è un file CSV standard, è necessario personalizzare un po' la chiamata.

In primo luogo, informiamo Panda dell'indisponibilità delle informazioni di intestazione nel GFF3 con header=None , quindi specifichiamo il nome esatto per ciascuna colonna con names=col_names .

Se l'argomento dei names non è specificato, Panda utilizzerà i numeri incrementali come nomi per ciascuna colonna.

sep='\t' dice a Pandas che le colonne sono separate da tabulazioni invece che da virgole. Il valore di sep= può effettivamente essere un'espressione regolare (regex). Questo diventa utile se il file a portata di mano utilizza separatori diversi per ogni colonna (oh sì, succede). comment='#' significa che le righe che iniziano con # sono considerate commenti e verranno ignorate.

compression='gzip' dice a Pandas che il file di input è compresso con gzip.

Inoltre, pandas.read_csv ha un ricco set di parametri che consente di leggere diversi tipi di formati di file simili a CSV.

Il tipo del valore restituito è un DataFrame , che è la struttura dati più importante in Pandas, utilizzata per rappresentare dati 2D.

Pandas ha anche una struttura di dati Series e Panel rispettivamente per i dati 1D e 3D. Fare riferimento alla documentazione per un'introduzione alle strutture dati di Pandas.

Diamo un'occhiata alle prime voci con il metodo .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

L'output è ben formattato in un formato tabulare con stringhe più lunghe nella colonna degli attributi parzialmente sostituite con ... .

Puoi impostare Panda per non omettere stringhe lunghe con pd.set_option('display.max_colwidth', -1) . Inoltre, Panda ha molte opzioni che possono essere personalizzate.

Successivamente, otteniamo alcune informazioni di base su questo dataframe con il metodo .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

Ciò mostra che GFF3 ha 2.601.848 righe annotate e ogni riga ha nove colonne.

Per ogni colonna, mostra anche i loro tipi di dati.

Quell'inizio e la fine sono di tipo int64, numeri interi che rappresentano posizioni nel genoma.

Le altre colonne sono tutte di tipo object , il che probabilmente significa che i loro valori sono costituiti da una combinazione di numeri interi, float e stringhe.

La dimensione di tutte le informazioni è di circa 178,7+ MB archiviati in memoria. Questo risulta essere più compatto del file non compresso, che sarà di circa 402 MB. Di seguito viene mostrata una rapida verifica.

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

Da una vista di alto livello, abbiamo caricato l'intero file GFF3 in un oggetto DataFrame in Python e tutte le nostre analisi successive saranno basate su questo singolo oggetto.

Ora, vediamo di cosa tratta la prima colonna seqid.

 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 è un modo per accedere ai dati delle colonne da un dataframe. Un altro modo è df['seqid'] , che è una sintassi più generale, perché se il nome della colonna è una parola chiave riservata Python (ad esempio class ) o contiene un . o carattere spazio, il primo modo ( df.seqid ) non funzionerà.

L'output mostra che ci sono 194 seqidi univoci, che includono i cromosomi da 1 a 22, X, Y e il DNA dei mitocondri (MT), nonché altri 169 seqidi.

I seqidi che iniziano con KI e GL sono sequenze di DNA - o scaffold - nel genoma che non sono state assemblate con successo nei cromosomi.

Per coloro che non hanno familiarità con la genomica, questo è importante.

Sebbene la prima bozza del genoma umano sia stata pubblicata più di 15 anni fa, l'attuale genoma umano è ancora incompleto. La difficoltà nell'assemblare queste sequenze è in gran parte dovuta a complesse regioni ripetitive nel genoma.

Quindi, diamo un'occhiata alla colonna sorgente.

Il README dice che il sorgente è un qualificatore di testo libero inteso a descrivere l'algoritmo o la procedura operativa che ha generato questa caratteristica.

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

Questo è un esempio dell'uso del metodo value_counts , estremamente utile per un conteggio rapido di variabili categoriali.

Dal risultato, possiamo vedere che ci sono sette valori possibili per questa colonna e la maggior parte delle voci nel file GFF3 provengono da havana, ensembl e ensembl_havana.

Puoi saperne di più sul significato di queste fonti e sulle relazioni tra loro in questo post.

Per semplificare le cose, ci concentreremo sulle voci dalle fonti GRCh38, havana, ensembl e ensembl_havan.a.

Quanto del genoma è incompleto?

Le informazioni su ogni intero cromosoma sono nelle voci dalla fonte GRCh38, quindi filtriamo prima il resto e assegniamo il risultato filtrato a una nuova variabile 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...

Il filtraggio è facile in Panda.

Se controlli il valore valutato dall'espressione df.source == 'GRCh38' , è una serie di valori True e False per ogni voce con lo stesso indice di df . Passandolo a df[] verranno restituite solo le voci in cui i valori corrispondenti sono True.

Ci sono 194 chiavi in df[] per le quali df.source == 'GRCh38' .

Come abbiamo visto in precedenza, ci sono anche 194 valori univoci nella colonna seqid , il che significa che ogni voce in gdf corrisponde a un particolare seqid.

Quindi selezioniamo casualmente 10 voci con il metodo di sample per dare un'occhiata più da vicino.

Si può notare che le sequenze non assemblate sono di tipo supercontig mentre le altre sono di cromosoma. Per calcolare la frazione incompleta del genoma, dobbiamo prima conoscere la lunghezza dell'intero genoma, che è la somma delle lunghezze di tutti i seqidi.

 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

Nello snippet sopra, abbiamo prima creato una copia di gdf con .copy() . Altrimenti, il gdf originale è solo una fetta di df e modificarlo direttamente risulterebbe in SettingWithCopyWarning (vedi qui per maggiori dettagli).

Quindi calcoliamo la lunghezza di ogni voce e la aggiungiamo a gdf come una nuova colonna denominata "lunghezza". La lunghezza totale risulta essere di circa 3,1 miliardi e la frazione di sequenze non assemblate è di circa lo 0,37%.

Ecco come funziona lo slicing negli ultimi due comandi.

Innanzitutto, creiamo un elenco di stringhe che copre tutti i seqidi di sequenze ben assemblate, che sono tutti cromosomi e mitocondri. Usiamo quindi il metodo isin per filtrare tutte le voci il cui seqid è nell'elenco chrs .

Un segno meno ( - ) viene aggiunto all'inizio dell'indice per invertire la selezione, perché in realtà vogliamo tutto ciò che non è nell'elenco (cioè vogliamo quelli smontati che iniziano con KI e GL)...

Nota: Poiché le sequenze assemblate e non assemblate si distinguono per la colonna del tipo, l'ultima riga può in alternativa essere riscritta come segue per ottenere gli stessi risultati.

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

Quanti geni ci sono?

Qui ci concentriamo sulle voci da source ensembl, havana e ensembl_havana poiché sono dove appartiene la maggior parte delle voci di annotazione.

 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

Il metodo isin viene utilizzato di nuovo per il filtraggio.

Quindi, un rapido conteggio dei valori mostra che la maggior parte delle voci sono esone, sequenza di codifica (CDS) e regione non tradotta (UTR).

Questi sono elementi del sottogene, ma stiamo principalmente cercando il conteggio dei geni. Come mostrato, sono 42.470, ma vogliamo saperne di più.

Nello specifico, quali sono i loro nomi e cosa fanno? Per rispondere a queste domande, dobbiamo esaminare attentamente le informazioni nella colonna degli attributi.

 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)

Sono formattati come elenchi separati da punto e virgola di coppie di valori di tag. Le informazioni che ci interessano di più sono il nome del gene, l'ID del gene e la descrizione, e le estrarremo con l'espressione regolare (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)

Innanzitutto, estraiamo i nomi dei geni.

Nella regex Name=(?P<gene_name>.+?); , +? viene utilizzato al posto di + perché vogliamo che non sia avido e che la ricerca si fermi al primo punto e virgola; in caso contrario, il risultato corrisponderà fino all'ultimo punto e virgola.

Inoltre, la regex viene prima compilata con re.compile invece di essere utilizzata direttamente come in re.search per prestazioni migliori perché la applicheremo a migliaia di stringhe di attributi.

extract_gene_name funge da funzione di supporto da utilizzare in pd.apply , che è il metodo da utilizzare quando è necessario applicare una funzione su ogni voce di un dataframe o serie.

In questo caso particolare, vogliamo estrarre il nome del gene per ogni voce in ndf.attributes e aggiungere nuovamente i nomi a ndf in una nuova colonna chiamata gene_name .

Gli ID e la descrizione del gene vengono estratti in modo simile.

 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 per RE_GENE_ID è un po' più specifica poiché sappiamo che ogni gene_id deve iniziare con ENSG , dove ENS significa ensembl e G significa gene.

Per le voci che non hanno alcuna descrizione, restituiremo una stringa vuota. Dopo che tutto è stato estratto, non useremo più la colonna degli attributi, quindi lasciamo perdere per mantenere le cose belle e pulite con il metodo .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...

Nella chiamata sopra, attributes indica la colonna specifica che vogliamo eliminare.

axis=1 significa che stiamo eliminando una colonna invece di una riga ( axis=0 per impostazione predefinita).

inplace=True significa che il rilascio viene eseguito sul DataFrame stesso invece di restituire una nuova copia con la colonna specificata eliminata.

Un rapido sguardo .head mostra che la colonna degli attributi è effettivamente scomparsa e sono state aggiunte tre nuove colonne: gene_name , gene_id e desc .

Per curiosità, vediamo se tutti gene_id e gene_name sono unici:

 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,)

Sorprendentemente, il numero dei nomi dei geni è inferiore a quello degli ID dei geni, indicando che alcuni nomi_gene devono corrispondere a più ID dei geni. Scopriamo quali sono.

 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

Raggruppiamo tutte le voci in base al valore di gene_name , quindi contiamo il numero di elementi in ciascun gruppo con .count() .

Se controlli l'output da ndf.groupby('gene_name').count() , tutte le colonne vengono conteggiate per ogni gruppo, ma la maggior parte di esse ha gli stessi valori.

Nota che i valori NA non verranno considerati durante il conteggio, quindi prendi solo il conteggio della prima colonna, seqid ( usiamo .ix[:, 0] per assicurarci che non ci siano valori NA).

Quindi ordina i valori di conteggio con .sort_values ​​e inverti l'ordine con .ix[::-1] .

Nel risultato, un nome di gene può essere condiviso con un massimo di sette ID di gene.

 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]

Uno sguardo più da vicino a tutti i geni SCARNA20 mostra che sono davvero tutti diversi.

Sebbene condividano lo stesso nome, si trovano in posizioni diverse del genoma.

Le loro descrizioni, tuttavia, non sembrano molto utili per distinguerli.

Il punto qui è sapere che i nomi dei geni non sono univoci per tutti gli ID dei geni e circa lo 0,15% dei nomi che sono condivisi da più geni.

Quanto dura un gene tipico?

Simile a quello che abbiamo fatto quando stavamo cercando di capire l'incompletezza del genoma, possiamo facilmente aggiungere una colonna di length a 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() calcola alcune semplici statistiche in base ai valori di lunghezza:

• La lunghezza media di un gene è di circa 36.000 basi

• La lunghezza mediana di un gene è di circa 5.200 basi

• Le lunghezze minime e massime dei geni sono rispettivamente di circa otto e 2,3 milioni di basi.

Poiché la media è molto più grande della mediana, implica che la distribuzione della lunghezza è distorta a destra. Per avere uno sguardo più concreto, tracciamo la distribuzione.

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

Pandas fornisce una semplice interfaccia per matplotlib per rendere la stampa molto utile con DataFrames o serie.

In questo caso, dice che vogliamo un grafico dell'istogramma ( kind='hist' ) con 50 bin e lascia che l'asse y sia su una scala logaritmica ( logy=True ).

Dall'istogramma, possiamo vedere che la maggior parte dei geni si trova all'interno del primo bin. Tuttavia, alcune lunghezze dei geni possono superare i due milioni di basi. Scopriamo insieme quali sono:

 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

Come puoi vedere, il gene più lungo si chiama CNTNAP2, che è l'abbreviazione di contactin-associated protein-like 2. Secondo la sua pagina di wikipedia,

Questo gene comprende quasi l'1,6% del cromosoma 7 ed è uno dei geni più grandi nel genoma umano.

Infatti! Lo abbiamo appena verificato noi stessi. Al contrario, che dire dei geni più piccoli? Si scopre che possono essere più brevi di otto basi.

 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

Le lunghezze dei due casi estremi sono distanti cinque ordini di grandezza (2,3 milioni contro 8), che è enorme e che può essere un'indicazione del livello di diversità della vita.

Un singolo gene può essere tradotto in molte proteine ​​diverse attraverso un processo chiamato splicing alternativo, qualcosa che non abbiamo esplorato. Tali informazioni sono anche all'interno del file GFF3, ma al di fuori dello scopo di questo post.

Distribuzione genica tra cromosomi

L'ultima cosa che vorrei discutere è la distribuzione del gene tra i cromosomi, che serve anche come esempio per l'introduzione del metodo .merge per combinare due DataFrame. Intuitivamente, i cromosomi più lunghi probabilmente ospitano più geni. Vediamo se è vero.

 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

Abbiamo preso in prestito la variabile chrs dalla sezione precedente e l'abbiamo usata per filtrare le sequenze non assemblate. Sulla base dell'output, il più grande cromosoma 1 ha effettivamente il maggior numero di geni. Sebbene il cromosoma Y abbia il minor numero di geni, non è il cromosoma più piccolo.

Si noti che sembra che non ci siano geni nel mitocondrio (MT), il che non è vero.

Un po' più di filtraggio sul primo df DataFrame restituito da pd.read_csv mostra che tutti i geni MT provengono dalla sorgente insdc (che è stata filtrata prima durante la generazione di edf dove abbiamo considerato solo le sorgenti di havana, ensembl o 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. Circa lo 0,37% del genoma umano è ancora incompleto anche se la prima bozza è stata pubblicata più di 15 anni fa.
2. Ci sono circa 42.000 geni nel genoma umano basati su questo particolare file GFF3 che abbiamo usato.
3. La lunghezza di un gene può variare da poche decine a oltre due milioni di basi.
4. I geni non sono distribuiti uniformemente tra i cromosomi. Nel complesso, più grande è il cromosoma, più geni ospita, ma per un sottoinsieme di cromosomi la correlazione può essere negativa.

Il file GFF3 è molto ricco di informazioni di annotazione e abbiamo appena scalfito la superficie. Se sei interessato ad ulteriori esplorazioni, ecco alcune domande con cui puoi giocare:

1. Quante trascrizioni ha in genere un gene? Quale percentuale di geni ha più di 1 trascrizione?
2. Quante isoforme ha in genere una trascrizione?
3. Quanti esoni, CDS e UTR ha in genere una trascrizione? Che taglie sono?
4. È possibile classificare i geni in base alla loro funzione come descritto nella colonna della descrizione?