Прогнозиране на протеинови функции с помощта на машинно обучение

Част I: Как да създадете референтни показатели и основни белтъчни предиктори, като използвате методологията за състезание CAFA2

TL;DR: Отидете в раздела Резултати, за да видите ефективността на базовите предсказатели на белтъчна анотация, като използвате набора от показатели, изграден в раздела Методи.

Въведение

Това е първата от поредица статии от няколко части, фокусирани върху разработването и анализа на ефективността на методи за прогнозиране на функцията на протеините. Серията ще бъде организирана, както следва:

• Част I: Тази статия ще предостави ръководство стъпка по стъпка за това как да създадете и оцените модел спрямо набор от показатели, използвайки същата методология като тази, използвана в състезанието CAFA2. Стъпките за постигане на тази цел са:

1. Създайте файлове с анотации за обучение, които ще служат като набор за обучение за предикторите.
2. Генерирайте поредица от файлове с анотации за бенчмарк, които ще служат като набор от бенчмаркове за оценка на предикторите.
3. Обучете два предиктора: базиран на BLAST и наивен; тези предиктори бяха използвани като базова линия в състезанието CAFA2, за да се сравнят с моделите на участниците.
4. Оценете предикторите, базирани на BLAST и Naive, спрямо бенчмарка, създаден в стъпка 2, и изчислете тяхната ефективност.

• Част II и III: В следващите статии ще разработя, обуча и оценя модели, базирани на задълбочено обучение, и ще сравня тяхното представяне с предикторите на базовата линия, както и с моделите, изпратени до състезанието CAFA2.

Наличност на изходния код

Целият код и файлове, използвани в тази статия, в момента са в частни хранилища в GitHub. Ако се интересувате от този материал и разглеждате по-задълбочено кода, не се колебайте да се свържете с мен в GitHub или LinkedIn. Ще се радвам да се свържа и да ви дам достъп до изходния код.

История на състезанието CAFA

Предизвикайте логистиката

Състезанието CAFA е базирано на времевата линия предизвикателство, организирано както следва:

• В t_₁ предизвикателството е обявено и организаторите на CAFA публикуват списъка с целеви последователности. Например CAFA2 беше обявен през септември 2013 г.
• t₀ е крайният срок за състезаващите се отбори да представят функционалните прогнози на целевите протеини. Крайният срок за CAFA2 беше януари 2014 г.
• t₀ бележи началото на период на „растеж на анотациите“, при който се натрупват нови експериментални анотации върху целевите протеини. Този период продължава около десет месеца.
• В t₁ (октомври 2014 г. за CAFA2) се събират целеви последователности с техните нови експериментални анотации и се създават набори от показатели.

Следващата фигура показва графика на състезанието CAFA2:

анотации

Състезанията CAFA използват Генната онтология (GO) като механизъм за анотиране на протеини. Онтологията е набор от контролиран речник, който се използва за точно дефиниране на въпросния термин. GO, по-специално, дефинира термини за три различни аспекта, свързани с протеини (известни още като пространства от имена): биологичният процес (BP), към който е свързан протеинът, клетъчният компонент (CC) и молекулярната функция (MF), с които е свързан или част на. GO е „винаги развиващ се“ набор, тъй като нови термини се включват, усъвършенстват или стават остарели.

Термините в онтологията на GO са организирани в дървовидна структура, където всеки от тях получава уникален идентификатор. Дървовидната структура на GO позволява свързване на термини, като се започне от по-общите на горните нива на дървото на GO до по-специфичните, разположени по-дълбоко в дървото. Например терминът на GO „GO:0032259“ представлява процеса на „метилиране“, който е вид метаболитен процес. Следната фигура, получена от уеб базирания инструмент QuickGO показва представянето на GO:0032259 в дървото на GO:

В резултат на изследователски експерименти учените могат да открият нови биологични процеси, места или функции, които протеинът изпълнява. Когато това се случи, протеините се анотират със съответния GO термин и тези анотации след това се включват в много бази данни за протеини като UniProtKB и UniProt-GOA. Например, следният фрагмент е записът за протеина „Putative methyltransferase 235L“ (ID: 235L_IIV6) в базата данни UniProtKB/Swiss-Prot (издание през април 2021 г.), анотиран с термина „GO:0032259“, споменат по-горе:

ID   235L_IIV6               Reviewed;         265 AA.
AC   Q91FT7;
DT   16-JUN-2009, integrated into UniProtKB/Swiss-Prot.
...
DR   GO; GO:0008168; F:methyltransferase activity; IEA:UniProtKB-KW.
DR   GO; GO:0032259; P:methylation; IEA:UniProtKB-KW.
...

Набор от показатели

Бенчмаркът е набор от протеини, които придобиват нови експериментални анотации по време на периода на „нарастване на анотациите“ (т.е. между времената t₀ и t₁). С други думи, бенчмаркът се генерира чрез записване на анотациите, присъстващи в изходните бази данни при t₁, но липсващи t₀. Организаторите на състезанието CAFA определят два вида показатели:

