В этом примере показано, как извлечь серию данных экспрессии генов (GSE) из NCBI Gene Expression Omnibus (GEO) и выполнить базовый анализ профилей экспрессии.
NCBI Gene Expression Omnibus (GEO) является крупнейшим публичным хранилищем высокопроизводительных экспериментальных данных микрочипов. Данные GEO имеют четыре типа объектов, включая GEO Platform (GPL), GEO Sample (GSM), GEO Series (GSE) и GUO DataSet (GDS).
Запись Платформы описывает список элементов на массиве в эксперименте (например, кДНК, олигонуклеотидные зонды). Каждой записи Платформы присваивается уникальный и стабильный регистрационный номер GEO (GPLxxx).
Запись образца описывает условия, в которых обрабатывался отдельный образец, манипуляции, которым он подвергся, и измерение плотности каждого элемента, полученного из него. Каждой записи образца присваивается уникальный и стабильный регистрационный номер GEO (GSMxxx).
Запись серии определяет группу связанных образцов и предоставляет фокус и описание всего исследования. Записи серии могут также содержать таблицы, описывающие извлеченные данные, итоговые выводы или анализы. Каждой записи серии присваивается уникальный регистрационный номер GEO (GSExxx).
Запись DataSet (GDSxxx) представляет собой обработанную коллекцию биологически и статистически сопоставимых образцов GEO. Наборы данных GEO (GDSxxx) являются управляемыми наборами данных GEO Sample.
Более подробную информацию о ГЕО можно найти в разделе Обзор ГЕО. Биоинформатика Toolbox™ предоставляет функции, которые могут извлекать и анализировать файлы данных формата GEO. Данные GSE, GSM, GSD и GPL могут быть получены с помощью getgeodata функция. getgeodata функция также может сохранять извлеченные данные в текстовом файле. Записи серии GEO доступны в файлах формата SOFT и в файлах формата текста с разделителями табуляции. Функция geoseriesread считывает текстовый файл формата серии GEO. geosoftread функция считывает обычно довольно большие файлы формата SOFT.
В этом примере выполняется извлечение набора GSE5847 данных из базы данных GEO и статистическое тестирование данных. GEO Series GSE5847 содержит экспериментальные данные исследования экспрессии генов опухолевой стромы и клеток эпителия из 15 случаев воспалительного рака молочной железы (IBC) и 35 случаев невоспалительного рака молочной железы (Boersma et al. 2008).
Функция getgeodata возвращает структуру, содержащую данные, полученные из базы данных GEO. Также можно сохранить возвращенные данные в исходном формате в локальной файловой системе для использования в последующих сессиях MATLAB ®. Следует отметить, что данные в публичных хранилищах часто обрабатываются и обновляются; поэтому результаты этого примера могут несколько отличаться при использовании актуальных наборов данных.
gseData = getgeodata('GSE5847', 'ToFile', 'GSE5847.txt')Используйте geoseriesread для разбора ранее загруженного файла формата текста GSE.
gseData = geoseriesread('GSE5847.txt')
gseData =
struct with fields:
Header: [1x1 struct]
Data: [22283x95 bioma.data.DataMatrix]
Возвращенная структура содержит Header поле, в котором хранятся метаданные данных серии, и Data поле, в котором хранятся данные матрицы серии.
Данные матрицы GSE5847 в Data сохраняются как объект DataMatrix. Объект DataMatrix представляет собой структуру данных, такую как массив MATLAB 2-D, но с дополнительными метаданными имен строк и столбцов. Доступ к свойствам DataMatrix можно получить, как и к другим объектам MATLAB.
get(gseData.Data)
Name: ''
RowNames: {22283x1 cell}
ColNames: {1x95 cell}
NRows: 22283
NCols: 95
NDims: 2
ElementClass: 'double'
Имена строк являются идентификаторами наборов зондов в массиве; имена столбцов представляют собой номера образцов GEO.
gseData.Data(1:5, 1:5)
ans =
GSM136326 GSM136327 GSM136328 GSM136329 GSM136330
1007_s_at 10.45 9.3995 9.4248 9.4729 9.2788
1053_at 5.7195 4.8493 4.7321 4.7289 5.3264
117_at 5.9387 6.0833 6.448 6.1769 6.5446
121_at 8.0231 7.8947 8.345 8.1632 8.2338
1255_g_at 3.9548 3.9632 3.9641 4.0878 3.9989
Метаданные серии хранятся в Header поле. Структура содержит информацию о серии в Header.Series и примеры информации в Header.Sample поле.
gseData.Header
ans =
struct with fields:
Series: [1x1 struct]
Samples: [1x1 struct]
Поле Series содержит название эксперимента и идентификатор платформы GEO.
gseData.Header.Series
ans =
struct with fields:
title: 'Tumor and stroma from breast by LCM'
geo_accession: 'GSE5847'
status: 'Public on Sep 30 2007'
submission_date: 'Sep 15 2006'
last_update_date: 'Nov 14 2012'
pubmed_id: '17999412'
summary: 'Tumor epithelium and surrounding stromal cells were isolated using laser capture microdissection of human breast cancer to examine differences in gene expression based on tissue types from inflammatory and non-inflammatory breast cancer...'
overall_design: 'We applied LCM to obtain samples enriched in tumor epithelium and stroma from 15 IBC and 35 non-IBC cases to study the relative contribution of each component to the IBC phenotype and to patient survival. '
type: 'Expression profiling by array'
contributor: 'Stefan,,Ambs...'
sample_id: 'GSM136326 GSM136327 GSM136328 GSM136329 GSM136330 GSM136331 GSM136332 GSM136333 GSM136334 GSM136335 GSM136336 GSM136337 GSM136338 GSM136339 GSM136340 GSM136341 GSM136342 GSM136343 GSM136344 GSM136345 GSM136346 GSM136347 GSM136348 GSM136349 GSM136350 GSM136351 GSM136352 GSM136353 GSM136354 GSM136355 GSM136356 GSM136357 GSM136358 GSM136359 GSM136360 GSM136361 GSM136362 GSM136363 GSM136364 GSM136365 GSM136366 GSM136367 GSM136368 GSM136369 GSM136370 GSM136371 GSM136372 GSM136373 GSM136374 GSM136375 GSM136376 GSM136377 GSM136378 GSM136379 GSM136380 GSM136381 GSM136382 GSM136383 GSM136384 GSM136385 GSM136386 GSM136387 GSM136388 GSM136389 GSM136390 GSM136391 GSM136392 GSM136393 GSM136394 GSM136395 GSM136396 GSM136397 GSM136398 GSM136399 GSM136400 GSM136401 GSM136402 GSM136403 GSM136404 GSM136405 GSM136406 GSM136407 GSM136408 GSM136409 GSM136410 GSM136411 GSM136412 GSM136413 GSM136414 GSM136415 GSM136416 GSM136417 GSM136418 GSM136419 GSM136420 '
contact_name: 'Stefan,,Ambs'
contact_laboratory: 'LHC'
contact_institute: 'NCI'
contact_address: '37 Convent Dr Bldg 37 Room 3050'
contact_city: 'Bethesda'
contact_state: 'MD'
contact_zip0x2Fpostal_code: '20892'
contact_country: 'USA'
supplementary_file: 'ftp://ftp.ncbi.nlm.nih.gov/pub/geo/DATA/supplementary/series/GSE5847/GSE5847_RAW.tar'
platform_id: 'GPL96'
platform_taxid: '9606'
sample_taxid: '9606'
relation: 'BioProject: http://www.ncbi.nlm.nih.gov/bioproject/97251'
gseData.Header.Samples
ans =
struct with fields:
title: {1x95 cell}
geo_accession: {1x95 cell}
status: {1x95 cell}
submission_date: {1x95 cell}
last_update_date: {1x95 cell}
type: {1x95 cell}
channel_count: {1x95 cell}
source_name_ch1: {1x95 cell}
organism_ch1: {1x95 cell}
characteristics_ch1: {2x95 cell}
molecule_ch1: {1x95 cell}
extract_protocol_ch1: {1x95 cell}
label_ch1: {1x95 cell}
label_protocol_ch1: {1x95 cell}
taxid_ch1: {1x95 cell}
hyb_protocol: {1x95 cell}
scan_protocol: {1x95 cell}
description: {1x95 cell}
data_processing: {1x95 cell}
platform_id: {1x95 cell}
contact_name: {1x95 cell}
contact_laboratory: {1x95 cell}
contact_institute: {1x95 cell}
contact_address: {1x95 cell}
contact_city: {1x95 cell}
contact_state: {1x95 cell}
contact_zip0x2Fpostal_code: {1x95 cell}
contact_country: {1x95 cell}
supplementary_file: {1x95 cell}
data_row_count: {1x95 cell}
data_processing содержит информацию о методах предварительной обработки, в данном случае процедуру RMA.
gseData.Header.Samples.data_processing(1)
ans =
1x1 cell array
{'RMA'}
source_name_ch1 поле содержит источник образца:
sampleSources = unique(gseData.Header.Samples.source_name_ch1);
sampleSources{:}
ans =
'human breast cancer stroma'
ans =
'human breast cancer tumor epithelium'
Область Header.Samples.characteristics_ch1 содержит характеристики образцов.
gseData.Header.Samples.characteristics_ch1(:,1)
ans =
2x1 cell array
{'IBC' }
{'Deceased'}
Определение меток КСГМГ и не-КСГМГ для образцов из Header.Samples.characteristics_ch1 поле в качестве меток группы.
sampleGrp = gseData.Header.Samples.characteristics_ch1(1,:);
Метаданные серии рассказали нам о платформе массива id: GPL96, которая представляет собой набор массивов Affymetrix ® GeneChip ® Human Genome U133 HG-U133A. Извлеките файл формата GPL96 SOFT из GEO с помощью getgeodata функция. Например, предположим, что вы использовали getgeodata функция, чтобы восстановить файл Платформы GPL96 и спасла его к файлу, такой как GPL96.txt. Используйте geosoftread для разбора этого файла формата SOFT.
gplData = geosoftread('GPL96.txt')
gplData =
struct with fields:
Scope: 'PLATFORM'
Accession: 'GPL96'
Header: [1x1 struct]
ColumnDescriptions: {16x1 cell}
ColumnNames: {16x1 cell}
Data: {22283x16 cell}
ColumnNames области gplData структура содержит имена столбцов для данных:
gplData.ColumnNames
ans =
16x1 cell array
{'ID' }
{'GB_ACC' }
{'SPOT_ID' }
{'Species Scientific Name' }
{'Annotation Date' }
{'Sequence Type' }
{'Sequence Source' }
{'Target Description' }
{'Representative Public ID' }
{'Gene Title' }
{'Gene Symbol' }
{'ENTREZ_GENE_ID' }
{'RefSeq Transcript ID' }
{'Gene Ontology Biological Process'}
{'Gene Ontology Cellular Component'}
{'Gene Ontology Molecular Function'}
Вы можете получить идентификаторы наборов зондов и символы генов для наборов зондов GPL69 платформы.
gplProbesetIDs = gplData.Data(:, strcmp(gplData.ColumnNames, 'ID')); geneSymbols = gplData.Data(:, strcmp(gplData.ColumnNames, 'Gene Symbol'));
Использование символов генов для маркировки генов в объекте DataMatrix gseData.Data. Имейте в виду, что зонд установлен IDs из файла платформы может быть не в том порядке, как в gseData.Data. В этом примере они находятся в одном порядке.
Измените имена строк данных экспрессии на символы генов.
gseData.Data = rownames(gseData.Data, ':', geneSymbols);
Отображение первых пяти строк и пяти столбцов данных экспрессии с именами строк в качестве символов генов.
gseData.Data(1:5, 1:5)
ans =
GSM136326 GSM136327 GSM136328 GSM136329 GSM136330
DDR1 10.45 9.3995 9.4248 9.4729 9.2788
RFC2 5.7195 4.8493 4.7321 4.7289 5.3264
HSPA6 5.9387 6.0833 6.448 6.1769 6.5446
PAX8 8.0231 7.8947 8.345 8.1632 8.2338
GUCA1A 3.9548 3.9632 3.9641 4.0878 3.9989
Инструментарий биоинформатики и Toolbox™ статистики и машинного обучения предлагают широкий спектр инструментов анализа и визуализации для анализа данных микрочипов. Однако, поскольку не является нашей основной целью объяснить методы анализа в этом примере, вы примените только несколько функций к профилю экспрессии генов из стромальных клеток. Более подробные примеры анализа экспрессии генов и доступных инструментов выбора признаков см. в разделах Изучение данных экспрессии генов микрочипов и Выбор функций для классификации высокоразмерных данных.
Эксперимент проводили на IBC и не-IBC образцах, полученных из стромальных клеток и эпителиальных клеток. В этом примере вы будете работать с профилем экспрессии генов из стромальных клеток. Определите индексы образцов для типа стромальных клеток из gseData.Header.Samples.source_name_ch1 поле.
stromaIdx = strcmpi(sampleSources{1}, gseData.Header.Samples.source_name_ch1);
Определите количество образцов из стромальных клеток.
nStroma = sum(stromaIdx)
nStroma =
47
stromaData = gseData.Data(:, stromaIdx); stromaGrp = sampleGrp(stromaIdx);
Определите количество образцов стромальных клеток IBC и без IBC.
nStromaIBC = sum(strcmp('IBC', stromaGrp))
nStromaIBC =
13
nStromaNonIBC = sum(strcmp('non-IBC', stromaGrp))
nStromaNonIBC =
34
Можно также пометить столбцы в stromaData с метками группы:
stromaData = colnames(stromaData, ':', stromaGrp);
Отображение гистограммы измерений нормализованной экспрессии гена указанного гена. Оси X представляют нормализованный уровень выражения. Например, проверьте распределение значений экспрессии генов этих генов.
fID = 331:339;
zValues = zscore(stromaData.(':')(':'), 0, 2); bw = 0.25; edges = -10:bw:10; bins = edges(1:end-1) + diff(edges)/2;
histStroma = histc(zValues(fID, :)', edges) ./ (stromaData.NCols*bw);
figure; for i = 1:length(fID) subplot(3,3,i); bar(edges, histStroma(:,i), 'histc') xlim([-3 3]) if i <= length(fID)-3 ax = gca; ax.XTickLabel = []; end title(sprintf('gene%d - %s', fID(i), stromaData.RowNames{fID(i)})) end sgtitle('Gene Expression Value Distributions')
Доступ к профилю экспрессии генов был получен с использованием Affymetrix GeneChip более чем 22000 признаков на небольшом количестве образцов (~ 100). Среди 47 образцов стромы опухоли есть 13 МКБ и 34 не-МКБ. Но не все гены дифференциально экспрессируются между опухолями IBC и не-IBC. Статистические тесты необходимы для идентификации сигнатуры экспрессии генов, которая отличает IBC от стромальных образцов без IBC.
Использовать genevarfilter для фильтрации генов с небольшой дисперсией по образцам.
[mask, stromaData] = genevarfilter(stromaData);
stromaData.NRows
ans =
20055
Примените t-статистику к каждому гену и сравните p-значения для каждого гена, чтобы найти значительно дифференциально экспрессируемые гены между группами IBC и не-IBC путем перестановки образцов (1000 раз для этого примера).
rng default [pvalues, tscores] = mattest(stromaData(:, 'IBC'), stromaData(:, 'non-IBC'),... 'Showhist', true', 'showplot', true, 'permute', 1000);
Выберите гены с заданным значением p.
sum(pvalues < 0.001)
ans =
52
Существует около 50 генов, отобранных непосредственно при значениях p < 0,001.
Сортируйте и перечисляйте 20 лучших генов:
testResults = [pvalues, tscores]; testResults = sortrows(testResults); testResults(1:20, :)
ans =
p-values t-scores
INPP5E 2.3318e-05 5.0389
ARFRP1 /// IGLJ3 2.7575e-05 4.9753
USP46 3.4336e-05 -4.9054
GOLGB1 4.7706e-05 -4.7928
TTC3 0.00010695 -4.5053
THUMPD1 0.00013164 -4.4317
0.00016042 4.3656
MAGED2 0.00017042 -4.3444
DNAJB9 0.0001782 -4.3266
KIF1C 0.00022122 4.2504
0.00022237 -4.2482
DZIP3 0.00022414 -4.2454
COPB1 0.00023199 -4.2332
PSD3 0.00024649 -4.2138
PLEKHA4 0.00026505 4.186
DNAJB9 0.0002767 -4.1708
CNPY2 0.0002801 -4.1672
USP9X 0.00028442 -4.1619
SEC22B 0.00030146 -4.1392
GFER 0.00030506 -4.1352
[1] Boersma, B.J., Reimers, M., Yi, M., Ludwig, J.A., et al. «Стромальная генная сигнатура, связанная с воспалительным раком молочной железы», International Journal of Cancer, 122 (6): 1324-32, 2008.