O introducere cuprinzătoare la genomul dvs. cu stiva SciPy

Bioinformatica este un domeniu interdisciplinar care dezvoltă metode și instrumente software pentru analiza și înțelegerea datelor biologice.

Mai simplu spus, te poți gândi pur și simplu la asta ca fiind știința datelor pentru biologie.

Dintre numeroasele tipuri de date biologice, datele genomice sunt una dintre cele mai analizate. În special odată cu progresul rapid al tehnologiilor de secvențiere a ADN-ului (NGS) de generație următoare, volumul datelor genomice a crescut exponențial. Potrivit lui Stephens, Zachary D și colab., achiziția datelor genomice este la scara exabyte-pe-an.

În această postare, demonstrez un exemplu de analiză a unui fișier GFF3 pentru genomul uman cu SciPy Stack. Formatul de caracteristici generice Versiunea 3 (GFF3) este formatul de fișier text standard actual pentru stocarea caracteristicilor genomice. În special, în această postare veți învăța cum să utilizați stiva SciPy pentru a răspunde la următoarele întrebări despre genomul uman:

1. Cât de mult din genom este incomplet?
2. Câte gene sunt în genom?
3. Cât durează o genă tipică?
4. Cum arată distribuția genelor între cromozomi?

Cel mai recent fișier GFF3 pentru genomul uman poate fi descărcat de aici. Fișierul README care vine în acest director oferă o scurtă descriere a acestui format de date și o specificație mai detaliată este găsită aici.

Vom folosi Pandas, o componentă majoră a stivei SciPy, care oferă structuri de date rapide, flexibile și expresive, pentru a manipula și înțelege fișierul GFF3.

Înființat

În primul rând, să instalăm un mediu virtual cu stiva SciPy instalată. Acest proces poate consuma mult timp dacă este construit din sursă manual, deoarece stiva implică multe pachete – dintre care unele depind de codul extern FORTRAN sau C. Aici, recomand să utilizați Miniconda, ceea ce face configurarea foarte ușoară.

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

Indicatorul -b de pe linia bash îi spune să se execute în modul batch. După ce comenzile de mai sus sunt folosite pentru a instala cu succes Miniconda, porniți un nou mediu virtual pentru genomică și apoi instalați stiva SciPy.

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

Rețineți că am specificat doar cele 3 pachete pe care le vom folosi în această postare. Dacă doriți toate pachetele listate în stiva SciPy, pur și simplu adăugați-le la sfârșitul comenzii conda create .

Dacă nu sunteți sigur de numele exact al unui pachet, încercați conda search . Să activăm mediul virtual și să pornim IPython.

 source activate venv/ ipython

IPython este un înlocuitor semnificativ mai puternic al interfeței implicite de interpret Python, astfel încât orice făceai în interpretul Python implicit poate fi făcut și în IPython. Recomand cu căldură fiecărui programator Python, care nu a folosit încă IPython, să încerce.

Descărcați fișierul de adnotare

Cu configurarea noastră acum finalizată, să descarcăm fișierul de adnotare a genomului uman în format GFF3.

Este de aproximativ 37 MB, un fișier foarte mic în comparație cu conținutul de informații al unui genom uman, care este de aproximativ 3 GB în text simplu. Asta pentru că fișierul GFF3 conține doar adnotarea secvențelor, în timp ce datele secvenței sunt de obicei stocate într-un alt format de fișier numit FASTA. Dacă sunteți interesat, puteți descărca FASTA de aici, dar nu vom folosi datele secvenței din acest tutorial.

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

Prefixul ! spune IPython că aceasta este o comandă shell în loc de o comandă Python. Cu toate acestea, IPython poate procesa și unele comenzi shell utilizate frecvent, cum ar fi ls , pwd , rm , mkdir , rmdir chiar și fără un prefix ! .

Aruncând o privire la capul fișierului GFF, veți vedea multe linii de metadate/pragma/directive care încep cu ## sau #! .

Conform README, ## înseamnă că metadatele sunt stabile, în timp ce #! înseamnă că este experimental.

Mai târziu, veți vedea, de asemenea, ### , care este o altă directivă cu un sens încă mai subtil bazat pe specificație.

Comentariile care pot fi citite de oameni ar trebui să fie după un singur # . Pentru simplitate, vom trata toate rândurile care încep cu # ca comentarii și le vom ignora pur și simplu în timpul analizei noastre.

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

Prima linie indică faptul că versiunea formatului GFF utilizată în acest fișier este 3.

