Kompleksowe wprowadzenie do Twojego genomu dzięki stosowi SciPy

Bioinformatyka to interdyscyplinarna dziedzina, która opracowuje metody i narzędzia programowe do analizy i zrozumienia danych biologicznych.

Mówiąc prościej, możesz po prostu myśleć o tym jako o nauce o danych dla biologii.

Wśród wielu rodzajów danych biologicznych, dane genomiczne są jednymi z najczęściej analizowanych. Szczególnie w związku z szybkim postępem technologii sekwencjonowania DNA nowej generacji (NGS), ilość danych genomicznych rośnie wykładniczo. Według Stephensa, Zachary'ego D i in., pozyskiwanie danych genomicznych odbywa się w skali eksabajtów rocznie.

W tym poście przedstawiam przykład analizy pliku GFF3 pod kątem ludzkiego genomu za pomocą SciPy Stack. Generic Feature Format Version 3 (GFF3) to obecnie standardowy format pliku tekstowego do przechowywania cech genomowych. W szczególności w tym poście dowiesz się, jak używać stosu SciPy, aby odpowiedzieć na następujące pytania dotyczące ludzkiego genomu:

1. Jaka część genomu jest niekompletna?
2. Ile genów jest w genomie?
3. Jak długi jest typowy gen?
4. Jak wygląda dystrybucja genów wśród chromosomów?

Najnowszy plik GFF3 dotyczący ludzkiego genomu można pobrać stąd. Plik README znajdujący się w tym katalogu zawiera krótki opis tego formatu danych, a dokładniejszą specyfikację można znaleźć tutaj.

Użyjemy Pand, głównego składnika stosu SciPy zapewniającego szybkie, elastyczne i ekspresyjne struktury danych, aby manipulować i zrozumieć plik GFF3.

Organizować coś

Po pierwsze, skonfigurujmy środowisko wirtualne z zainstalowanym stosem SciPy. Ten proces może być czasochłonny, jeśli zostanie zbudowany ręcznie ze źródeł, ponieważ stos obejmuje wiele pakietów – niektóre z nich zależą od zewnętrznego kodu FORTRAN lub C. Tutaj polecam używanie Minicondy, dzięki czemu konfiguracja jest bardzo łatwa.

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

Flaga -b w wierszu bash nakazuje mu wykonanie w trybie wsadowym. Po użyciu powyższych poleceń do pomyślnej instalacji Minicondy, uruchom nowe środowisko wirtualne dla genomiki, a następnie zainstaluj stos SciPy.

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

Zauważ, że określiliśmy tylko 3 pakiety, których będziemy używać w tym poście. Jeśli chcesz, aby wszystkie pakiety zostały wymienione na stosie SciPy, po prostu dołącz je na końcu polecenia conda create .

Jeśli nie masz pewności co do dokładnej nazwy pakietu, spróbuj conda search . Aktywujmy środowisko wirtualne i uruchommy IPythona.

 source activate venv/ ipython

IPython jest znacznie potężniejszym zamiennikiem domyślnego interfejsu interpretera Pythona, więc wszystko, co robiłeś w domyślnym interpreterze Pythona, można również zrobić w IPython. Gorąco polecam każdemu programiście Pythona, który jeszcze nie używał IPythona, aby spróbował.

Pobierz plik adnotacji

Po zakończeniu naszej konfiguracji pobierzmy plik z adnotacjami ludzkiego genomu w formacie GFF3.

Jest to około 37 MB, bardzo mały plik w porównaniu do zawartości informacyjnej ludzkiego genomu, czyli około 3 GB w postaci zwykłego tekstu. Dzieje się tak, ponieważ plik GFF3 zawiera tylko adnotacje sekwencji, podczas gdy dane sekwencji są zwykle przechowywane w innym formacie pliku o nazwie FASTA. Jeśli jesteś zainteresowany, możesz pobrać FASTA tutaj, ale nie będziemy używać danych sekwencji w tym samouczku.

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