• Без знание (NK): Този бенчмарк е съставен от протеини, които нямат експериментална анотация в нито едно от „пространствата от имена на GO“ и са натрупали експериментална анотация между t₀ и t₁.
• Ограничени познания (LK): Този бенчмарк е направен от протеини, които съдържат експериментални анотации в едно или две пространства от имена на GO и придобити нови анотации между t₀ и t₁ в пространство от имена без предишни анотации. Например, ако протеин с експериментални анотации в пространството от имена MF при t₀ придобие експериментални анотации в пространството от имена CC между t₀ и t1, тогава той е включен в LK бенчмарка.

Предвид тези класификации има общо шест показателя: два (т.е. NK и LK) за всяко от пространствата от имена на онтологията GO: BP, CC и MF. На практика файловете за сравнение са файлове със стойности, разделени с раздели (TSV) с две колони: ID на протеина и новия GO термин, с който е анотиран. Следното е фрагмент от файла за сравнение на BP-NK, използван за този анализ:

A0A098D6U0  GO:0106150
A0A098D8A0  GO:0106150
A0A163UT06  GO:0051568
A0A1R3RGK0  GO:1900818
A0A1U8QH20  GO:0062032
A0A286LEZ8  GO:0140380
A0A5B9  GO:0002250
A0A5B9  GO:0046631
A0A5B9  GO:0050852
A0JLT2  GO:0032968
...

Конкуренцията на CAFA като проблем за прогнозиране на връзката на графиката

Един от начините да разберете състезанието CAFA е като проблем с прогнозиране на връзката на графиката. В този контекст можем да мислим за двустранна графика, където има два типа възли: (1) генни продукти и (2) GO термини. След това анотациите са ръбовете на графиката. Въз основа на структурата на графиката и свойствата на протеините и GO термините, проблемът с конкуренцията CAFA може да бъде намален до предвиждане на липсващи връзки/анотации в графиката:

Методи

Този раздел ще премине през стъпките за генериране на данни за обучение, файлове за сравнение и предиктори. Необходимите инструменти са Anaconda, python, git, DVC, BLAST+ пакет и MATLAB.

Първоначалната стъпка за започване на изграждането на файловете с анотации за обучение и бенчмарка е да изберете времевите точки t₀ и t₁, тъй като те определят версиите на версиите на базите данни, които ще се използват. По времето, когато беше извършен този анализ, последните издания на базите данни бяха от юли 2022 г. Времевите точки за бенчмарка бяха избрани така, че продължителността на периода на „растеж на анотацията“ да е равна на тази, използвана в CAFA2 предизвикателство (което беше приблизително 10 месеца; януари до септември 2014 г.). Така за този анализ:

• t₀: ноември 2021 г
• t₁: юли 2022 г

Изтегляне на изходния код и настройка на средата

Изтегляме изходния код и настройваме нашата среда за изпълнение на необходимите скриптове:

1. Клонирайте CAFA2 хранилището и подмодулите:

git clone --recurse-submodules [email protected]:balvisio/CAFA2.git
cd CAFA2
dvc pull
pushd protein-function-prediction/protein-function-prediction-data/
dvc pull
popd

2. Създайте и активирайте Conda среда:

conda create -n <env-name> python=3.10
conda activate <env-name>
pip install -r requirements.txt

Изтегляне на изходните бази данни

За да изградим набора от бенчмаркове, трябва да извлечем информация от следните 4 бази данни:

1. UniProtKB/Swiss-Prot: е база данни, която съдържа записи на протеинови последователности заедно с ръчно анотирана функционална информация, като например GO термини. Тези записи са извлечени от проекти за секвениране на генома, както и от изследователската литература².

Най-новата версия на базата данни UniProtKB и нейните предишни версии могат да бъдат изтеглени от тук. За този проект използвах следните версии:

• Версия 2021_04: Изданието от април 2021 г., тъй като това беше наличната версия, пусната преди и най-близо до t₀. Връзката за изтегляне е тук и съдържа данни/текст, FASTA и XML формати. Сумата md5 на файла .tar.gz е bd2b1b6b0a0027f674017fe5b41fadcb.
• Версия 2022_02: Версията на изданието от февруари 2022 г., която е последната версия през юли 2022 г. (t₁), може да бъде изтеглена от тук. md5 на файла .tar.gz е c3e8dd5b138f0795aeba2b701982ac2c.

След дестартиране и разархивиране на изтеглените бази данни получихме uniprot_sprot.dat файловете, които ще са необходими за генериране на набора от данни.

2. UniProtKB/TrEMBL: също съдържа записи на протеинови последователности, но за разлика от базата данни Swiss-Prot, анотациите не се преглеждат и се генерират от изчислителни автоматични методи². Базата данни във формат FASTA може да бъде изтеглена от тук.

3. База данни за генна онтология (GO): съдържа записи за термините на GO, които съставляват GO. „GO представлява набор от знания за биологичната област. Онтологии като тази обикновено се състоят от набор от класове (или термини или концепции) с „отношения“, които действат между тях.“³ Както беше споменато, GO е подразделено на три пространства от имена: молекулярна функция, клетъчен компонент и биологичен процес. Следният фрагмент се показва като примерен запис от базата данни GO:

...
[Term]
id: GO:0032259
name: methylation
namespace: biological_process
def: "The process in which a methyl group is covalently attached to a molecule." [GOC:mah]
subset: goslim_chembl
xref: Wikipedia:Methylation
is_a: GO:0008152 ! metabolic process
...

