Bioinformatikk er et tverrfaglig felt som utvikler metoder og programvareverktøy for å analysere og forstå biologiske data.

Mer enkelt sagt, du kan bare tenke på det som datavitenskap for biologi.

Blant de mange typene biologiske data er genomikkdata en av de mest analyserte. Spesielt med den raske utviklingen av neste generasjons DNA-sekvenseringsteknologi (NGS), har volumet av genomikkdata vokst eksponentielt. I følge Stephens, Zachary D et al. , er datainnsamling av genomikk i skalaen exabyte per år.

SciPy samler mange Python-moduler for vitenskapelig databehandling, noe som er ideelt for mange bioinformatiske behov. kvitring

I dette innlegget demonstrerer jeg et eksempel på å analysere en GFF3-fil for det menneskelige genomet med SciPy Stack. Generisk funksjonsformat versjon 3 (GFF3) er det nåværende standard tekstfilformatet for lagring av genomiske funksjoner. Spesielt i dette innlegget vil du lære hvordan du bruker SciPy-stakken for å svare på følgende spørsmål om det menneskelige genomet:

1. Hvor mye av genomet er ufullstendig?
2. Hvor mange gener er det i genomet?
3. Hvor lenge er et typisk gen?
4. Hvordan ser genfordeling blant kromosomer ut?

Den siste GFF3-filen for det menneskelige genomet kan lastes ned fra her . De LES LES filen som kommer i denne katalogen gir en kort beskrivelse av dette dataformatet, og en grundigere spesifikasjon er funnet her .

Vi vil bruke Pandaer , en hovedkomponent av SciPy-stakken som gir raske, fleksible og uttrykksfulle datastrukturer, for å manipulere og forstå GFF3-filen.

Oppsett

Først og fremst, la oss sette opp et virtuelt miljø med SciPy-stakken installert. Denne prosessen kan være tidkrevende hvis den bygges fra kilden manuelt, da stakken involverer mange pakker - hvorav noen avhenger av ekstern FORTRAN- eller C-kode. Her anbefaler jeg å bruke Miniconda , noe som gjør oppsettet veldig enkelt.

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

-b flagg på bash-linjen forteller at den skal kjøres i batch-modus. Etter at kommandoene ovenfor er brukt til å installere Miniconda, starter du et nytt virtuelt miljø for genomikk og installerer deretter SciPy-stakken.

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

Merk at vi bare har spesifisert de 3 pakkene vi skal bruke i dette innlegget. Hvis du vil at alle pakkene skal være oppført i SciPy-stakken, kan du bare legge dem til på slutten av conda create kommando.

Hvis du er usikker på det eksakte navnet på en pakke, kan du prøve conda search. La oss aktivere det virtuelle miljøet og starte IPython.

source activate venv/ ipython

IPython er en betydelig kraftigere erstatning til standard Python-tolkegrensesnitt , så hva du pleide å gjøre i standard pythontolker kan også gjøres i IPython. Jeg anbefaler alle Python-programmerere, som ikke har brukt IPython ennå, å prøve.

Last ned kommentarfilen

Når oppsettet vårt nå er fullført, la oss laste ned filen for menneskelig genomkommentar i GFF3-format.

Det handler om 37 MB, en veldig liten fil sammenlignet med informasjonsinnholdet i et menneskelig genom, som er omtrent 3 GB i ren tekst. Det er fordi GFF3-filen bare inneholder merknaden av sekvensene, mens sekvensdataene vanligvis lagres i et annet filformat kalt FASTA . Hvis du er interessert, kan du laste ned FASTA her , men vi bruker ikke sekvensdataene i denne opplæringen.

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

Det prefikset ! forteller IPython at dette er en skallkommando i stedet for en Python-kommando. Imidlertid kan IPython også behandle noen ofte brukte skallkommandoer som ls, pwd, rm, mkdir, rmdir selv uten et prefiks !.

Når du ser på hodet til GFF-filen, vil du se mange metadata / pragmas / direktelinjer som begynner med ## eller #!.

Ifølge LES LES , ## betyr at metadataene er stabile, mens #! betyr at det er eksperimentelt.