Przedrostek ! mówi IPythonowi, że jest to polecenie powłoki zamiast polecenia Pythona. Jednak IPython może również przetwarzać niektóre często używane polecenia powłoki, takie jak ls , pwd , rm , mkdir , rmdir , nawet bez prefiksu ! .

Patrząc na nagłówek pliku GFF, zobaczysz wiele linii metadanych/prag/dyrektyw zaczynających się od ## lub #! .

Według README ## oznacza, że ​​metadane są stabilne, a #! oznacza, że ​​jest eksperymentalny.

Później zobaczysz także ### , która jest kolejną dyrektywą o jeszcze bardziej subtelnym znaczeniu opartym na specyfikacji.

Komentarze czytelne dla ludzi powinny znajdować się po jednym # . Dla uproszczenia potraktujemy wszystkie wiersze zaczynające się od # jako komentarze i po prostu zignorujemy je podczas naszej analizy.

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

Pierwsza linia wskazuje, że wersja formatu GFF użyta w tym pliku to 3.

Poniżej znajdują się podsumowania wszystkich regionów sekwencji. Jak zobaczymy później, takie informacje można znaleźć również w głównej części pliku.

Linie zaczynające się od #! pokazuje informacje o konkretnej budowie genomu, GRCh38.p7, której dotyczy ten plik adnotacji.

Genome Reference Consortium (GCR) to międzynarodowe konsorcjum, które nadzoruje aktualizacje i ulepszenia kilku zestawów genomów referencyjnych, w tym dla człowieka, myszy, danio pręgowanego i kurczaka.

Przeglądając ten plik, oto kilka pierwszych linii adnotacji.

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

Kolumny to seqid , source , type , start , end , score , strand , phase , atrybuty . Niektóre z nich są bardzo łatwe do zrozumienia. Weźmy jako przykład pierwszą linię:

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

Jest to adnotacja pierwszego chromosomu z sekwencją 1, która zaczyna się od pierwszej zasady do 24 895 622 zasady.

Innymi słowy, pierwszy chromosom ma długość około 25 milionów zasad.

Nasza analiza nie będzie wymagała informacji z trzech kolumn o wartości . (tj. wynik, pasmo i faza), więc na razie możemy je po prostu zignorować.

Ostatnia kolumna atrybutów mówi, że Chromosome 1 ma również trzy aliasy, a mianowicie CM000663.2, chr1 i NC_000001.11. W zasadzie tak wygląda plik GFF3, ale nie będziemy sprawdzać ich linia po linii, więc czas załadować cały plik do Pand.

Pandas dobrze nadaje się do radzenia sobie z formatem GFF3, ponieważ jest to plik rozdzielany tabulatorami, a Pandas ma bardzo dobrą obsługę odczytu plików podobnych do CSV.

Zauważ, że jednym wyjątkiem od formatu rozdzielanego tabulatorami jest sytuacja, w której GFF3 zawiera ##FASTA .

Zgodnie ze specyfikacją, ##FASTA wskazuje koniec fragmentu adnotacji, po którym nastąpi jedna lub więcej sekwencji w formacie FASTA (nie rozdzielanym tabulatorami). Ale tak nie jest w przypadku pliku GFF3, który będziemy analizować.

 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)

Ostatni wiersz powyżej ładuje cały plik GFF3 za pomocą metody pandas.read_csv .

Ponieważ nie jest to standardowy plik CSV, musimy nieco dostosować wywołanie.

Najpierw informujemy Pandy o niedostępności informacji nagłówka w GFF3 za pomocą header=None , a następnie podajemy dokładną nazwę dla każdej kolumny za pomocą names=col_names .

Jeśli argument names nie zostanie określony, Pandy użyją liczb przyrostowych jako nazw dla każdej kolumny.

sep='\t' mówi Pandzie, że kolumny są oddzielone tabulatorami zamiast przecinkami. Wartość sep= może być w rzeczywistości wyrażeniem regularnym (regex). Jest to przydatne, jeśli plik pod ręką używa różnych separatorów dla każdej kolumny (o tak, tak się dzieje). comment='#' oznacza, że ​​wiersze zaczynające się od # są uważane za komentarze i będą ignorowane.