GO версията, използвана за генериране на набора от данни, беше тази, пусната на 26 октомври 2021 г. и може да бъде изтеглена от тук. Сумата md5 на файла е 6757c819642e79e1406cad3ffcb6ea3d.

4. База данни UniProt-GOA: е база данни, която свързва описаната по-горе GO база данни с генни продукти (т.е. гени и всякакви единици, кодирани от гена, като протеин или функционални РНК). Тези връзки са това, което Консорциумът по генна онтология (GOC) нарича анотации. Източниците на тези анотации са съвместни бази данни и те трябва да съдържат доказателствен код, който описва неговия произход. Има различни категории кодове на доказателства, като например експериментален, изчислителен анализ или електронен³. За по-задълбочено обяснение и документация на базите данни на GO, моля, посетете уебсайта на GOC.

GOC публикува два типа бази данни UniProt-GOA:

• Филтрирано: „не съдържа анотации за онези видове, при които различна група от Консорциум е основно отговорна за предоставянето на GO анотации и също така изключва анотации, направени с помощта на автоматизирани методи.“⁴
• Нефилтрирано: Съдържа анотации от по-голям брой видове и включва автоматично генерирани анотации.

Използвана е филтрираната версия на базата данни UniProt-GOA, тъй като наборите от показатели се отнасят до генни продукти, които са експериментално анотирани. Базите данни UniProt-GOA могат да бъдат получени от:

• Издаване от 26 октомври 2021 г.: Сумата md5 на .gz файла е: 6d2be8b48dd0f2505092a8a2cadabd6e .
• Издаване от 1 юли 2022 г.: Сумата md5 на .gz файла е: 6f778bdc0b0ea0978e7d75f5e5289064 .

Поради „бъг“ в процеса на изграждане на тези бази данни (проблемът беше докладван „тук“), заглавките на GAF липсват. По този начин, след изтегляне на базите данни, заглавката !gaf-version: 2.0 беше добавена в началото на файловете.

Филтриране и обединяване на изходните бази данни

След като базите данни бъдат изтеглени, е необходимо да се хармонизират тези различни източници на информация. Стъпките на обработка са:

1. Базата данни Uniprot-GOA съдържа протеини както от Swiss-Prot, така и от TrEMBL. За да се подобри ефективността на скриптовете, които ще генерират обучителните набори от данни и бенчмаркове, базата данни TrEMBL се филтрира така, че да се запазят само термините, присъстващи в UniProt-GOA. Скриптът protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py прави това филтриране и приема като входни данни:

• GAF файл: UniProt-GOA база данни в момент t₀, която ще се използва за генериране на данни за обучение.
• FASTA файл: за нашия конкретен случай на употреба, файлът на базата данни TrEMBL, който трябва да бъде филтриран.

Скриптът за филтриране може да бъде извикан, както следва:

python protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py \
    <path/to/uniprot-goa-db> \
    <path/to/trembl-db> \
    <path/to/output-filtered-trembl-db>

# Example
python protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py \
    uniprot-goa/2021-10-26/goa_uniprot_all_noiea.gaf \
    ~/trembl/uniprot_trembl.fasta \
    uniprotkb/trembl/trembl_seqs_found_in_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.fasta

Полученият файл се намира в uniprotkb/trembl/trembl_seqs_found_in_goa_2021_10_26.fasta .

2. Както е обяснено в CAFA2 хартия⁵, бенчмаркът е изграден от обединението на двете бази данни UniProtKB/Swiss-Prot и UniProt-GOA. За да обединим анотациите, налични в тези бази данни, ние използвахме някои леко модифицирани скриптове, публикувани в хранилището на CAFA-Toolset.

Стъпките за обединяване на анотациите от двете бази данни са:

2.1. Изпълнете скрипта Mergedb, за да обедините базите данни UniProtKB/Swiss-Prot и UniProt-GOA в момент t₀:

# Run the next two commands only if the file has not been previously uncompressed
tar -xzvf uniprotkb/swiss-prot/2021-04/uniprot_sprot-only2021_04.tar.gz -C uniprotkb/swiss-prot/2021-04/
gunzip uniprotkb/swiss-prot/2021-04/uniprot_sprot.dat.gz

# Mergedb script
python CAFA-Toolset/Mergedb \
  -I1 <path/to/uniprotkb> \
  -I2 <path/to/uniprot-goa> \
  -O <path/to/merged-db>

# Example
rm -r workspace/
python CAFA-Toolset/Mergedb \
  -I1 uniprotkb/swiss-prot/2021-04/uniprot_sprot.dat \
  -I2 uniprot-goa/2021-10-26/goa_uniprot_all_noiea.gaf \
  -O uniprotkb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.gaf

2.2. Изпълнете скрипта Mergedb, за да обедините базите данни UniProtKB/Swiss-Prot и UniProt-GOA в момент t₁. Например:

# Run the next two commands only if the file has not been previously uncompressed
tar -zxvf uniprotkb/swiss-prot/2022-02/uniprot_sprot-only2022_02.tar.gz -C uniprotkb/swiss-prot/2022-02/
gunzip uniprotkb/swiss-prot/2022-02/uniprot_sprot.dat.gz