Senere vil du også se ###, som er et annet direktiv med enda mer subtil betydning basert på spesifikasjon .

Menneskelige lesbare kommentarer skal være etter en enkelt #. For enkelhets skyld vil vi behandle alle linjene som begynner med # som kommentarer, og bare ignorere dem under analysen.

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

Den første linjen indikerer at versjonen av GFF-format som brukes i denne filen er 3.

Deretter følger oppsummeringer av alle sekvensregioner. Som vi vil se senere, kan slik informasjon også bli funnet i hoveddelen av filen.

Linjene som begynner med #! viser informasjon om den spesielle bygningen av genomet, GRCh38.p7, som denne kommentarfilen gjelder for.

Genomreferansekonsortium (GCR) er et internasjonalt konsortium som fører tilsyn med oppdateringer og forbedringer av flere referansegenomsamlinger, inkludert de for mennesker, mus, sebrafisk og kylling.

Skanning gjennom denne filen, her er de første få kommentarlinjene.

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

Kolonnene er seqid , kilde , type , start , slutt , score , Strand , fase , attributter . Noen av dem er veldig enkle å forstå. Ta den første linjen som et eksempel:

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

Dette er merknaden av det første kromosomet med en seqid på 1, som starter fra den første basen til den 24 895 622. basen.

Med andre ord, det første kromosomet er omtrent 25 millioner baser langt.

Analysen vår trenger ikke informasjon fra de tre kolonnene med verdien . (dvs. poengsum, streng og fase), så vi kan bare ignorere dem for nå.

Den siste attributtkolonnen sier at Kromosom 1 også har tre aliasnavn, nemlig CM000663.2, chr1 og NC_000001.11. Slik ser en GFF3-fil ut, men vi vil ikke inspisere dem linje for linje, så det er på tide å laste hele filen inn i Pandas.

Pandas passer godt for å håndtere GFF3-format fordi det er en tabulatoravgrenset fil, og Pandas har veldig god støtte for å lese CSV-lignende filer.

Merk at et unntak fra det tabulatoravgrensede formatet er når GFF3 inneholder ##FASTA .

Ifølge spesifikasjon , ##FASTA angir slutten på en kommenteringsdel, som vil bli fulgt med en eller flere sekvenser i FASTA (et ikke-tabulatoravgrenset) format. Men dette er ikke tilfelle for GFF3-filen vi skal analysere.

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)

Den siste linjen over laster hele GFF3-filen med pandas.read_csv metode.

Siden det ikke er en standard CSV-fil, må vi tilpasse samtalen litt.

Først informerer vi Pandas om utilgjengeligheten av topptekstinformasjon i GFF3 med header=None, og deretter spesifiserer vi det nøyaktige navnet for hver kolonne med names=col_names.

Hvis names argument er ikke spesifisert, vil Pandas bruke trinnnummer som navn for hver kolonne.

sep=' ' forteller Panda at kolonnene er tabu-skilt i stedet for komma-skilt. Verdien til sep= kan faktisk være et vanlig uttrykk (regex). Dette blir nyttig hvis filen for hånden bruker forskjellige separatorer for hver kolonne (åh, det skjer). comment='#' betyr linjer som starter med # regnes som kommentarer og vil bli ignorert.

compression='gzip' forteller Pandas at inndatafilen er en gzip-komprimert.

I tillegg, pandas.read_csv har et rikt sett med parametere som gjør det mulig å lese forskjellige typer CSV-lignende filformater.

Typen av returnert verdi er DataFrame, som er den viktigste datastrukturen i Pandas, brukt til å representere 2D-data.

Pandas har også en Series og Panel datastruktur for henholdsvis 1D og 3D-data. Vennligst referer til dokumentasjon for en introduksjon til Pandas ’datastrukturer.

La oss ta en titt på de første få oppføringene med .head metode.

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

Utdataene er pent formatert i tabellformat med lengre strenger i attributtkolonnen delvis erstattet med ....

Du kan angi at Pandaer ikke skal utelate lange strenger med pd.set_option('display.max_colwidth', -1). I tillegg har Pandas mange alternativer som kan tilpasses.