În continuare sunt rezumate ale tuturor regiunilor secvenței. După cum vom vedea mai târziu, astfel de informații pot fi găsite și în corpul fișierului.

Rândurile care încep cu #! arată informații despre construcția particulară a genomului, GRCh38.p7, căreia i se aplică acest fișier de adnotare.

Genome Reference Consortium (GCR) este un consorțiu internațional, care supraveghează actualizările și îmbunătățirile mai multor ansambluri de genom de referință, inclusiv cele pentru om, șoarece, pește-zebra și pui.

Scanând prin acest fișier, iată primele câteva rânduri de adnotare.

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

Coloanele sunt seqid , source , type , start , end , score , strand , phase , attributes . Unele dintre ele sunt foarte ușor de înțeles. Luați prima linie ca exemplu:

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

Aceasta este adnotarea primului cromozom cu un seqid de 1, care începe de la prima bază până la baza 24.895.622.

Cu alte cuvinte, primul cromozom are o lungime de aproximativ 25 de milioane de baze.

Analiza noastră nu va avea nevoie de informații din cele trei coloane cu valoarea . (adică scor, fir și fază), așa că pur și simplu le putem ignora pentru moment.

Ultima coloană cu atribute spune că Cromozomul 1 are și trei nume de alias, și anume CM000663.2, chr1 și NC_000001.11. Practic, așa arată un fișier GFF3, dar nu le vom inspecta linie cu linie, așa că este timpul să încărcăm întregul fișier în Pandas.

Pandas este potrivit pentru a face față formatului GFF3, deoarece este un fișier delimitat de tabulatori, iar Pandas are un suport foarte bun pentru citirea fișierelor asemănătoare CSV.

Rețineți că o excepție de la formatul delimitat de tabulatori este atunci când GFF3 conține ##FASTA .

Conform specificației, ##FASTA indică sfârșitul unei porțiuni de adnotare, care va fi urmată de una sau mai multe secvențe în format FASTA (un format nedelimitat de tabulatori). Dar nu este cazul fișierului GFF3 pe care îl vom analiza.

 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)

Ultima linie de mai sus încarcă întregul fișier pandas.read_csv .

Deoarece nu este un fișier CSV standard, trebuie să personalizăm puțin apelul.

Mai întâi, informăm Pandas despre indisponibilitatea informațiilor de antet în GFF3 cu header=None și apoi specificăm numele exact pentru fiecare coloană cu names=col_names .

Dacă argumentul names nu este specificat, Pandas va folosi numere incrementale ca nume pentru fiecare coloană.

sep='\t' îi spune lui Pandas că coloanele sunt separate prin tabulatori în loc să fie separate prin virgulă. Valoarea pentru sep= poate fi de fapt o expresie regulată (regex). Acest lucru devine util dacă fișierul la îndemână folosește separatori diferiți pentru fiecare coloană (oh, da, asta se întâmplă). comment='#' înseamnă că liniile care încep cu # sunt considerate comentarii și vor fi ignorate.

compression='gzip' îi spune lui Pandas că fișierul de intrare este comprimat cu gzip.

În plus, pandas.read_csv are un set bogat de parametri care permite citirea diferitelor tipuri de formate de fișiere asemănătoare CSV.

Tipul valorii returnate este un DataFrame , care este cea mai importantă structură de date din Pandas, utilizată pentru reprezentarea datelor 2D.

Pandas are, de asemenea, o structură de date Series și Panel date 1D și, respectiv, 3D. Vă rugăm să consultați documentația pentru o introducere în structurile de date ale Pandas.

Să aruncăm o privire la primele câteva intrări cu metoda .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

Ieșirea este frumos formatată într-un format tabelar cu șiruri mai lungi în coloana de atribute înlocuite parțial cu ... .

Puteți seta Pandas să nu omite șiruri lungi cu pd.set_option('display.max_colwidth', -1) . În plus, Pandas are multe opțiuni care pot fi personalizate.

În continuare, să obținem câteva informații de bază despre acest cadru de date cu metoda .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

Aceasta arată că GFF3 are 2.601.848 de linii adnotate și fiecare linie are nouă coloane.

Pentru fiecare coloană, arată și tipurile de date ale acestora.

Acel început și sfârșit sunt de tip int64, numere întregi reprezentând poziții în genom.

Celelalte coloane sunt toate de tip object , ceea ce înseamnă probabil că valorile lor constau dintr-un amestec de numere întregi, flotanți și șiruri de caractere.