rm -r workspace/
python CAFA-Toolset/Mergedb \
  -I1 uniprotkb/swiss-prot/2022-02/uniprot_sprot.dat \
  -I2 uniprot-goa/2022-07-01/goa_uniprot_all_noiea.gaf \
  -O uniprotkb_2022_02_goa_2022_07_01.gaf

Ако предпочитате да пропуснете предишните стъпки, двата GAF файла, създадени по-горе, се записват в директорията uniprotkb-goa-merged на хранилището CAFA2.

Генериране на файлове с анотации на обучение

Данните за обучение за базовите предиктори са файлове с анотации, които съдържат асоциации между генни продукти и GO термини. Това са TSV файлове с две колони: ID на протеиновата последователност и свързания GO термин:

<sequence ID> <GO term ID>
...

Логично, тези анотации трябва да бъдат изградени от данни, които идват от времева точка преди или равна на t₀. Това се прави, за да се предотврати „изтичането“ на данни или с други думи, да се даде на предикторите информация, която не е била налична в момента на подаване t₀.

На високо ниво създаването на данни за обучение включва множество логически стъпки:

1. Анализирайте базата данни на Gene Ontology и генерирайте графика, където всеки възел представлява GO термин.
2. Анализирайте UniProtKB (Swiss-Prot и TrEMBL) и генерирайте представяне на всеки генен продукт, открит в базите данни.
3. Извлечете анотациите от UniProtKB и ги обединете с анотациите, намерени в базата данни UniProt-GOA на t₀.
4. Прегледайте всички анотации, намерени в предишните стъпки, и създайте логическа връзка между генни продукти от стъпка 2 и термините GO от стъпка 1.
5. Разпространете всички асоциации, открити в стъпка 4, към всички предшественици на термина GO, така че ако даден генен продукт е бил анотиран с конкретен термин GO, също да има асоциация с всеки от родителите / също по-общи термини на GO.

Кодът трябваше да обработва случаи като:

• Термините на GO могат да имат алтернативни идентификатори и някои анотации са съпоставени с тези. По този начин бяха необходими допълнителни съпоставяния, за да се съпостави с основен GO термин.
• Протеините в UniProtKB могат да имат вторични идентификатори за присъединяване, така че е необходимо допълнително картографиране, за да се отнася винаги към първичния.
• UniProt-GOA съдържа записи, които принадлежат към базите данни „ComplexPortal“ и „RNACentral“. Тези записи бяха филтрирани и не бяха включени във файловете с анотации за обучение.
• UniProt-GOA съдържа анотационни записи, които се отнасят до протеини, които са или остарели, или са били преместени в други бази данни като „UniRef“ или „UniParc“. Анотациите бяха пропуснати и не бяха включени във файловете с анотации за обучение.
• Уверете се, че файловете с пояснения не съдържат остарели GO термини.

Всички анотации се считат за:

• електронен, когато кодът на доказателството е „IEA“.
• ръководство

Читателите се насърчават да прочетат изходния код в protein-function-prediction/deepred/preprocessing/dataset.py, за да разберат по-задълбочено логиката зад генерирането на файловете с анотации.

Стъпките за създаване на файлове с анотации за обучение са:

1. Създаване на файл „Набор от данни“, който съдържа графично представяне на GO, генните продукти и техните анотации:

python protein-function-prediction/preprocessing-scripts/create_dataset.py \
    <path/to/GO Ontology/file> \
    <path/to/Swiss-Prot-db-dat-file at t0> \
    <path/to/Swiss-Prot-db-dat-file at t1> \
    <path/to/filtered/TrEMBL-db-fasta file> \
    <path/to/UniProt-GOA-db-gaf-file at t0> \
    <output/dir/to/save/serialized/dataset>

# Example  
python protein-function-prediction/preprocessing-scripts/create_dataset.py \
    protein-function-prediction/protein-function-prediction-data/data/ontology/t0/go_20211026.obo \
    protein-function-prediction/protein-function-prediction-data/data/uniprotkb/t0/uniprot_sprot.dat \
    protein-function-prediction/protein-function-prediction-data/data/uniprotkb/t1/uniprot_sprot.dat \
    protein-function-prediction/protein-function-prediction-data/data/trembl/trembl_seqs_found_in_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.fasta \
    protein-function-prediction/protein-function-prediction-data/data/goa/t0/goa_uniprot_all_noiea.gaf \
    protein-function-prediction/protein-function-prediction-data/datasets/generic/u-2021-04-g-2021-10-26/

Скриптът ще генерира файл с име dataset.bin в указаната изходна директория.

2. Създайте файловете с анотации, като използвате сериализирания набор от данни, създаден в стъпка 1, извиквайки скрипта create_annotation_files.py, както следва:

python protein-function-prediction/preprocessing-scripts/create_annotation_files.py \
    -f <path/to/dataset/file> \
    <output/dir/to/save/training/annotation/files>

# Example
python protein-function-prediction/preprocessing-scripts/create_annotation_files.py \
    -f protein-function-prediction/protein-function-prediction-data/datasets/generic/u-2021-04-g-2021-10-26/dataset.bin \
    protein-function-prediction/data/annotations/annotations_u_2021_04_g_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26