La oss deretter få grunnleggende informasjon om denne datarammen med .info metode.

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

Dette viser at GFF3 har 2 601 848 merkede linjer, og hver linje har ni kolonner.

For hver kolonne viser den også datatypene.

Start og slutt er av int64-type, heltall som representerer posisjoner i genomet.

De andre kolonnene er alle av typen object, noe som sannsynligvis betyr at deres verdier består av en blanding av heltall, flyter og strenger.

Størrelsen på all informasjon er omtrent 178,7+ MB lagret i minnet. Dette viser seg å være mer kompakt enn den ukomprimerte filen, som vil være omtrent 402 MB. En rask bekreftelse er vist nedenfor.

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

Fra et høyt nivå har vi lastet hele GFF3-filen inn i et DataFrame-objekt i Python, og alle våre følgende analyser vil være basert på dette enkelt objektet.

La oss nå se hva den første kolonnen seqid handler om.

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 er en måte å få tilgang til kolonnedata fra en dataramme. En annen måte er df['seqid'], som er mer generell syntaks, for hvis kolonnenavnet er et Python-reservert nøkkelord (f.eks. class) Eller inneholder et . eller mellomrom, fungerer den første måten (df.seqid) ikke.

Resultatet viser at det er 194 unike seqider, som inkluderer kromosomer 1 til 22, X, Y og mitokondrion (MT) DNA, samt 169 andre seqider.

Seqidene som begynner med KI og GL er DNA-sekvenser - eller stillaser - i genomet som ikke har blitt vellykket samlet i kromosomene.

For de som ikke er kjent med genomikk, er dette viktig.

Selv om det første menneskelige genomutkastet kom ut for mer enn 15 år siden, er det nåværende menneskelige genomet fortsatt ufullstendig. Vanskeligheten med å samle disse sekvensene skyldes i stor grad komplekse repeterende regioner i genomet .

La oss ta en titt på kildekolonnen.

De LES LES sier at kilden er en fritekstkvalifiseringsprogram beregnet på å beskrive algoritmen eller driftsprosedyren som genererte denne funksjonen.

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

Dette er et eksempel på bruk av value_counts metode, som er ekstremt nyttig for en rask telling av kategoriske variabler.

Fra resultatet kan vi se at det er syv mulige verdier for denne kolonnen, og flertallet av oppføringene i GFF3-filen kommer fra havana, ensembl og ensembl_havana.

Du kan lære mer om hva disse kildene betyr og forholdet mellom dem i denne posten .

For å holde ting enkelt, vil vi fokusere på oppføringer fra kilder GRCh38, havana, ensembl og ensembl_havan.a.

Hvor mye av genomet er ufullstendig?

Informasjonen om hvert hele kromosom er i oppføringene fra kilden GRCh38, så la oss først filtrere ut resten, og tildele det filtrerte resultatet til en ny variabel 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...

Det er enkelt å filtrere i Pandas.

Hvis du inspiserer verdien evaluert fra uttrykket df.source == 'GRCh38', er det en serie med True og False verdier for hver oppføring med samme indeks som df. Gi den videre til df[] vil bare returnere de oppføringene der de tilsvarende verdiene er sanne.

Det er 194 nøkler i df[] for hvilken df.source == 'GRCh38'.

Som vi har sett tidligere, er det også 194 unike verdier i seqid kolonne, som betyr hver oppføring i gdf tilsvarer en bestemt seqid.

Deretter velger vi tilfeldig 10 oppføringer med sample metode for å se nærmere på.

Du kan se at de umonterte sekvensene er av typen superkontig mens de andre er av kromosom. For å beregne brøkdelen av genomet som er ufullstendig, må vi først vite lengden på hele genomet, som er summen av lengden på alle sekvider.

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

I utdraget ovenfor laget vi en kopi av gdf med .copy(). Ellers er originalen gdf er bare et stykke av df, og endring av det direkte vil resultere i SettingWithCopyWarning (se her for flere detaljer).

Vi beregner deretter lengden på hver oppføring og legger den tilbake til gdf som en ny kolonne kalt “lengde”. Den totale lengden viser seg å være omtrent 3,1 milliarder kroner, og brøkdelen av umonterte sekvenser er omtrent 0,37%.

