В этом примере показано, как обогатить данные об экспрессии гена микромассивов с помощью Генных отношений Онтологии.
Генная Онтология является управляемым методом для описания терминов, связанных с генами в любом организме. Когда больше генных данных получено из организмов, они аннотируются с помощью Генной Онтологии. Генная Онтология сделана из трех меньших онтологий или аспектов: Молекулярная Функция, Биологический Процесс и Сотовый Компонент. Каждая из этих онтологий содержит термины, которые организованы в направленном графе без петель с этими тремя терминами как корни. Корни являются самыми широкими терминами, относящимися к генам. Условия еще дальше от корней становятся более конкретными. Для этого примера вы будете использовать микроданные массива из Аналитического примера Профиля Экспрессии гена, чтобы посмотреть на значение интересных генов и Генных терминов Онтологии, которые сопоставлены с зондами микромассивов. А именно, вы далее займетесь расследованиями, чтобы определить, если набор генов, что кластер вместе также вовлечен в общую молекулярную функцию.
База данных Gene Ontology загружается в объект MATLAB® использование Bioinformatics Toolbox geneont функция.
GO = geneont('live',true); % this step takes a while get(GO)
default_namespace: 'gene_ontology'
format_version: '1.2'
data_version: 'releases/2020-07-16'
version: ''
date: ''
saved_by: ''
auto_generated_by: ''
subsetdef: {17×1 cell}
import: ''
synonymtypedef: 'systematic_synonym "Systematic synonym" EXACT'
idspace: ''
default_relationship_id_prefix: ''
id_mapping: ''
remark: 'Includes Ontology(OntologyID(OntologyIRI(<http://purl.obolibrary.org/obo/go/never_in_taxon.owl>))) [Axioms: 18 Logical Axioms: 0]'
typeref: ''
unrecognized_tag: {'ontology' 'go'}
Terms: [47203×1 geneont.term]
Каждый Генный термин Онтологии имеет инвентарный номер, который является семизначным числом, которому предшествует префикс, 'ПОЙДИТЕ':. чтобы посмотреть на конкретный термин, вы сначала создаете подонтологию путем индексирования объекта и затем смотрите 'термины' свойство. Хеш-таблица реализована в объекте GO для того, чтобы эффективно искать идентификаторы термина.
GO(5840).terms % the ribosome Gene Ontology term
id: 5840
name: 'ribosome'
ontology: 'cellular component'
definition: '"An intracellular organelle, about 200 A in diameter, consisting of RNA and protein. It is the site of protein biosynthesis resulting from translation of messenger RNA (mRNA). It consists of two subunits, one large and one small, each containing only protein and RNA. Both the ribosome and its subunits are characterized by their sedimentation coefficients, expressed in Svedberg units (symbol: S). Hence, the prokaryotic ribosome (70S) comprises a large (50S) subunit and a small (30S) subunit, while the eukaryotic ribosome (80S) comprises a large (60S) subunit and a small (40S) subunit. Two sites on the ribosomal large subunit are involved in translation, namely the aminoacyl site (A site) and peptidyl site (P site). Ribosomes from prokaryotes, eukaryotes, mitochondria, and chloroplasts have characteristically distinct ribosomal proteins." [ISBN:0198506732]'
comment: ''
synonym: {4×2 cell}
is_a: 43232
part_of: [0×1 double]
obsolete: 0
Этот термин представляет 'рибосому', и это имеет поля для is_a и part_of. Эти поля представляют отношения между Генными терминами Онтологии. Генные термины Онтологии могут рассматриваться как узлы в графе без петель. Можно пересечь такие отношения с методами getancestors, getdescendants, getrelatives, и getmatrix. Например, getancestors метод возвращает любой менее конкретный термин, чем 'рибосома' (i.e., его родительские элементы в графике).
ancestors = getancestors(GO,5840) riboanc = GO(ancestors)
ancestors =
5575
5622
5840
43226
43228
43229
43232
110165
Gene Ontology object with 8 Terms.
Чтобы визуализировать эти отношения, используйте biograph функционируйте и getmatrix метод от Bioinformatics Toolbox. getmatrix метод возвращает квадратную матрицу отношений данного объекта Gene Ontology. Этот график иногда называется 'вызванным' графиком.
cm = getmatrix(riboanc);
BG = biograph(cm,get(riboanc.Terms,'name'))
view(BG)
Biograph object with 8 nodes and 9 edges.
Чтобы показать, как Генная информация об Онтологии полезна, вы посмотрите на микроданные массива из Аналитического примера Профиля Экспрессии гена. Эти данные имеют 6 400 генов на микромассиве, которые связаны со многими различными аспектами экспрессии гена дрожжей. Небольшая часть их может показать интересное поведение для этого эксперимента микромассивов. Вы будете использовать Генную Онтологию, чтобы лучше изучить, если и как эти гены связаны в ячейке. Полные данные о дрожжах могут быть найдены в Веб-сайте NCBI.
load yeastdata whos yeastvalues genes
Name Size Bytes Class Attributes genes 6400x1 755322 cell yeastvalues 6400x7 358400 double
Пример "Анализ Профиля Экспрессии гена" показывает несколько способов кластеризировать данные из эксперимента. В этом примере кластеризация K-средних значений используется, чтобы выбрать группу приблизительно из 240 генов для исследования.
Во-первых, данным нужна чистка. На чипе существуют некоторые пустые места, и некоторые гены имеют отсутствующие значения. В этом примере пустые пятна удалены из набора данных и knnimpute функция приписывает отсутствующие значения (отмеченный NaNs).
% Remove data for empty spots emptySpots = strcmp('EMPTY',genes); yeastvalues = yeastvalues(~emptySpots,:); genes = genes(~emptySpots); fprintf('Number of genes after removing empty spots is %d.\n',numel(genes)) % Impute missing values yeastvalues = knnimpute(yeastvalues);
Number of genes after removing empty spots is 6314.
Затем функциональный genelowvalfilter удаляет гены с низким полным выражением. В этом примере довольно большое значение используется, чтобы отфильтровать гены. Аналитический пример Профиля Экспрессии гена показывает альтернативные методы фильтрации.
mask = genelowvalfilter(yeastvalues,genes,'absval',log2(3.5)); highexpIdx = find(mask); yeastvalueshighexp = yeastvalues(highexpIdx,:); fprintf('Number of genes with high expression is %d.\n',numel(highexpIdx))
Number of genes with high expression is 613.
kmeans функция от Statistics and Machine Learning Toolbox™ группирует данные в четыре кластера с помощью корреляции в качестве метрики расстояния.
rng default [cidx, ctrs] = kmeans(yeastvalueshighexp, 4, 'dist','corr', 'rep',20); for c = 1:4 subplot(2,2,c); plot(times,yeastvalueshighexp((cidx == c),:)'); axis tight end sgtitle('K-Means Clustering of Profiles');
Графики показывают четыре довольно различных кластера. Путем рассмотрения центроидов кластеров вы видите ясно, как они отличаются.
figure for c = 1:4 subplot(2,2,c); plot(times,ctrs(c,:)','linewidth',4,'color','k'); axis tight axis off % turn off the axis end sgtitle('Centroids of K-Means Clustering of Profiles');
Первый кластер в левом верхнем угле представляет гены, которые отрегулированы с их уровнями экспрессии, уменьшающимися немного в итоговом чипе. Гены в этом кластере будут подмножеством, используемым для остатка от этого эксперимента.
clusterIdx = highexpIdx(cidx==1);
fprintf('Number of genes in the first cluster is %d.\n',numel(clusterIdx))
Number of genes in the first cluster is 140.
Много Проектов Генома взаимодействуют с Генным Консорциумом Онтологии при аннотировании генов. Генные аннотации для нескольких организмов могут быть найдены в Генном Веб-сайте Онтологии. Кроме того, аннотации для отдельных организмов могут быть найдены в их соответствующих веб-сайтах (таких как база данных Yeast Genome). Эти аннотации часто обновляются и обычно курируются членами группы генома для каждого организма. NCBI также имеет коллективный список генных аннотаций, которые относятся к их базе данных Entrez Gene. Эти файлы аннотации состоят из больших списков генов и их связанных Генных терминов Онтологии. Эти файлы следуют за структурой, заданной Генным Консорциумом Онтологии. Функциональный goannotread проанализирует эти несжатые файлы и поместит информацию в структуру MATLAB. Файл yeastgenes.sgd был получен из Генного сайта Аннотации Онтологии.
Для этого анализа вы посмотрите на гены, которые аннотируются как молекулярная функция (т.е. поле 'Aspect' установлено в 'F'). Однако вы могли расширить этот анализ, чтобы видеть, вовлечены ли кластеризованные гены в общие процессы ('P') или совмещены в том же отсеке ячейки ('C'). Поля, которые представляют интерес, являются названием гена и сопоставленным ID. В ИДУТ файлы Аннотации, они имеют имена полей DB_Object_Symbol и GOid, соответственно. Создайте карту для эффективного поиска Названия гена. Заметьте, что не каждый ген от этих 6 314 генов на микромассиве аннотируется. Обратите внимание на то, что данные в общедоступных репозиториях часто курируются и обновляются; поэтому результаты этого примера могут немного отличаться, когда вы используете актуальные наборы данных. Также возможно, что вы получаете предупреждения о недопустимых или устаревших идентификаторах из-за устаревшего yeastgenes.sgd генный файл аннотации.
SGDann = goannotread('yeastgenes.sgd','Aspect','F',... 'Fields',{'DB_Object_Symbol','GOid'}); SGDmap = containers.Map(); for i = 1:numel(SGDann) key = SGDann(i).DB_Object_Symbol; if isKey(SGDmap,key) SGDmap(key) = [SGDmap(key) SGDann(i).GOid]; else SGDmap(key) = SGDann(i).GOid; end end fprintf('Number of annotated genes related to molecular function is %d.\n',SGDmap.Count) fprintf('Number of unique GO terms associated to annotated genes is %d.\n',numel(unique([SGDann.GOid]))) fprintf('Number of gene-GO term associations is %d.\n',numel(SGDann))
Number of annotated genes related to molecular function is 6428. Number of unique GO terms associated to annotated genes is 2017. Number of gene-GO term associations is 31352.
Для каждого термина в Генной Онтологии вы считаете следующие два элемента:
Количество генов, которые аннотируются к термину.
Количество под - или сверхописанные гены, которые аннотируются к термину.
На основе этой информации можно статистически определить, как часто Генные термины Онтологии сверх - или underrepresentated в эксперименте микромассивов.
Существуют некоторые случаи, где генные наборы точно не аннотируются. Поэтому вы создаете эти счета путем распространения к соседним терминам. По определению более определенные Генные термины Онтологии включены в термины предка. Таким образом, если ген аннотируется к конкретному термину, можно также хотеть увеличить счет для некоторых или всех его терминов предка. Используйте getrelatives, getdescendants, или getancestors протестировать различные схемы распространения. В этом примере вы используете getrelatives получить предков и потомков каждого термина до 1 уровня, который является значением по умолчанию, в основном чтобы преодолеть неточно аннотируемый генный набор.
m = GO.Terms(end).id; % gets the last term id geneschipcount = zeros(m,1); % a vector of GO term counts for the entire chip. genesclustercount = zeros(m,1); % a vector of GO term counts for interesting genes. for i = 1:numel(genes) if isKey(SGDmap,genes{i}) goid = getrelatives(GO,SGDmap(genes{i})); % update vector counts geneschipcount(goid) = geneschipcount(goid) + 1; if (any(i == clusterIdx)) genesclustercount(goid) = genesclustercount(goid) + 1; end end end
Warning: Invalid or obsolete IDs: 15238 15238 15238 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 36459 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 44212 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 1077 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 1077 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 982 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 4004 Check that the annotation file is up to date. Warning: Invalid or obsolete IDs: 982 Check that the annotation file is up to date.
Можно найти старшие значащие аннотируемые термины путем рассмотрения вероятностей, что термины считаются случайно. Используйте гипергеометрическую функцию распределения вероятностей (hygepdf), который вычисляет статистическое значение того, что вытянул определенное количество успехов из общего количества ничьих от населения. Расчетное p-значение является вероятностью получения такой тестовой статистики, которая является количествами счета интересных генов в этом примере. Эта функция возвращает p-значение, сопоставленное в каждый термин, и можно создать список старшего значащего, ИДУТ термины путем упорядоченного расположения p-значений.
pvalues = hygepdf(genesclustercount,max(geneschipcount),... max(genesclustercount),geneschipcount); [dummy,idx] = sort(pvalues); % create a report report = sprintf('GO Term p-val counts definition\n'); for i = 1:10 term = idx(i); report = sprintf('%s%s\t%-1.4f\t%-d / %-d\t%s...\n', report, ... char(num2goid(term)), pvalues(term),... genesclustercount(term),geneschipcount(term),... GO(term).Term.definition(2:min(end,60))); end disp(report);
GO Term p-val counts definition GO:0000104 0.0141 1 / 1 Catalysis of the reaction: succinate + acceptor = fumarate ... GO:0003995 0.0141 1 / 1 Catalysis of the reaction: acyl-CoA + acceptor = 2,3-dehydr... GO:0008177 0.0141 1 / 1 Catalysis of the reaction: succinate + ubiquinone = fumarat... GO:0009046 0.0141 1 / 1 Catalysis of the cleavage of the D-alanyl-D-alanine bond in... GO:0009917 0.0141 1 / 1 Catalysis of the removal of a C-5 double bond in the B ring... GO:0009918 0.0141 1 / 1 Catalysis of the reaction: 5-dehydroepisterol = 24-methylen... GO:0016166 0.0141 1 / 1 Catalysis of the dehydrogenation of phytoene to produce a c... GO:0016628 0.0141 1 / 1 Catalysis of an oxidation-reduction (redox) reaction in whi... GO:0016632 0.0141 1 / 1 Catalysis of an oxidation-reduction (redox) reaction in whi... GO:0016634 0.0141 1 / 1 Catalysis of an oxidation-reduction (redox) reaction in whi...
Можно использовать методы, описанные ранее в этом примере, чтобы узнать больше о терминах, которые кажутся высокими в этом списке.
Сначала посмотрите на предков главного элемента в списке.
topItem = idx(1);
GO(topItem).terms % the most significant gene
topItemAncestors = getancestors(GO,topItem)
id: 104
name: 'succinate dehydrogenase activity'
ontology: 'molecular function'
definition: '"Catalysis of the reaction: succinate + acceptor = fumarate + reduced acceptor." [GOC:kd]'
comment: ''
synonym: {11×2 cell}
is_a: 16627
part_of: [0×1 double]
obsolete: 0
topItemAncestors =
104
3674
3824
16491
16627
Посмотрите на список для второго элемента, чтобы видеть некоторых из тех же предков.
secondItem = idx(2);
GO(secondItem).terms % the second most significant gene
secondItemAncestors = getancestors(GO,secondItem)
id: 3995
name: 'acyl-CoA dehydrogenase activity'
ontology: 'molecular function'
definition: '"Catalysis of the reaction: acyl-CoA + acceptor = 2,3-dehydroacyl-CoA + reduced acceptor." [EC:1.3.99.3]'
comment: ''
synonym: {14×2 cell}
is_a: 16627
part_of: [0×1 double]
obsolete: 0
secondItemAncestors =
3674
3824
3995
16491
16627
Можно теперь создать подонтологию, которая включает предков десяти (например), старших значащих терминов, и визуализируйте это использование biograph функция.
subGO = GO(getancestors(GO,idx(1:10)))
[cm,acc,rels] = getmatrix(subGO);
BG = biograph(cm,get(subGO.Terms,'name'))
Gene Ontology object with 23 Terms. Biograph object with 23 nodes and 26 edges.
Используйте p-значения, вычисленные прежде, чтобы присвоить цвет вершинам графика. В этом примере используется произвольная карта цветов, где яркий красный является старшим значащим и ярким зеленым, младше значащее.
for i = 1:numel(acc) pval = pvalues(acc(i)); color = [(1-pval).^(10),pval.^(1/10),0.3]; BG.Nodes(i).Color = color; BG.Nodes(i).Label = num2str(acc(i)); % add info to datatips end view(BG);
[1] http://www.geneontology.org /
[2] http://www.yeastgenome.org /
[3] Джентльмен, R. 'Основной ИДУТ Использование'. Биопроводниковая виньетка 16 мая 2005 http://bioconductor.org/docs/vignettes.html
[4] Джентльмен, R. 'Используя ИДУТ для Статистических анализов'. Биопроводниковая виньетка 16 мая 2005 http://bioconductor.org/docs/vignettes.html