compression='gzip' informuje Pandy, że plik wejściowy jest skompresowany w formacie gzip.

Ponadto pandas.read_csv ma ​​bogaty zestaw parametrów, które umożliwiają odczyt różnych formatów plików podobnych do CSV.

Typ zwracanej wartości to DataFrame , która jest najważniejszą strukturą danych w Pandas, używaną do reprezentowania danych 2D.

Pandas ma również strukturę danych Series i Panel odpowiednio dla danych 1D i 3D. Zapoznaj się z dokumentacją, aby uzyskać wprowadzenie do struktur danych Pand.

Przyjrzyjmy się kilku pierwszym wpisom z metodą .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

Dane wyjściowe są ładnie sformatowane w formacie tabelarycznym z dłuższymi ciągami w kolumnie atrybutów częściowo zastąpionymi przez ... .

Możesz ustawić Pandy, aby nie pomijały długich ciągów za pomocą pd.set_option('display.max_colwidth', -1) . Ponadto Pandas ma wiele opcji, które można dostosować.

Następnie zdobądźmy podstawowe informacje o tej ramce danych za pomocą metody .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

To pokazuje, że GFF3 ma 2 601 848 wierszy z adnotacjami, a każdy wiersz ma dziewięć kolumn.

Dla każdej kolumny pokazuje również ich typy danych.

Ten początek i koniec są typu int64, liczb całkowitych reprezentujących pozycje w genomie.

Wszystkie pozostałe kolumny są typu object , co prawdopodobnie oznacza, że ​​ich wartości składają się z kombinacji liczb całkowitych, zmiennoprzecinkowych i ciągów.

Rozmiar wszystkich informacji to około 178,7+ MB przechowywanych w pamięci. Okazuje się, że jest bardziej kompaktowy niż nieskompresowany plik, który będzie miał około 402 MB. Szybka weryfikacja jest pokazana poniżej.

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

Z widoku wysokiego poziomu załadowaliśmy cały plik GFF3 do obiektu DataFrame w Pythonie i cała nasza następująca analiza będzie oparta na tym pojedynczym obiekcie.

Zobaczmy teraz, o co chodzi w pierwszej kolumnie 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 to jeden ze sposobów dostępu do danych kolumn z ramki danych. Innym sposobem jest df['seqid'] , który jest bardziej ogólną składnią, ponieważ jeśli nazwa kolumny jest zastrzeżonym słowem kluczowym Pythona (np. class ) lub zawiera . lub spacji, pierwszy sposób ( df.seqid ) nie zadziała.

Wyniki pokazują, że istnieje 194 unikalnych sekwencji, które obejmują chromosomy od 1 do 22, X, Y i DNA mitochondriów (MT), a także 169 innych sekwencji.

Seqid zaczynające się od KI i GL to sekwencje DNA – lub rusztowania – w genomie, które nie zostały pomyślnie połączone w chromosomy.

Dla tych, którzy nie są zaznajomieni z genomiką, jest to ważne.

Chociaż pierwszy projekt ludzkiego genomu pojawił się ponad 15 lat temu, obecny ludzki genom jest wciąż niekompletny. Trudność w składaniu tych sekwencji jest w dużej mierze spowodowana złożonymi, powtarzającymi się regionami w genomie.

Następnie spójrzmy na kolumnę źródłową.

README mówi, że źródło jest kwalifikatorem dowolnego tekstu, który ma opisywać algorytm lub procedurę operacyjną, która wygenerowała tę funkcję.

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

Jest to przykład użycia metody value_counts , która jest niezwykle przydatna do szybkiego zliczania zmiennych kategorialnych.

Z wyniku widzimy, że istnieje siedem możliwych wartości dla tej kolumny, a większość wpisów w pliku GFF3 pochodzi z havana, ensemble i ensembl_havana.

Możesz dowiedzieć się więcej o znaczeniu tych źródeł i relacjach między nimi w tym poście.