Slik fungerer oppskjæringen i de to siste kommandoene.

Først oppretter vi en liste over strenger som dekker alle sekvider av godt sammensatte sekvenser, som alle er kromosomer og mitokondrier. Vi bruker deretter isin metode for å filtrere alle oppføringer hvis seqid er i chrs liste.

Et minustegn (-) legges til begynnelsen av indeksen for å reversere valget, fordi vi faktisk vil ha alt som er ikke i listen (dvs. vi vil ha de umonterte som starter med KI og GL) ...

Merk: Siden de sammensatte og umonterte sekvensene er preget av typekolonnen, kan den siste linjen alternativt skrives om som følger for å oppnå de samme resultatene.

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

Hvor mange gener er det?

Her fokuserer vi på oppføringene fra kildeensembl, havana og ensembl_havana siden de er der flertallet av kommentaroppføringene hører hjemme.

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

isin metoden brukes igjen for filtrering.

Deretter viser en rask verditelling at flertallet av oppføringene er exon, kodende sekvens (CDS) og ikke-oversatt region (UTR).

Dette er undergenelementer, men vi ser hovedsakelig etter genantallet. Som vist er det 42.470, men vi vil vite mer.

Spesielt hva heter de, og hva gjør de? For å svare på disse spørsmålene, må vi se nøye på informasjonen i attributtkolonnen.

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)

De er formatert som en semikolonseparert liste over tag-verdi-par. Informasjonen vi er mest interessert i er gennavn, gen-ID og beskrivelse, og vi vil trekke dem ut med vanlig uttrykk (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)

Først trekker vi ut gennavnene.

I regex Name=(?P.+?); , +? brukes i stedet for + fordi vi vil at den skal være ikke grådig og la søket stoppe ved første semikolon; Ellers vil resultatet matche opp til det siste semikolonet.

Regex blir også først kompilert med re.compile i stedet for å bli brukt direkte som i re.search for bedre ytelse fordi vi vil bruke den på tusenvis av attributtstrenger.

extract_gene_name fungerer som en hjelperfunksjon som skal brukes i pd.apply, som er metoden som skal brukes når en funksjon må brukes på hver oppføring i en dataramme eller serie.

I dette spesielle tilfellet ønsker vi å trekke ut gennavnet for hver oppføring i ndf.attributes, og legge navnene tilbake til ndf i en ny kolonne kalt gene_name.

Gen-ID og beskrivelse trekkes ut på en lignende måte.

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)

Regex for RE_GENE_ID er litt mer spesifikk siden vi vet at hver gene_id må starte med ENSG, der ENS midler ensembl og G betyr gen.

For oppføringer som ikke har noen beskrivelse, returnerer vi en tom streng. Etter at alt er hentet ut, bruker vi ikke attributtkolonnen lenger, så la oss slippe den for å holde ting pent og rent med metoden .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...

I ovennevnte samtale, attributes indikerer den spesifikke kolonnen vi vil slippe.

axis=1 betyr at vi slipper en kolonne i stedet for en rad (axis=0 som standard).

inplace=True betyr at slipp opereres på selve DataFrame i stedet for å returnere en ny kopi med spesifisert kolonne falt.

En rask .head utseende viser at attributtkolonnen faktisk er borte, og tre nye kolonner: gene_name, gene_id og desc er lagt til.

Av nysgjerrighet, la oss se om alle gene_id og gene_name er unike:

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

Overraskende nok er antallet gennavn mindre enn for gen-ID-er, noe som indikerer at noe gennavn må tilsvare flere gen-ID-er. La oss finne ut hva de er.

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

Vi grupperer alle oppføringene etter verdien gene_name, og teller deretter antall elementer i hver gruppe med .count().

Hvis du inspiserer utdataene fra ndf.groupby('gene_name').count(), telles alle kolonnene for hver gruppe, men de fleste av dem har de samme verdiene.

Vær oppmerksom på at NA-verdier ikke blir vurdert når du teller, så ta bare tellingen av den første kolonnen seqid (vi bruker .ix[:, 0] for å sikre at det ikke er noen NA-verdier).

