Этот пример показывает, как обогатить данные экспрессии генов микрочипов с использованием отношений генной онтологии.
Генная онтология - контролируемый метод описания терминов, связанных с генами в любом организме. Поскольку из организмов получают больше данных о генах, они аннотируются с использованием Gene Ontology. Генная онтология состоит из трех меньших онтологий или аспектов: молекулярная функция, биологический процесс и клеточный компонент. Каждая из этих онтологий содержит члены, которые организованы в направленный ациклический граф с этими тремя членами в качестве корней. Корни являются широчайшими терминами, относящимися к генам. Термины дальше от корней становятся более конкретными. Для этого примера вы будете использовать данные микрочипов из примера анализа профиля экспрессии генов, чтобы посмотреть на значимость интересных генов и терминов генной онтологии, которые связаны с зондами микрочипов. Более конкретно, вы дополнительно изучите, чтобы определить, участвует ли набор генов, которые объединяются вместе, в общей молекулярной функции.
База данных Gene Ontology загружается в объект MATLAB ® с помощью панели инструментов биоинформатики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:'. Чтобы просмотреть конкретный термин, необходимо сначала создать субонтологию путем подскрипта объекта, а затем проверить свойство «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. Эти поля представляют взаимосвязи между терминами генной онтологии. Термины генной онтологии можно рассматривать как узлы в ациклическом графе. Вы можете пройти такие отношения с методами getancestors, getdescendants, getrelatives, и getmatrix. Например, getancestors метод возвращает любой менее конкретный член, чем «рибосома» (т.е. его родители в графе).
ancestors = getancestors(GO,5840) riboanc = GO(ancestors)
ancestors =
5575
5622
5840
43226
43228
43229
43232
110165
Gene Ontology object with 8 Terms.
Чтобы визуализировать эти отношения, используйте biograph функции и getmatrix метод из панели инструментов биоинформатики. getmatrix метод возвращает квадратную матрицу отношений данного объекта Gene Ontology. Этот граф иногда называют «индуцированным» графом.
cm = getmatrix(riboanc);
BG = biograph(cm,get(riboanc.Terms,'name'))
view(BG)
Biograph object with 8 nodes and 9 edges.
Чтобы показать, как полезна информация о генной онтологии, вы посмотрите на данные микрочипов из примера анализа профиля экспрессии генов. Эти данные имеют 6400 генов на микрочипе, которые участвуют во многих различных аспектах экспрессии генов дрожжей. Небольшая часть из них может показать интересное поведение для этого эксперимента с микрочипами. Вы будете использовать генную онтологию, чтобы лучше понять, связаны ли эти гены в клетке. Полные данные о дрожжах можно найти на веб-сайте 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 функция из 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. Кроме того, аннотации для отдельных организмов можно найти на их соответствующих веб-сайтах (таких как база данных генома дрожжей). Эти аннотации часто обновляются и обычно обрабатываются членами группы генома для каждого организма. NCBI также имеет сводный список аннотаций генов, которые относятся к их базе данных Entrez Gene. Эти файлы аннотаций состоят из больших списков генов и связанных с ними терминов Gene Ontology. Эти файлы соответствуют структуре, определенной консорциумом Gene Ontology. Функция goannotread разберет эти несжатые файлы и поместит информацию в структуру MATLAB. Файл yeastgenes.sgd был получен с сайта аннотации генной онтологии.
Для этого анализа вы посмотрите на гены, которые аннотированы как молекулярная функция (то есть поле «Аспект» установлено на «F»). Тем не менее, вы можете расширить этот анализ, чтобы увидеть, участвуют ли кластерные гены в общих процессах («P») или находятся в одном клеточном компартменте («C»). Представляющими интерес полями являются символ гена и связанный идентификатор. В файлах аннотаций GO они имеют имена полей DB_Object_Symbol и GOid соответственно. Создайте карту для эффективного поиска символа гена. Обратите внимание, что не каждый ген из генов 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. Таким образом, если ген аннотирован к определенному термину, можно также увеличить количество для некоторых или всех его предков. Использовать 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, связанное с каждым термином, и можно создать список наиболее значимых терминов 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] Джентльмен, Р. 'Базовое использование ГО'. Биокондукторная виньетка 16 мая 2005 года http://bioconductor.org/docs/vignettes.html
[4] Джентльмен, Р. «Использование ГО для статистического анализа». Биокондукторная виньетка 16 мая 2005 года http://bioconductor.org/docs/vignettes.html