Работа с данными серии GEO
Этот пример показывает, как извлечь серию данных экспрессии генов (GSE) из NCBI Gene Expression Omnibus (GEO) и выполнить базовый анализ профилей экспрессии.
Введение
NCBI Gene Expression Omnibus (GEO) является крупнейшим публичным хранилищем экспериментальных данных высокопроизводительных микрочипов. Данные GEO имеют четыре типа сущностей, включая GEO Platform (GPL), GEO Sample (GSM), GEO Series (GSE) и курированный GEO DataSet (GDS).
Запись Platform описывает список элементов массива в эксперименте (например, кДНК, олигонуклеотидные наборы проб). Каждой записи Платформы присваивается уникальный и стабильный номер присоединения к ГЕО (GPLxxx).
Запись Sample описывает условия, при которых обрабатывался индивидуум Выборки, манипуляции, которым он подвергся, и измерение численности каждого выведенного из него элемента. Каждой выборочной записи присваивается уникальный и стабильный номер присоединения GEO (GSMxxx).
Запись Series определяет группу связанных выборок и предоставляет координационную точку и описание всего исследования. Последовательные записи могут также содержать таблицы, описывающие извлеченные данные, итоговые выводы или анализы. Каждой записи серии присваивается уникальный номер присоединения GEO (GSExxx).
Запись DataSet (GDSxxx) представляет собой кураторский набор биологически и статистически сопоставимых выборок GEO. Наборы данных GEO (GDSxxx) являются кураторскими наборами выборочных данных GEO.
Более подробную информацию о GEO можно найти в обзоре GEO. Bioinformatics Toolbox™ обеспечивает функции, которые могут извлекать и анализировать файлы данных формата GEO. Данные GSE, GSM, GSD и GPL могут быть получены при помощи getgeodata функция. The getgeodata функция может также сохранять извлеченные данные в текстовом файле. Записи GEO Series доступны в файлах формата SOFT и в файлах текстового формата с разделителями табуляций. Функция geoseriesread считывает текстовый файл формата GEO Series. The geosoftread функция считывает обычно довольно большие файлы формата SOFT.
В этом примере будет получен набор GSE5847 данных из базы данных GEO и выполнена статистическая проверка данных. GEO Series GSE5847 содержит экспериментальные данные исследования экспрессии генов опухолевой стромы и камер из 15 случаев воспалительного рака молочной железы (IBC) и 35 случаев невоспалительного рака молочной железы (Boersma et al. 2008).
Получение данных серии GEO
Функция 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 поле, в котором хранятся матричные данные Series.
Исследование данных GSE
Матричные данные 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
Метаданные Series хранятся в Header поле. Структура содержит информацию Series в 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}
The data_processing поле содержит информацию о методах предварительной обработки, в этом случае о процедуре Robust Мультимассива Среднего значения (RMA).
gseData.Header.Samples.data_processing(1)
ans =
1x1 cell array
{'RMA'}
The 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'}
Определите метки IBC и не-IBC для выборок из Header.Samples.characteristics_ch1 поле как метки группы.
sampleGrp = gseData.Header.Samples.characteristics_ch1(1,:);
Получение данных платформы GEO (GPL)
Метаданные Series сообщили нам, что платформа массива id: GPL96, которая является набором U133 массива Affymetrix ® GeneChip ® Human Genome HG-U133A. Загрузите GPL96 файл формата SOFT из GEO с помощью getgeodata функция. Например, предположим, что вы использовали getgeodata функция для извлечения файла GPL96 Platform и сохранения его в файл, например 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}
The 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
Анализ данных
Bioinformatics Toolbox и Statistics and Machine Learning 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 IBC и 34 не-IBC. Но не все гены дифференциально экспрессируются между опухолями 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.