Sorter deretter telleverdiene med .sort_values og snu rekkefølgen med .ix[::-1].

I resultatet kan et gennavn deles med opptil syv gen-IDer.

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]

En nærmere titt på alle SCARNA20-genene viser at de virkelig er forskjellige.

Mens de har samme navn, ligger de på forskjellige posisjoner i genomet.

Deres beskrivelser virker imidlertid ikke veldig nyttige for å skille dem ut.

Poenget her er å vite at gennavn ikke er unikt for alle gen-IDer, og omtrent 0,15% av navnene som deles av flere gener.

Hvor lang er et typisk gen?

I likhet med hva vi gjorde da vi prøvde å forstå genomets ufullstendighet, kan vi enkelt legge til et length kolonne til 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() beregner enkel statistikk basert på lengdeverdiene:

• Gjennomsnittlig lengde på et gen er ca 36.000 baser

• Medianlengden på et gen er omtrent 5 200 baser lang

• Minimum og maksimal genlengde er henholdsvis om lag åtte og 2,3 millioner baser.

Fordi gjennomsnittet er mye større enn medianen, innebærer det at lengdefordelingen er skjev til høyre. For å få et mer konkret utseende, la oss plotte fordelingen.

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

Pandas gir et enkelt grensesnitt til matplotlib for å gjøre plotting veldig nyttig med DataFrames eller serier.

I dette tilfellet står det at vi vil ha et histogramdiagram (kind='hist') med 50 søppel, og la y-aksen være på en loggskala (logy=True).

Fra histogrammet kan vi se at flertallet av gener er innenfor den første bin. Imidlertid kan noen genlengder være mer enn to millioner baser. La oss finne ut hva de er:

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

Som du kan se, heter det lengste genet CNTNAP2, som er en forkortelse for kontaktin-assosiert proteinlignende 2. I henhold til dets wikipedia-side ,

Dette genet omfatter nesten 1,6% av kromosom 7 og er et av de største genene i det menneskelige genomet.

Faktisk! Vi har bare bekreftet det selv. I kontrast, hva med de minste genene? Det viser seg at de kan være så korte som åtte baser.

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

Lengden på de to ekstreme tilfellene er fem størrelsesordener fra hverandre (2,3 millioner mot 8), noe som er enormt og som kan være en indikasjon på nivået på livets mangfold.

Et enkelt gen kan oversettes til mange forskjellige proteiner via en prosess som kalles alternativ spleising, noe vi ikke har utforsket. Slik informasjon er også inne i GFF3-filen, men utenfor omfanget av dette innlegget.

Genfordeling blant kromosomer

Det siste jeg vil diskutere er genfordeling blant kromosomer, som også fungerer som et eksempel for innføring av .merge metode for å kombinere to DataFrames. Intuitivt er det lengre kromosomer som sannsynligvis er vert for flere gener. La oss se om det er sant.

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

Vi lånte chrs variabel fra forrige seksjon, og brukte den til å filtrere ut de umonterte sekvensene. Basert på produksjonen har faktisk det største kromosom 1 flest gener. Mens kromosom Y har det minste antallet gener, er det ikke det minste kromosomet.

Merk at det ser ut til å være ingen gener i mitokondrion (MT), noe som ikke er sant.