Aby uprościć sprawę, skupimy się na wpisach ze źródeł GRCh38, havana, ensemble i ensembl_havan.a.

Jaka część genomu jest niekompletna?

Informacje o każdym całym chromosomie znajdują się we wpisach ze źródła GRCh38, więc najpierw odfiltrujmy resztę i przypiszmy odfiltrowany wynik do nowej zmiennej 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...

Filtrowanie w Pandas jest łatwe.

Jeśli sprawdzisz wartość obliczoną na podstawie wyrażenia df.source == 'GRCh38' , jest to seria wartości True i False dla każdego wpisu o tym samym indeksie co df . Przekazanie go do df[] zwróci tylko te wpisy, których odpowiadające im wartości są True.

W df[] znajdują się 194 klucze, dla których df.source == 'GRCh38' .

Jak widzieliśmy wcześniej, w kolumnie seqid znajdują się również 194 unikalne wartości, co oznacza, że ​​każdy wpis w gdf odpowiada konkretnemu seqidowi.

Następnie losowo wybieramy 10 wpisów metodą sample , aby przyjrzeć się bliżej.

Widać, że niezmontowane sekwencje są typu superkontig, podczas gdy inne są typu chromosomowego. Aby obliczyć ułamek genomu, który jest niekompletny, najpierw musimy znać długość całego genomu, która jest sumą długości wszystkich sekwid.

 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

W powyższym fragmencie najpierw stworzyliśmy kopię gdf za pomocą .copy() . W przeciwnym razie oryginalny gdf jest tylko kawałkiem df , a jego bezpośrednia modyfikacja spowoduje SettingWithCopyWarning (więcej szczegółów znajdziesz tutaj).

Następnie obliczamy długość każdego wpisu i dodajemy go z powrotem do gdf jako nową kolumnę o nazwie „długość”. Całkowita długość okazuje się wynosić około 3,1 miliarda, a udział niezmontowanych sekwencji wynosi około 0,37%.

Oto jak działa krojenie w dwóch ostatnich poleceniach.

Najpierw tworzymy listę ciągów, która obejmuje wszystkie seqidy dobrze złożonych sekwencji, którymi są wszystkie chromosomy i mitochondria. Następnie używamy metody isin do filtrowania wszystkich wpisów, których seqid znajdują się na liście chrs .

Znak minus ( - ) jest dodawany na początku indeksu, aby odwrócić zaznaczenie, ponieważ tak naprawdę chcemy wszystkiego, czego nie ma na liście (tj. Chcemy, aby niezmontowane zaczynały się od KI i GL)…

Uwaga: Ponieważ złożone i niezmontowane sekwencje są rozróżniane według kolumny typu, ostatni wiersz można alternatywnie przepisać w następujący sposób, aby uzyskać te same wyniki.

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

Ile jest genów?

Tutaj skupiamy się na wpisach ze źródła ensemble, havana i ensembl_havana, ponieważ to tam należy większość wpisów adnotacji.

 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 jest ponownie używana do filtrowania.

Następnie szybkie zliczenie wartości pokazuje, że większość wpisów to ekson, sekwencja kodująca (CDS) i region niepodlegający translacji (UTR).

Są to elementy podgenów, ale głównie szukamy liczby genów. Jak pokazano, jest ich 42 470, ale chcemy wiedzieć więcej.

W szczególności, jakie są ich imiona i co robią? Aby odpowiedzieć na te pytania, musimy uważnie przyjrzeć się informacjom w kolumnie atrybutów.

 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)

Są one sformatowane jako rozdzielona średnikami lista par tag-wartość. Informacje, które nas najbardziej interesują, to nazwa genu, identyfikator i opis genu, które wyodrębnimy za pomocą wyrażenia regularnego (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)

Najpierw wyodrębniamy nazwy genów.

W wyrażeniu regularnym Name=(?P<gene_name>.+?); , +? jest używany zamiast + , ponieważ chcemy, aby nie był zachłanny i aby wyszukiwanie zatrzymywało się na pierwszym średniku; w przeciwnym razie wynik zostanie dopasowany do ostatniego średnika.