Dimensiunea tuturor informațiilor este de aproximativ 178,7+ MB stocate în memorie. Acesta se dovedește a fi mai compact decât fișierul necomprimat, care va avea aproximativ 402 MB. O verificare rapidă este prezentată mai jos.

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

Dintr-o vedere la nivel înalt, am încărcat întregul fișier GFF3 într-un obiect DataFrame în Python și toate analizele noastre următoare se vor baza pe acest singur obiect.

Acum, să vedem despre ce este vorba despre prima coloană 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 este o modalitate de a accesa datele coloanei dintr-un cadru de date. O altă modalitate este df['seqid'] , care este o sintaxă mai generală, deoarece dacă numele coloanei este un cuvânt cheie rezervat Python (de exemplu class ) sau conține un . sau caracter de spațiu, prima modalitate ( df.seqid ) nu va funcționa.

Rezultatele arată că există 194 de seqide unice, care includ cromozomii 1 până la 22, X, Y și ADN mitocondrial (MT), precum și alte 169 de seqide.

Seqidele care încep cu KI și GL sunt secvențe de ADN - sau schele - din genom care nu au fost asamblate cu succes în cromozomi.

Pentru cei care nu sunt familiarizați cu genomica, acest lucru este important.

Deși prima schiță a genomului uman a apărut acum mai bine de 15 ani, genomul uman actual este încă incomplet. Dificultatea de a asambla aceste secvențe se datorează în mare măsură regiunilor complexe repetitive din genom.

În continuare, să aruncăm o privire la coloana sursă.

README spune că sursa este un calificativ de text liber destinat să descrie algoritmul sau procedura de operare care a generat această caracteristică.

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

Acesta este un exemplu de utilizare a metodei value_counts , care este extrem de utilă pentru o numărare rapidă a variabilelor categoriale.

Din rezultat, putem vedea că există șapte valori posibile pentru această coloană, iar majoritatea intrărilor din fișierul GFF3 provin din havana, ensembl și ensembl_havana.

Puteți afla mai multe despre ce înseamnă aceste surse și relațiile dintre ele în această postare.

Pentru a menține lucrurile simple, ne vom concentra pe intrările din sursele GRCh38, havana, ensembl și ensembl_havan.a.

Cât de mult din genom este incomplet?

Informațiile despre fiecare cromozom întreg se află în intrările din sursa GRCh38, așa că haideți mai întâi să filtrăm restul și să atribuim rezultatul filtrat unei noi 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...

Filtrarea este ușoară în Pandas.

Dacă inspectați valoarea evaluată din expresia df.source == 'GRCh38' , este o serie de valori True și False pentru fiecare intrare cu același index ca df . Transmiterea lui df[] va returna doar acele intrări în care valorile corespunzătoare sunt adevărate.

Există 194 de chei în df[] pentru care df.source == 'GRCh38' .

După cum am văzut anterior, există și 194 de valori unice în coloana seqid , ceea ce înseamnă că fiecare intrare din gdf corespunde unui anumit seqid.

Apoi selectăm aleatoriu 10 intrări cu metoda sample pentru a ne uita mai atent.

Puteți vedea că secvențele neasamblate sunt de tip supercontig în timp ce celelalte sunt de cromozom. Pentru a calcula fracțiunea de genom care este incompletă, trebuie mai întâi să cunoaștem lungimea întregului genom, care este suma lungimilor tuturor seqidelor.

 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

În fragmentul de mai sus, am făcut mai întâi o copie a gdf cu .copy() . În caz contrar, gdf original este doar o porțiune de df , iar modificarea directă ar duce la SettingWithCopyWarning (vezi aici pentru mai multe detalii).

Apoi calculăm lungimea fiecărei intrări și o adăugăm înapoi la gdf ca o nouă coloană numită „lungime”. Lungimea totală se dovedește a fi de aproximativ 3,1 miliarde, iar fracțiunea secvențelor neasamblate este de aproximativ 0,37%.

Iată cum funcționează tăierea în ultimele două comenzi.

În primul rând, creăm o listă de șiruri care acoperă toate secvențele de secvențe bine asamblate, care sunt toți cromozomi și mitocondrii. Apoi folosim metoda isin pentru a filtra toate intrările al căror seqid se află în lista chrs .

Se adaugă un semn minus ( - ) la începutul indexului pentru a inversa selecția, pentru că de fapt vrem tot ceea ce nu este în listă (adică vrem ca cele neasamblate începând cu KI și GL)...