Litt mer filtrering på den første DataFrame df returnert av pd.read_csv viser at alle MT-gener er fra kilde insdc (som ble filtrert ut før når vi genererte edf der vi bare vurderte kilder til havana, ensembl eller 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...

Dette eksemplet viser også hvordan man kombinerer to forhold under filtrering med &; den logiske operatoren for “eller” ville være |.

Merk at parentesene rundt hver tilstand er påkrevd, og denne delen av syntaksen i Pandas er forskjellig fra Python, som ville ha vært bokstavelig and og or.

La oss deretter låne gdf DataFrame fra forrige avsnitt som kilde for lengden på hvert kromosom:

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

Kolonnene som ikke er relevante for analysen, slippes for klarhetens skyld.

Ja, .drop kan også ta en liste over kolonner og slippe dem helt i en operasjon.

Merk at raden med en seqid av MT fremdeles er der; vi kommer tilbake til det senere. Den neste operasjonen vi skal utføre er å slå sammen de to datasettene basert på verdiene til 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

Siden chr_gene_counts fortsatt er et serieobjekt, som ikke støtter en sammenslåing, må vi først konvertere det til et DataFrame-objekt med .to_frame

.reset_index() konverterer den opprinnelige indeksen (dvs. seqid) til en ny kolonne og tilbakestiller gjeldende indeks som 0-baserte trinn.

Resultatet fra cdf.head(2) viser hvordan det ser ut. Deretter brukte vi .merge metode for å kombinere de to DataFrame i seqid-kolonnen (on='seqid').

Etter sammenslåing gdf og cdf, MT oppføringen mangler fortsatt. Dette skyldes at .merge som standard driver en indre skjøt, mens venstre sammenføyning, høyre sammenføyning eller ytre sammenføyning er tilgjengelig ved å stille inn how parameter.

Vennligst referer til dokumentasjon for flere detaljer.

Senere kan du oppdage at det også er et beslektet .join metode. .merge og .join er like, men har forskjellige API-er.

I følge tjenestemannen dokumentasjon sier

Den relaterte DataFrame.join-metoden bruker fusjon internt for indeks-på-indeks og indeks-på-kolonne (r), men blir med på indekser som standard i stedet for å prøve å bli med på vanlige kolonner (standardadferd for sammenslåing). Hvis du blir med på indeksen, kan det være lurt å bruke DataFrame.join for å spare deg for å skrive.

I utgangspunktet, .merge er mer generelt formål og brukes av .join.

Til slutt er vi klare til å beregne sammenhengen mellom kromosom length og 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

Som standard .corr beregner Pearson-korrelasjon mellom alle kolonnepar i en dataramme.

Men vi har bare et enkelt par kolonner i dette tilfellet, og korrelasjonen viser seg å være positiv - 0,73.

Med andre ord, jo større kromosom, desto mer sannsynlig er det å ha flere gener. La oss også plotte de to kolonnene etter at vi har sortert verdiparene etter 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))

Som vist på bildet ovenfor, selv om det generelt er en positiv korrelasjon, holder det ikke for alle kromosomer. Spesielt for kromosom 17, 16, 15, 14, 13 er korrelasjonen faktisk negativ, noe som betyr at antallet gener på kromosomet avtar etter hvert som kromosomstørrelsen øker.

Funn og fremtidig forskning

Det avslutter opplæringen vår om manipulering av en kommentarfil for menneskelig genom i GFF3-format med SciPy-stakken. Verktøyene vi hovedsakelig har brukt inkluderer IPython, Pandas og matplotlib. I løpet av opplæringen har vi ikke bare lært noen av de vanligste og mest nyttige operasjonene i Pandaer, vi har også svart på noen veldig interessante spørsmål om genomet vårt. Oppsummert:

1. Omtrent 0,37% av det menneskelige genomet er fremdeles ufullstendig, selv om det første utkastet kom ut for over 15 år siden.
2. Det er omtrent 42.000 gener i det menneskelige genomet basert på dette GFF3-fil vi brukte.
3. Lengden på et gen kan variere fra noen titalls til over to millioner baser.
4. Gener er ikke jevnt fordelt mellom kromosomene. Samlet sett er jo større kromosom, jo ​​flere gener er det vert, men for en delmengde av kromosomene, kan korrelasjonen være negativ.

GFF3-filen er veldig rik på merknadsinformasjon, og vi har nettopp skrapet overflaten. Hvis du er interessert i videre utforskning, er det noen spørsmål du kan leke med:

1. Hvor mange transkripsjoner har et gen vanligvis? Hvor stor prosentandel av gener har mer enn 1 transkripsjon?
2. Hvor mange isoformer har en transkripsjon vanligvis?
3. Hvor mange eksoner, CDS og UTR-er har en transkripsjon vanligvis? Hvilke størrelser er de?
4. Er det mulig å kategorisere genene basert på deres funksjon som beskrevet i beskrivelseskolonnen?