Этот пример показывает, как обогатить данные экспрессии генов микромассивов, используя отношения Gene Ontology.
Генная онтология является контролируемым методом для описания терминов, связанных с генами в любом организме. Поскольку больше данных генов получают от организмов, это аннотируется с помощью Gene Ontology. Генная онтология состоит из трех небольших онтологий или аспектов: молекулярной функции, биологического процесса и клеточного компонента. Каждая из этих онтологий содержит термины, которые организованы в ориентированном ациклическом графике с этими тремя терминами в качестве корней. Корни являются самыми широкими терминами, относящимися к генам. Термины дальше от корней становятся более конкретными. В этом примере вы будете использовать данные микромассивов из примера Gene Expression Profile Analysis, чтобы посмотреть на значимость интересных генов и терминов Gene Ontology, которые связаны с зондами микромассивов. Более конкретно, вы дополнительно исследуете, чтобы определить, участвует ли набор генов, которые объединяются вместе, в общей молекулярной функции.
База данных 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]
Каждый термин Gene Ontology имеет номер присоединения, который является семизначным числом, предшествующим префиксу 'GO:'. Чтобы взглянуть на определенный термин, сначала создайте субонтологию путем индексирования объекта, а затем проверьте свойство 'terms'. Хеш-таблица реализована в объекте 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. Эти поля представляют отношения между терминами Gene Ontology. Термины Gene Ontology можно рассматривать как узлы в ациклическом графике. Можно обойти такие отношения с помощью методов getancestors, getdescendants, getrelatives, и getmatrix. Для примера, getancestors метод возвращает любой менее специфический термин, чем 'ribosome' (т.е. его родительские элементы в графике).
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. The getmatrix метод возвращает квадратную матрицу отношений данного объекта Gene Ontology. Этот график иногда называют 'индуцированным' графиком.
cm = getmatrix(riboanc);
BG = biograph(cm,get(riboanc.Terms,'name'))
view(BG)
Biograph object with 8 nodes and 9 edges.
Чтобы показать, как информация о генной онтологии полезна, вы посмотрите на данные микромассивов из примера анализа профиля экспрессии генов. Эти данные содержат 6400 генов на микромассивы, которые участвуют во многих различных аспектах экспрессии генов дрожжей. Небольшой фрагмент из них может показать интересное поведение для этого микромассива эксперимента. Вы будете использовать Gene Ontology, чтобы лучше понять, связаны ли эти гены в камере и как они связаны. Полные данные о дрожжах можно найти на веб-сайте 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.
The 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.
Многие проекты генома взаимодействуют с консорциумом Gene Ontology при аннотировании генов. Аннотации генов для нескольких организмов можно найти на сайте Gene Ontology. В сложении аннотации для отдельных организмов можно найти на их соответствующих веб-сайтах (таких как база данных генома дрожжей). Эти аннотации часто обновляются и обычно кураторятся представителями группы генома для каждого организма. NCBI также имеет сводный список аннотаций генов, которые относятся к их базе данных Entrez Gene. Эти файлы аннотаций состоят из больших списков генов и связанных с ними терминов Gene Ontology. Эти файлы соответствуют структуре, заданной консорциумом Gene Ontology. Функция goannotread проанализирует эти несжатые файлы и поместит информацию в структуру MATLAB. Файл yeastgenes.sgd был получен из сайта аннотации генной онтологии.
Для этого анализа вы рассмотрим гены, которые аннотируются как молекулярная функция (то есть, поле 'Aspect' установлено на 'F'). Однако можно расширить этот анализ, чтобы увидеть, участвуют ли кластеризованные гены в общих процессах ('P') или находятся вместе в том же клеточном отсеке ('C'). Поля, которые представляют интерес, являются символом гена и связанным с ним ID. В файлах аннотации GO они имеют имена полей DB_Object_Symbol и GOid, соответственно. Создайте карту для эффективного поиска Gene Symbol. Заметьте, что не каждый ген из генов 6314 в микромассиве аннотирован. Обратите внимание, что данные в общедоступных хранилищах часто кураторятся и обновляются; поэтому результаты этого примера могут несколько отличаться при использовании современных наборов данных. Также возможно, что вы получаете предупреждения о недопустимых или устаревших идентификаторах из-за устаревшего 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.
Для каждого термина в генной онтологии вы подсчитываете следующие два элементов:
Количество генов, которые аннотированы к термину.
Количество недостаточно или сверхэкспрессированных генов, которые аннотируются к термину.
Основываясь на этой информации, вы можете статистически определить, как часто условия Gene Ontology чрезмерно или недостаточно представлены в эксперименте с микромассивами.
Существуют некоторые случаи, когда наборы генов не точно аннотируются. Поэтому вы создаете эти талисманы путем распространения на соседние условия. По определению, более специфические термины Gene Ontology включены в предковые термины. Таким образом, если ген аннотирован к конкретному термину, вы также можете захотеть увеличить количество для некоторых или всех его предковых терминов. Использование 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-значение является вероятностью получения такой тестовой статистики, которая является подсчетом tally интересных генов в этом примере. Эта функция возвращает p-значение, сопоставленное с каждым членом, и можно создать список наиболее значимых членов GO, упорядочив 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] Джентльмен, Р. 'Basic GO Usage'. Биокондукторная виньетка 16 мая 2005 http://bioconductor.org/docs/vignettes.html
[4] Джентльмен, Р. 'Использование GO для статистического анализа'. Биокондукторная виньетка 16 мая 2005 http://bioconductor.org/docs/vignettes.html