Трите окончателни файла с анотации, по един за всяко пространство от имена на онтология (т.е. BP, CC и MF) могат да бъдат намерени в директорията protein-function-prediction/data/annotations/annotations_u_2021_04_g_2021_10<output_file_name>.benchmark_{LK|NK}_{bpo|cco|mfo}{bpo|cco|mfo}.tsv. Конвенцията за именуване на тези файлове е: annotations_u_<UniProtKB_release_date>_g_<UniProt-GOA_release_date>_<ontology>.tsv . В този конкретен проект:

• u_2021_04: Издаване на UniProt април 2021 г.
• g_2021_10_26: Издаване на UniProt-GOA октомври, 2021 г.

Може да се стори странно на читателите, че изграждането на данните за обучение изисква като вход базата данни Swiss-Prot на t₁ (и точно така!). Причината за това изискване е, че тъй като UniProt-GOA беше пуснат след UniProt/Swiss-Prot (т.е. съответно октомври и април 2021 г.), определени последователности в UniProt-GOA може да не присъстват в изданието на използваната база данни Swiss-Prot. По този начин Swiss-Prot при t₁ беше използван за извличане на тези последователности, присъстващи в UniProt-GOA, но липсващи в Swiss-Prot при t₀.

Изграждане на бенчмарка

Използвайки обединените бази данни UniProtKB/Swiss-Prot и UniProt-GOA във времеви точки t₀ и t₁ (вижте раздела „Филтриране и обединяване на изходните бази данни“), ние извикваме скрипта Benchmark в хранилището на CAFA-Toolset, за да изградим бенчмарка. Този скрипт приема двата GAF файла с обединените анотации като входни данни и извежда файловете с набор от показатели във формат TSV. Този скрипт предоставя няколко опции за персонализиране на параметрите, използвани за изграждане на набора от показатели. Един от най-важните са кодовете за доказателства, които се считат за експериментални. В документа CAFA2 се посочва, че кодовете на доказателства, които се считат за експериментални, са: „EXP“, „IDA“, „IPI“, „IMP“, „IGI“, „IEP“, „TAS“ и „IC“.⁵ За създаване бенчмарк, възможно най-идентичен с този, използван в състезанието CAFA2, същите кодове на доказателства и техните двойници с висока производителност (т.е. „HTP“, „HDA“, „HMP“, „HGI“ и „HEP“) бяха предадени на скрипта Benchmark, както следва:

rm -r workspace/
python CAFA-Toolset/Benchmark \
    -I1 <path/to/t0/merged/swiss-prot-and-uniprot-goa-file> \
    -I2 <path/to/t1/merged/swiss-prot-and-uniprot-goa-file> \
    --evidence EXP IDA IPI IMP IGI IEP TAS IC HTP HDA HMP HGI HEP \
    -O <output_file_name>

# Example
rm -r workspace/
python CAFA-Toolset/Benchmark \
    -I1 uniprotkb-goa-merged/uniprotkb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.gaf \
    -I2 uniprotkb-goa-merged/uniprotkb_2022_02_goa_2022_07_01.gaf \
    --evidence EXP IDA IPI IMP IGI IEP TAS IC HTP HDA HMP HGI HEP \
    -O unfiltered_benchmark_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01

Скриптът Benchmark ще създаде директория с име workspace в директорията, от която се изпълнява скриптът, и вътре в нея ще генерира шестте бенчмарк файла с име: <output_file_name>.benchmark_{LK|NK}_{bpo|cco|mfo}.

Като се има предвид, че данните за обучение съдържат само термини на GO, присъстващи само в t₀, скриптът benchmark/filter_annotation_file.py беше използван за филтриране на термините на GO, създадени между t₀ и t₁ от генерираните файлове за сравнение. С други думи, термините на GO, които са създадени в периода на „растеж на анотациите“, не са включени в бенчмарка. Скриптът benchmark/filter_annotation_file.py приема като вход сериализирания набор от данни, генериран в раздела Генериране на файлове с анотация за обучение, както и файла за сравнение от последната стъпка. Обърнете внимание, че този скрипт трябва да се изпълни за всеки от трите файла за сравнение (т.е. bco, cco и mfo), както следва:

python benchmark/filter_annotation_file.py \
    <path/to/serialized/dataset/file> \
    <path/to/unfiltered/benchmark/file> \
    <output/path/to/save/filtered/benchmark/file>

# Example
python benchmark/filter_annotation_file.py \
    protein-function-prediction/protein-function-prediction-data/datasets/generic/u-2021-04-g-2021-10-26/dataset.bin \
    benchmark/u-2021-04-g-2021-10-26-u-2022-02-g-2022-07-01/unfiltered_benchmark_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo \
    benchmark/u-2021-04-g-2021-10-26-u-2022-02-g-2022-07-01/benchmark_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo

Филтрираните бенчмарк файлове могат да бъдат намерени в директорията benchmark/u-2021-04-g-2021-10-26-u-2022-02-g-2022-07-01 на CAFA2 хранилището.

Изградете списъка с последователности на заявки

Използвайки бенчмарк файловете, създаваме следните два файла:

1. Списък на последователности на заявки: Текстов файл, съдържащ първичните достъпи на последователностите, открити във файловете за сравнение. Това е списъкът с последователности на заявки, спрямо които ще бъдат оценени предикторите.
2. Файл FASTA на последователност от заявки: Съответният FASTA файл с последователностите от заявки. Този файл ще се използва като вход към базирания на BLAST предиктор за генериране на прогнози и оценка на неговата ефективност.

Скриптът query-sequences/generate_query_sequences.sh в хранилището на CAFA2 генерира горните два файла. Той прави следното:

1. Чете протеиновите последователности от референтен файл.
2. Филтрира повтарящите се термини.
3. Ако има идентификатор на последователност, който е вторично присъединяване, той го преобразува в основния.
4. Създава съответния FASTA файл.

Може да се извика по следния начин:

./query-sequences/generate_query_sequences.sh \
    <path/to/benchmark/file> \
    <path/to/dataset/file> \
    <path/to/save/query/sequence/list/and/fasta/files>

# Example
./query-sequences/generate_query_sequences.sh \
    benchmark/u-2021-04-g-2021-10-26-u-2022-02-g-2022-07-01/benchmark_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo \
    protein-function-prediction/protein-function-prediction-data/datasets/generic/u-2021-04-g-2021-10-26/dataset.bin \
    query-sequences/query_sequences_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo

Генерираните файлове с последователности от заявки за всяка от онтологиите могат да бъдат намерени в директорията query-sequences.

Изграждане на базовите предиктори

За да се сравни ефективността на предикторите, стандартна практика е да се създаде основен модел, който в много случаи се счита за логичен, но може би прост или наивен. Например, ако се опитваме да създадем модел, който прогнозира външната температура един ден в бъдещето, базовият модел може да бъде да предскаже, че температурата на следващия ден ще бъде средната температура за текущия ден.

Състезанието CAFA2 установи два основни предиктора:

• Базиран на BLAST предиктор: Той използва резултатите от BLAST търсене, за да генерира прогнози. Този предиктор може да бъде персонализиран да използва получения процент на идентични параметри (пидент) или BLAST E-стойност (оценка). За по-подробна информация, моля, вижте файла matlab/pfp_blast.m в хранилището на CAFA2. В тази статия стойността „pident“ ще бъде използвана в предиктора BLAST.
• Наивен предиктор: Той просто предвижда GO термин (за всички последователности от заявки) с честотата, с която се появява в данните за обучение.

Следващите раздели ще опишат стъпките за изграждане на двата предиктора и тяхната оценка. Предсказателите ще използват изходните бази данни, както и данните за обучение, които са създадени в предишните части на тази статия. За напомняне за тази статия ще използвам пространството от имена на онтологията на биологичния процес в примерните команди, но стъпките за другите две онтологии са идентични.

Заредете онтологията, анотацията на онтологията и последователностите от заявки в MATLAB

След като отворите конзолата MATLAB:

1. Добавяме поддиректорията matlab на хранилището CAFA2 към пътя на MATLAB:

>> addpath('<path/to/CAFA2/repo>/matlab');

2. Изграждаме онтологията:

>> ont = pfp_ontbuild(<path/to/ontology/at/t0>);
>>
>> % Example
>> ont = pfp_ontbuild('ontology/t0/go_20211026.obo');

3. Използваме онтологията от предишната стъпка и файла с анотация за обучение, за да изградим структурата за анотация на онтологията за обучение:

>> oa = pfp_oabuild(ont{X}, <path/to/training/annotation/file>);
>>
>> % The ontology structure contains 3 cells, one for each ontology:
>> % Biological Process (BP): ont{1}
>> % Cellular Component (CC): ont{2}
>> % Molecular Function (MF): ont{3}
>> %
>> % For example, for the BP ontology annotation:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  :
>> oa = pfp_oabuild(ont{1}, 'protein-function-prediction/protein-function-prediction-data/data/annotations/annotations_u_2021_04_g_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_bpo.tsv');

4. Зареждаме текстовия файл със списъка с последователности на заявки като клетъчен масив:

>> qseqid = readcell(<path/to/query/sequence/list/file>, 'Filetype', 'text');
>>
>> % Example
>> qseqid = readcell('query-sequences/query_sequences_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo', 'Filetype', 'text');

Предсказател, базиран на BLAST

Следващите подраздели ще опишат стъпките, необходими за изграждане на базирания на BLAST предиктор и генериране на прогнози.

Създайте базата данни BLAST

За да изградим предиктора BLAST, първо трябва да изградим базата данни BLAST. Това може да се постигне, като първо „изтеглите пакета BLAST+ от уебсайта на NCBI“. Пакетът съдържа няколко инструмента, включително makeblastdb, който изгражда база данни, използвана от BLAST при извършване на търсене за сходство на последователности. Файлът FASTA, който използвахме за създаване на базата данни BLAST, беше конкатенация на базата данни Swiss-Prot и филтрираната база данни TrEMBL, която създадохме в предишния раздел:

# Run the following command only if the file has not been previously uncompressed
gunzip uniprotkb/swiss-prot/2021-04/uniprot_sprot.fasta.gz