Ponadto wyrażenie regularne jest najpierw kompilowane za pomocą re.compile , zamiast być używane bezpośrednio, jak w przypadku re.search , aby uzyskać lepszą wydajność, ponieważ zastosujemy je do tysięcy ciągów atrybutów.

extract_gene_name służy jako funkcja pomocnicza do użycia w pd.apply , która jest metodą używaną, gdy funkcja musi być zastosowana do każdego wpisu ramki danych lub serii.

W tym konkretnym przypadku chcemy wyodrębnić nazwę genu dla każdego wpisu w ndf.attributes i dodać nazwy z powrotem do ndf w nowej kolumnie nazwanej gene_name .

Identyfikatory i opis genów są wyodrębniane w podobny sposób.

 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)

Wyrażenie regularne dla RE_GENE_ID jest nieco bardziej szczegółowe, ponieważ wiemy, że każdy gene_id musi zaczynać się od ENSG , gdzie ENS oznacza zespół, a G oznacza gen.

W przypadku wpisów, które nie mają żadnego opisu, zwrócimy pusty ciąg. Po wypakowaniu wszystkiego nie będziemy już używać kolumny atrybutów, więc usuńmy ją, aby zachować porządek za pomocą metody .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...

W powyższym wywołaniu attributes wskazują konkretną kolumnę, którą chcemy usunąć.

axis=1 oznacza, że ​​upuszczamy kolumnę zamiast wiersza (domyślnie axis=0 ).

inplace=True oznacza, że ​​upuszczanie jest wykonywane na samej ramce DataFrame zamiast zwracania nowej kopii z upuszczoną określoną kolumną.

Szybki wygląd .head pokazuje, że kolumna atrybutów rzeczywiście desc i dodano trzy nowe kolumny: gene_name , identyfikator_genu i gene_id .

Z ciekawości zobaczmy, czy wszystkie gene_name gene_id unikalne:

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

Co zaskakujące, liczba nazw genów jest mniejsza niż identyfikatorów genów, co wskazuje, że jakaś nazwa_genu musi odpowiadać wielu identyfikatorom genów. Dowiedzmy się, jakie one są.

 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

Grupujemy wszystkie wpisy według wartości gene_name , a następnie zliczamy liczbę elementów w każdej grupie za pomocą .count() .

Jeśli sprawdzisz dane wyjściowe z ndf.groupby('gene_name').count() , wszystkie kolumny są liczone dla każdej grupy, ale większość z nich ma te same wartości.

Zwróć uwagę, że wartości NA nie będą brane pod uwagę podczas liczenia, więc weź tylko liczbę z pierwszej kolumny, seqid ( używamy .ix[:, 0] , aby upewnić się, że nie ma wartości NA).

Następnie posortuj wartości licznika za pomocą .sort_values ​​i odwróć kolejność za pomocą .ix[::-1] .

W rezultacie nazwa genu może być współdzielona z maksymalnie siedmioma identyfikatorami genów.

 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]

Bliższe spojrzenie na wszystkie geny SCARNA20 pokazuje, że rzeczywiście wszystkie są różne.

Chociaż mają tę samą nazwę, znajdują się w różnych pozycjach genomu.

Ich opisy nie wydają się jednak zbyt pomocne w ich rozróżnianiu.

Chodzi o to, aby wiedzieć, że nazwy genów nie są unikalne dla wszystkich identyfikatorów genów i około 0,15% nazw wspólnych dla wielu genów.

Jak długi jest typowy gen?

Podobnie do tego, co zrobiliśmy, gdy próbowaliśmy zrozumieć niekompletność genomu, możemy łatwo dodać kolumnę length do 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() oblicza kilka prostych statystyk na podstawie wartości długości:

• Średnia długość genu wynosi około 36 000 zasad

• Mediana długości genu wynosi około 5200 zasad

• Minimalna i maksymalna długość genów to odpowiednio około 8 i 2,3 miliona zasad.