Notă: Deoarece secvențele asamblate și neasamblate se disting prin coloana tip, ultima linie poate fi rescrisă alternativ după cum urmează pentru a obține aceleași rezultate.

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

Câte gene sunt?

Aici ne concentrăm asupra intrărilor din sursă ensembl, havana și ensembl_havana, deoarece acestea sunt locul unde aparțin majoritatea intrărilor de adnotare.

 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

Metoda isin este folosită din nou pentru filtrare.

Apoi, o contorizare rapidă a valorii arată că majoritatea intrărilor sunt exon, secvență de codificare (CDS) și regiune netradusă (UTR).

Acestea sunt elemente sub-gene, dar căutăm în principal numărul de gene. După cum se arată, sunt 42.470, dar vrem să știm mai multe.

Mai exact, care sunt numele lor și ce fac ei? Pentru a răspunde la aceste întrebări, trebuie să ne uităm îndeaproape la informațiile din coloana atribute.

 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)

Sunt formatate ca o listă de perechi etichetă-valoare, separate prin punct și virgulă. Informația care ne interesează cel mai mult este numele genei, ID-ul genei și descrierea și le vom extrage cu expresie regulată (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)

În primul rând, extragem numele genelor.

În regex Name=(?P<gene_name>.+?); , +? este folosit în loc de + pentru că vrem să fie non-lacom și să lăsăm căutarea să se oprească la primul punct și virgulă; în caz contrar, rezultatul se va potrivi până la ultimul punct și virgulă.

De asemenea, regex este mai întâi compilat cu re.compile în loc să fie folosit direct ca în re.search pentru o performanță mai bună, deoarece îl vom aplica la mii de șiruri de atribute.

extract_gene_name servește ca o funcție de ajutor pentru a fi utilizată în pd.apply , care este metoda de utilizat atunci când o funcție trebuie să fie aplicată la fiecare intrare a unui cadru de date sau a unei serii.

În acest caz particular, dorim să extragem numele genei pentru fiecare intrare din ndf.attributes și să adăugăm numele înapoi la ndf într-o nouă coloană numită gene_name .

ID-urile genelor și descrierea sunt extrase într-un mod similar.

 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)

Regex pentru RE_GENE_ID este puțin mai specific, deoarece știm că fiecare gene_id trebuie să înceapă cu ENSG , unde ENS înseamnă ansamblu și G înseamnă genă.

Pentru intrările care nu au nicio descriere, vom returna un șir gol. După ce totul este extras, nu vom mai folosi coloana de atribute, așa că haideți să o renunțăm pentru a păstra lucrurile frumoase și curate cu metoda .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...

În apelul de mai sus, attributes indică coloana specifică pe care dorim să o renunțăm.

axis=1 înseamnă că aruncăm o coloană în loc de un rând ( axis=0 în mod implicit).

inplace=True înseamnă că drop-ul este operat pe DataFrame în sine, în loc să returneze o nouă copie cu coloana specificată eliminată.

O privire rapidă .head arată că coloana de atribute a dispărut într-adevăr și au fost adăugate trei coloane noi: gene_name , gene_id și desc .

De curiozitate, să vedem dacă toate gene_id și gene_name sunt unice:

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

În mod surprinzător, numărul de nume de gene este mai mic decât cel al ID-urilor genelor, ceea ce indică faptul că unele nume_gene trebuie să corespundă cu mai multe ID-uri de gene. Să aflăm care sunt.

 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

Grupăm toate intrările după valoarea gene_name , apoi numărăm numărul de articole din fiecare grup cu .count() .

Dacă inspectați rezultatul de la ndf.groupby('gene_name').count() , toate coloanele sunt numărate pentru fiecare grup, dar cele mai multe dintre ele au aceleași valori.

Rețineți că valorile NA nu vor fi luate în considerare la numărare, așa că luați numai numărul primei coloane, seqid (folosim .ix[:, 0] pentru a ne asigura că nu există valori NA).

Apoi sortați valorile numărului cu .sort_values ​​și inversați ordinea cu .ix[::-1] .

În rezultat, un nume de genă poate fi partajat cu până la șapte ID-uri de 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]

O privire mai atentă la toate genele SCARNA20 arată că acestea sunt într-adevăr toate diferite.

Deși au același nume, ele sunt localizate în poziții diferite ale genomului.

Descrierile lor, însă, nu par foarte utile pentru a le distinge.

