Этот пример показывает, как обогатить данные об экспрессии гена микромассивов с помощью Генных отношений Онтологии.
Генная Онтология является управляемым методом для описания условий, связанных с генами в любом организме. Когда больше генных данных получено из организмов, они аннотируются с помощью Генной Онтологии. Генная Онтология сделана из трех меньших онтологий или аспектов: Молекулярная Функция, Биологический Процесс и Сотовый Компонент. Каждая из этих онтологий содержит условия, которые организованы в направленном графе без петель с этими тремя условиями с должности корней. Корни являются самыми широкими условиями, относящимися к генам. Условия еще дальше от корней становятся более конкретными. Для этого примера вы будете использовать микроданные массива от Аналитического примера Профиля Экспрессии гена, чтобы посмотреть на значение интересных генов и Генных условий Онтологии, которые сопоставлены с зондами микромассивов. А именно, вы далее займетесь расследованиями, чтобы определить, если набор генов, что кластер вместе также вовлечен в общую молекулярную функцию.
База данных Gene Ontology загружается в объект MATLAB® использование Bioinformatics Toolbox функция geneont.
GO = geneont('live',true); % this step takes a while get(GO)
default_namespace: 'gene_ontology'
format_version: '1.0'
data_version: ''
version: ''
date: '23:03:2018 01:38'
saved_by: 'jenkins-slave'
auto_generated_by: ''
subsetdef: {21×1 cell}
import: ''
synonymtypedef: ''
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: {3×2 cell}
Terms: [47113×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: [3×1 double]
part_of: [0×1 double]
obsolete: 0
Этот термин представляет 'рибосому', и это имеет поля для is_a и part_of. Эти поля представляют отношения между Генными условиями Онтологии. Генные условия Онтологии могут рассматриваться как узлы в графе без петель. Можно пересечь такие отношения с методами getancestors, getdescendants, getrelatives и getmatrix. Например, метод getancestors возвращает любой менее конкретный термин, чем 'рибосома' (т.е. ее родительские элементы в графике).
ancestors = getancestors(GO,5840) riboanc = GO(ancestors)
ancestors =
5575
5622
5623
5737
5840
32991
43226
43228
43229
43232
44424
44444
44464
1990904
Gene Ontology object with 14 Terms.
Чтобы визуализировать эти отношения, используйте функцию biograph и метод getmatrix от Bioinformatics Toolbox. Метод getmatrix возвращает квадратную матрицу отношений данного объекта Gene Ontology. Этот график иногда называется 'вызванным' графиком.
cm = getmatrix(riboanc);
BG = biograph(cm,get(riboanc.Terms,'name'))
view(BG)
Biograph object with 14 nodes and 20 edges.
Чтобы показать, как Генная информация об Онтологии полезна, вы посмотрите на микроданные массива от Аналитического примера Профиля Экспрессии гена. Эти данные имеют 6 400 генов на микромассиве, которые связаны со многими различными аспектами экспрессии гена дрожжей. Небольшая часть их может показать интересное поведение для этого эксперимента микромассивов. Вы будете использовать Генную Онтологию, чтобы лучше понять, если и как эти гены связаны в ячейке. Полные данные о дрожжах могут быть найдены в Веб-сайте NCBI.
load yeastdata whos yeastvalues genes
Name Size Bytes Class Attributes genes 6400x1 806522 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 suptitle('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 suptitle('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
Можно найти старшие значащие аннотируемые условия путем рассмотрения вероятностей, что условия считаются случайно. Используйте гипергеометрическую функцию распределения вероятностей (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." [EC:1.3.99.1, 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 24 Terms. Biograph object with 24 nodes and 27 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