Ponieważ średnia jest znacznie większa niż mediana, oznacza to, że rozkład długości jest przechylony w prawo. Aby uzyskać bardziej konkretny wygląd, narysujmy rozkład.

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

Pandas zapewnia prosty interfejs do matplotlib, dzięki czemu kreślenie jest bardzo przydatne za pomocą DataFrames lub serii.

W tym przypadku mówi, że chcemy wykresu histogramu ( kind='hist' ) z 50 przedziałami i niech oś y będzie na skali logarytmicznej ( logy=True ).

Z histogramu widać, że większość genów znajduje się w pierwszym przedziale. Jednak niektóre długości genów mogą przekraczać dwa miliony zasad. Dowiedzmy się, jakie one są:

 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

Jak widać, najdłuższy gen nazywa się CNTNAP2, co jest skrótem od białka związanego z kontaktyną 2. Według strony wikipedii,

Gen ten obejmuje prawie 1,6% chromosomu 7 i jest jednym z największych genów w ludzkim genomie.

Rzeczywiście! Po prostu sami to sprawdziliśmy. A co z najmniejszymi genami? Okazuje się, że mogą mieć nawet osiem podstawek.

 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

Długości dwóch skrajnych przypadków są oddalone od siebie o pięć rzędów wielkości (2,3 miliona w porównaniu z 8), co jest ogromne i może wskazywać na poziom różnorodności życia.

Pojedynczy gen może zostać przetłumaczony na wiele różnych białek w procesie zwanym alternatywnym splicingiem, czego jeszcze nie zbadaliśmy. Taka informacja znajduje się również w pliku GFF3, ale poza zakresem tego postu.

Dystrybucja genów wśród chromosomów

Ostatnią rzeczą, którą chciałbym omówić, jest dystrybucja genów między chromosomami, która służy również jako przykład wprowadzenia metody .merge do łączenia dwóch ramek DataFrame. Intuicyjnie, dłuższe chromosomy prawdopodobnie zawierają więcej genów. Zobaczmy, czy to prawda.

 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

Pożyczyliśmy zmienną chrs z poprzedniej sekcji i użyliśmy jej do odfiltrowania niezłożonych sekwencji. Na podstawie wyników można stwierdzić, że największy Chromosom 1 rzeczywiście ma najwięcej genów. Chociaż chromosom Y ma najmniejszą liczbę genów, nie jest najmniejszym chromosomem.

Zauważ, że wydaje się, że w mitochondrium (MT) nie ma genów, co nie jest prawdą.

Nieco więcej filtrowania na pierwszym df DataFrame zwróconym przez pd.read_csv pokazuje, że wszystkie geny MT pochodzą ze źródła insdc (które zostały wcześniej odfiltrowane podczas generowania edf , gdzie braliśmy pod uwagę tylko źródła havana, ensemble lub 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. Około 0,37% ludzkiego genomu jest nadal niekompletne, mimo że pierwszy szkic pojawił się ponad 15 lat temu.
2. W ludzkim genomie jest około 42 000 genów opartych na tym konkretnym pliku GFF3, którego użyliśmy.
3. Długość genu może wynosić od kilkudziesięciu do ponad dwóch milionów zasad.
4. Geny nie są równomiernie rozmieszczone w chromosomach. Ogólnie rzecz biorąc, im większy chromosom, tym więcej genów zawiera, ale w przypadku podzbioru chromosomów korelacja może być ujemna.

Plik GFF3 jest bardzo bogaty w informacje o adnotacjach, a my właśnie zarysowaliśmy powierzchnię. Jeśli jesteś zainteresowany dalszą eksploracją, oto kilka pytań, z którymi możesz się bawić:

1. Ile transkryptów zazwyczaj ma gen? Jaki procent genów ma więcej niż 1 transkrypt?
2. Ile izoform ma zazwyczaj transkrypt?
3. Ile eksonów, CDS i UTR zazwyczaj zawiera transkrypcja? Jakie są rozmiary?
4. Czy można kategoryzować geny na podstawie ich funkcji, jak opisano w kolumnie opisu?