Ideea aici este să știm că numele genelor nu sunt unice pentru toate ID-urile genelor și aproximativ 0,15% dintre numele care sunt împărtășite de mai multe gene.

Cât durează o genă tipică?

Similar cu ceea ce am făcut când am încercat să înțelegem caracterul incomplet al genomului, putem adăuga cu ușurință o coloană de length la 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() calculează câteva statistici simple pe baza valorilor lungimii:

• Lungimea medie a unei gene este de aproximativ 36.000 de baze

• Lungimea medie a unei gene este de aproximativ 5.200 de baze

• Lungimile minime și maxime ale genelor sunt de aproximativ opt și, respectiv, 2,3 milioane de baze.

Deoarece media este mult mai mare decât mediana, aceasta implică faptul că distribuția lungimii este înclinată spre dreapta. Pentru a avea o privire mai concretă, să diagramăm distribuția.

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

Pandas oferă o interfață simplă pentru matplotlib pentru a face plotarea foarte la îndemână cu DataFrames sau serii.

În acest caz, se spune că vrem o diagramă histogramă ( kind='hist' ) cu 50 de binuri și să lăsăm axa y pe o scară logaritmică ( logy=True ).

Din histogramă, putem vedea că majoritatea genelor se află în primul bin. Cu toate acestea, lungimile unor gene pot fi mai mari de două milioane de baze. Să aflăm care sunt acestea:

 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

După cum puteți vedea, cea mai lungă genă se numește CNTNAP2, care este prescurtarea de la contactin-like protein-like 2. Conform paginii sale wikipedia,

Această genă cuprinde aproape 1,6% din cromozomul 7 și este una dintre cele mai mari gene din genomul uman.

Într-adevăr! Tocmai am verificat noi înșine. În schimb, cum rămâne cu cele mai mici gene? Se dovedește că pot fi de până la opt baze.

 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

Lungimile celor două cazuri extreme sunt între cinci ordine de mărime (2,3 milioane vs. 8), ceea ce este enorm și care poate fi un indiciu al nivelului de diversitate a vieții.

O singură genă poate fi tradusă în multe proteine ​​diferite printr-un proces numit splicing alternativ, ceva ce nu am explorat. Astfel de informații se află și în interiorul fișierului GFF3, dar în afara domeniului acestui post.

Distribuția genelor între cromozomi

Ultimul lucru pe care aș dori să-l discut este distribuția genelor între cromozomi, care servește și ca exemplu pentru introducerea metodei .merge pentru combinarea a două DataFrames. În mod intuitiv, cromozomii mai lungi găzduiesc probabil mai multe gene. Să vedem dacă este adevărat.

 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

Am împrumutat variabila chrs din secțiunea anterioară și am folosit-o pentru a filtra secvențele neasamblate. Pe baza rezultatelor, cel mai mare cromozom 1 are într-adevăr cele mai multe gene. În timp ce cromozomul Y are cel mai mic număr de gene, nu este cel mai mic cromozom.

Rețineți că nu par să existe gene în mitocondrie (MT), ceea ce nu este adevărat.

Un pic mai multă filtrare pe primul DataFrame df returnat de pd.read_csv arată că toate genele MT provin din sursă insdc (care au fost filtrate înainte la generarea edf unde am luat în considerare doar sursele de havana, ensembl sau 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. Aproximativ 0,37% din genomul uman este încă incomplet, chiar dacă prima versiune a apărut acum peste 15 ani.
2. Există aproximativ 42.000 de gene în genomul uman pe baza acestui fișier GFF3 pe care l-am folosit.
3. Lungimea unei gene poate varia de la câteva zeci până la peste două milioane de baze.
4. Genele nu sunt distribuite uniform între cromozomi. În general, cu cât cromozomul este mai mare, cu atât găzduiește mai multe gene, dar pentru un subset de cromozomi, corelația poate fi negativă.

Fișierul GFF3 este foarte bogat în informații de adnotare și tocmai am zgâriat suprafața. Dacă sunteți interesat de explorări suplimentare, iată câteva întrebări cu care vă puteți juca:

1. Câte transcrieri are de obicei o genă? Ce procent de gene au mai mult de o transcriere?
2. Câte izoforme are de obicei o transcriere?
3. Câți exoni, CDS și UTR-uri are de obicei o transcriere? Ce dimensiuni au?
4. Este posibilă clasificarea genelor în funcție de funcția lor, așa cum este descris în coloana de descriere?