python uniprotkb/merged-swiss-prot-trembl/merge_dbs_for_blast.py \
    <path/to/swiss-prot-db/fasta/file> \
    <path/to/filtered-trembl-db/fasta/file> \
    <path/to/save/the/merged/databases>

# Example
python uniprotkb/merged-swiss-prot-trembl/merge_dbs_for_blast.py \
    uniprotkb/swiss-prot/2021-04/uniprot_sprot.fasta \
    uniprotkb/trembl/trembl_seqs_found_in_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.fasta \
    uniprotkb/merged-swiss-prot-trembl/merged_swiss_prot_2021_04_filtered_trembl_from_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.fasta

Конкатенираната база данни може да бъде намерена в uniprotkb/merged-swiss-prot-trembl/merged_swiss_prot_2021_04_filtered_trembl_from_goa_2021_10_26.fasta. След това, за да създадете базата данни BLAST:

makeblastdb \
    -in <path/to/fasta/file> \
    -title <some_title> \
    -dbtype prot \
    -out <output/path/to/BLAST/database>

# Example
makeblastdb \
    -in uniprotkb/merged-swiss-prot-trembl/merged_swiss_prot_2021_04_filtered_trembl_from_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26.fasta \
    -title "Merged SwissProt 2021-04 and filtered TrEMBL GOA 2021-10-26" \
    -dbtype prot \
    -out blast/db/merged_swiss_prot_2021_04_filtered_trembl_from_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26

Базата данни BLAST може да бъде намерена в директорията blast/db; съставен е от 8 файла с разширения: pdb, pjs , ptf , phr , pot , pto , pin и psq .

Подгответе BLAST резултати

След като изградихме базата данни BLAST от обучаващите последователности и последователностите на заявки FASTA файлове, имаме всичко необходимо, за да оценим BLAST предиктора. Първата стъпка е да взривите последователностите на заявките срещу последователностите за обучение, като използвате следната команда:

blastp \
    -db <path/to/blast/BLAST/database/file> \
    -query <path/to/query/sequence/list/fasta/file> \
    -outfmt "6 qseqid sseqid evalue length pident nident" \
    -out <path/to/save/blast/search/results>

# Example
blastp \
    -db blast/db/merged_swiss_prot_2021_04_filtered_trembl_from_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26 \
    -query query-sequences/query_sequences_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo.fasta \
    -outfmt "6 qseqid sseqid evalue length pident nident" \
    -out blast/blastp-results/blastp_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK_bpo.blast

Резултатите от BLAST се намират в директорията blast/blastp-results.

След като разрушим последователностите, ги зареждаме в MATLAB:

>> B = pfp_importblastp(<path/to/blast/search/results/file>);
>>
>> % Example
>> B = pfp_importblastp('blast/blastp-results/blastp_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK_bpo.blast');

и генерираме прогнозите:

>> blast = pfp_blast(qseqid, B, oa);

Наивен предсказател

Единствените елементи, необходими за създаване на наивен предиктор, са клетъчният масив qseqid и структурата на онтологичните анотации, които бяха създадени и заредени в предишните стъпки. Прогнозите на наивния предиктор се генерират от:

>> naive = pfp_naive(qseqid, oa);

Създаване на файлове за изпращане от структурите за прогнозиране

След като имаме структурата на наивния предиктор, ние запазваме резултатите във форматиран файл за изпращане, като използваме функцията prediction_to_submission_file() , така че след това да може да бъде зареден за получаване на оценките за оценка:

>> % BLAST
>> prediction_to_submission_file(blast, '<path/to/save/naive/predictions/submission/file>');
>> % Example
>> prediction_to_submission_file(blast, 'blast/blast-predictions/blast_predictions_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK_bpo.txt');
>> %
>> %
>> % Naive
>> prediction_to_submission_file(naive, '<path/to/save/naive/predictions/submission/file>');
>> % Example
>> prediction_to_submission_file(naive, 'naive/naive_predictions_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK_bpo.txt');

Файловете за подаване за предикторите BLAST и Naive могат да бъдат намерени съответно в директориите blast/blast-predictions/и naive/.

Оценяване на BLAST и наивните предиктори

След като генерираме файловете за изпращане на CAFA, най-накрая можем да изчислим резултатите на предикторите, следвайки стъпките, описани във файла README на хранилището. Тук ще покажем как да изчислим показателя fmax резултат на наивния предиктор за онтологията на BP. За предиктора BLAST, други онтологии и показатели командите са почти идентични:

• Заредете онтологията в MATLAB:

>> ont = pfp_ontbuild(<path/to/ontology/at/t0>);
>>
>> % Example
>> ont = pfp_ontbuild('ontology/go_20211026.obo');

• Изградете анотациите за сравнителен анализ:

>> oa = pfp_oabuild(ont{1}, '<path/to/benchmark/annotation/file>);
>>
>> % Example
>> oa = pfp_oabuild(ont{1}, 'benchmark/u-2021-04-g-2021-10-26-u-2022-02-g-2022-07-01/benchmark_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo');

• Заредете прогнозите:

>> pred = cafa_import(<path/to/prediction/file>, ont{X}, false);
>>
>> % Example
>> pred = cafa_import('naive/naive_predictions_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK_bpo.txt', ont{1}, false);

• Заредете идентификаторите на бенчмарк протеин:

>> benchmark = pfp_loaditem(<path/to/query/sequences/list/text/file>, 'char');
>>
>> % Example
>> benchmark = pfp_loaditem('query-sequences/query_sequences_ukb_2021_04_goa_2021protein-function-prediction/preprocessing-scripts/filter_fasta_from_goa.py26_ukb_2022_02_goa_2022_07_01_NK.bpo', 'char');

• Оценете модела:

>> fmax = pfp_seqmetric(benchmark, pred, oa, 'fmax');

Стъпките в този раздел са кодирани в скрипта matlab/evaluate_baseline_predictors.m.

Резултати

Кодът за оценка CAFA2 предоставя различни показатели за измерване на ефективността на предикторите, като максимална F-мярка (fmax), осреднена крива на прецизност-припомняне (pr) и минимално семантично разстояние (smin). За по-подробна информация, моля, вижте документацията в matlab/pfp_seqmetric.m .

За тази статия ние оценихме и двата базови предиктора, като използвахме резултата fmax на бенчмарка за липса на знания (NK). Следващите графики показват ефективността на предикторите за трите различни онтологии. Доверителните интервали от 95% бяха изчислени нормално, приблизително с помощта на метода за първоначално зареждане на набора от бенчмаркове със 100 итерации (т.е. B=100):

Ефективността на базовите предиктори варира значително в различните онтологии. Както е обяснено в документа CAFA2, тези разлики могат да бъдат резултат от (1) топологична разлика между GO дърветата, като техния размер и дълбочина, (2) факта, че термините на BP онтологията могат да се считат за по-абстрактни от, за например термините CC, които описват физическото местоположение на генен продукт, и (3) броя, както и отклоненията на анотациите в различните онтологии.

Следващите графики сравняват ефективността на предикторите Naive и BLAST, които са създадени с помощта на стъпките, описани в раздела Методи в тази статия и тези, публикувани в документа CAFA2:

Ефективността на наивните предиктори в оригиналния набор от показатели CAFA2 и този, изграден в тази статия, са много сходни. Максималната fmax разлика беше приблизително 0,05. Докато предикторът Naive 2021 превъзхожда предиктора CAFA2 в онтологията CC, той се представя по-слабо в тези за BP и MF. Това показва, че допълнителните анотации не повишават ефективността на наивния предиктор. Резултатът се очаква, като се има предвид опростената логика на предиктора, базирана само на честотата на термините GO, за да се предвидят анотации.

Предсказателят BLAST, изграден в тази статия (т.е. BLAST 2021), превъзхожда предиктора BLAST на публикацията CAFA2 съответно с 0,18 и 0,10 в онтологиите CC и BP. Разликата в производителността може да се отдаде както на: (1) допълнителни анотации, добавени между 2014 г. и 2021 г. към изходните бази данни, (2) разлики в процеса на изграждане и параметрите, използвани при генерирането на базата данни BLAST, така и на набора от бенчмаркове. Въпреки че е необходим допълнителен анализ, за ​​да се потвърди причината за тази наблюдавана разлика, този набор от показатели може да се използва за сравняване на ефективността на по-сложни предиктори.

Заключение

В тази статия преминахме през стъпките за създаване на набор от бенчмаркове и оценка на ефективността на предсказателите на протеинови анотации, използвайки насоките на състезанието CAFA2. Резултатите от ефективността на базовите предсказатели, използващи по-новите версии на изходните бази данни, са много подобни на тези, публикувани в документа CAFA2. Тези резултати показват, че допълнителните анотации, включени през годините (от 2014 до 2021) в базите данни, не осигуряват тласък в ефективността на тези прости предиктори.

Тази статия предоставя много ценна рамка за последователно сравняване на ефективността на различни методи за прогнозиране. Описаните стъпки могат да бъдат повторени, като се използват бъдещи версии на изходните бази данни, за да се оцени промяната в ефективността на предиктора, тъй като генните продукти са анотирани с повече и по-нови GO термини.

В следващите статии ще разработя предиктори, базирани на задълбочено обучение, и ще оценя тяхното представяне с помощта на бенчмарка, създаден в тази статия. Препоръчвам ви да щракнете върху бутона „Следване“, за да бъдете уведомени, когато публикувам Част II и III от тази серия. Също така имам планове да публикувам още страхотни статии, свързани с Data Science и ML, в близко бъдеще. Останете на линия!

Благодарности

Бих искал да благодаря на създателите на следните две хранилища, на които се основава голяма част от работата в тази статия:

• CAFA2: https://github.com/yuxjiang/CAFA2
• CAFA-Toolset: https://github.com/arkatebi/CAFA-Toolset

Препратки

[1] https://biofunctionprediction.org/cafa/

[2] https://en.wikipedia.org/wiki/UniProt

[3] http://geneontology.org/docs/ontology-documentation/

[4] http://geneontology.org/docs/faq/

[5] Jiang, Y., Oron, T., Clark, W. et al. Разширената оценка на методите за прогнозиране на протеиновата функция показва подобрение в точността. Genome Biol 17, 184 (2016). https://doi.org/10.1186/s13059-016-1037-6

Забележка относно изображенията в статията: Освен ако не е отбелязано друго, всички изображения са създадени от автора на тази статия.