Этот пример показывает, как анализировать сводные данные по экспрессии гена Illumina BeadChip с помощью функций MATLAB ® и Bioinformatics Toolbox™.
В этом примере показано, как импортировать и проанализировать данные экспрессии генов с платформы микромассивов Illumina BeadChip. Данные, представленные в этом примере, взяты из исследования профилей экспрессии генов сперматогенеза человека Platts et al. [1]. Уровни экспрессии измеряли на Illumina Sentrix Human 6 (WG6) BeadChips.
Данные большинства экспериментов по экспрессии генов микромассивов обычно состоят из четырех компонентов: значения данных эксперимента, информация о выборке, аннотации функций и информация об эксперименте. Вы будете работать с данными микромассивов, создавать каждый из четырех компонентов, собирать их в ExpressionSet объект, и найти дифференциально экспрессированные гены. Для получения дополнительных примеров о ExpressionSet класс см. «Работа с объектами для данных эксперимента с микромассивами».
Выборки гибридизовали с тремя Human 6 (WG6) BeadChips Illumina Sentrix. Получите данные серии GEO GSE6967 используя getgeodata функция.
TNGEOData = getgeodata('GSE6967')
TNGEOData =
struct with fields:
Header: [1×1 struct]
Data: [47293×13 bioma.data.DataMatrix]
The TNGEOData структура содержит Header и Data поля. The Header поле имеет два поля, Series и Samples, содержащего описание эксперимента и выборок соответственно. The Data поле содержит DataMatrix нормализованных суммарных уровней экспрессии из эксперимента.
Определите количество выборок в эксперименте.
nSamples = numel(TNGEOData.Header.Samples.geo_accession)
nSamples =
13
Смотрите заголовки образцов из Header.Samples поле.
TNGEOData.Header.Samples.title'
ans =
13×1 cell array
{'Teratozoospermic individual: Sample T2'}
{'Teratozoospermic individual: Sample T1'}
{'Teratozoospermic individual: Sample T6'}
{'Teratozoospermic individual: Sample T4'}
{'Teratozoospermic individual: Sample T8'}
{'Normospermic individual: Sample N11' }
{'Teratozoospermic individual: Sample T3'}
{'Teratozoospermic individual: Sample T7'}
{'Teratozoospermic individual: Sample T5'}
{'Normospermic individual: Sample N6' }
{'Normospermic individual: Sample N12' }
{'Normospermic individual: Sample N5' }
{'Normospermic individual: Sample N1' }
Для простоты извлеките метки выборок из заголовков образцов.
sampleLabels = cellfun(@(x) char(regexp(x, '\w\d+', 'match')),... TNGEOData.Header.Samples.title, 'UniformOutput',false)
sampleLabels =
1×13 cell array
Columns 1 through 7
{'T2'} {'T1'} {'T6'} {'T4'} {'T8'} {'N11'} {'T3'}
Columns 8 through 13
{'T7'} {'T5'} {'N6'} {'N12'} {'N5'} {'N1'}
Загрузите дополнительный файл GSE6967_RAW.tar и разархивируйте файл, чтобы получить доступ к 13 текстовым файлам, созданным программным обеспечением BeadStudio, которые содержат необработанные, ненормализованные суммарные значения шариков.
Необработанные файлы данных называются с их номерами доступа GSM. В данном примере создайте имена файлов текстовых файлов данных с помощью пути, в котором расположены текстовые файлы.
rawDataFiles = cell(1,nSamples); for i = 1:nSamples rawDataFiles {i} = [TNGEOData.Header.Samples.geo_accession{i} '.txt']; end
Установите переменную rawDataPath в путь и директорию, в который вы извлекли файлы данных.
rawDataPath = 'C:\Examples\illuminagedemo\GSE6967';
Используйте ilmnbsread функция для чтения первого из сводных файлов и хранения содержимого в структуре.
rawData = ilmnbsread(fullfile(rawDataPath, rawDataFiles{1}))
rawData =
struct with fields:
Header: [1×1 struct]
TargetID: {47293×1 cell}
ColumnNames: {1×8 cell}
Data: [47293×8 double]
TextColumnNames: {}
TextData: {}
Просмотрите имена столбцов в rawData структура.
rawData.ColumnNames'
ans =
8×1 cell array
{'MIN_Signal-1412091085_A' }
{'AVG_Signal-1412091085_A' }
{'MAX_Signal-1412091085_A' }
{'NARRAYS-1412091085_A' }
{'ARRAY_STDEV-1412091085_A'}
{'BEAD_STDEV-1412091085_A' }
{'Avg_NBEADS-1412091085_A' }
{'Detection-1412091085_A' }
Определите количество целевых зондов.
nTargets = size(rawData.Data, 1)
nTargets =
47293
Считайте ненормализованные значения выражения (Avg_Signal), пределы доверия обнаружения и идентификаторы микросхем Sentrix из файлов сводных данных. Значения экспрессии генов идентифицируют с идентификаторов зонда Illumina. Можно задать столбцы для чтения из файла данных.
rawMatrix = bioma.data.DataMatrix(zeros(nTargets, nSamples),... rawData.TargetID, sampleLabels); detectionConf = bioma.data.DataMatrix(zeros(nTargets, nSamples),... rawData.TargetID, sampleLabels); chipIDs = categorical([]); for i = 1:nSamples rawData =ilmnbsread(fullfile(rawDataPath, rawDataFiles{i}),... 'COLUMNS', {'AVG_Signal', 'Detection'}); chipIDs(i) = regexp(rawData.ColumnNames(1), '\d*', 'match', 'once'); rawMatrix(:, i) = rawData.Data(:, 1); detectionConf(:,i) = rawData.Data(:,2); end
В эксперименте используются три Sentrix BeadChips. Смотрите идентификаторы бисера Illumina Sentrix в chipIDs для определения количества выборок, гибридизованных на каждом чипе.
summary(chipIDs) samplesChip1 = sampleLabels(chipIDs=='1412091085') samplesChip2 = sampleLabels(chipIDs=='1412091086') samplesChip3 = sampleLabels(chipIDs=='1477791158')
1412091085 1412091086 1477791158
6 4 3
samplesChip1 =
1×6 cell array
{'T2'} {'T1'} {'T6'} {'T4'} {'T8'} {'N11'}
samplesChip2 =
1×4 cell array
{'T3'} {'T7'} {'T5'} {'N6'}
samplesChip3 =
1×3 cell array
{'N12'} {'N5'} {'N1'}
Шесть выборок (T2, T1, T6, T4, T8 и N11) были скрещены к шести массивам на первом чипе, четыре выборки (T3, T7, T5 и N6) на втором чипе и трех выборках (N12, N5 и N1) на третьем чипе.
Используйте boxplot, чтобы просмотреть необработанные уровни экспрессии каждой выборки в эксперименте.
logRawExprs = log2(rawMatrix); maboxplot(logRawExprs, 'Orientation', 'horizontal') ylabel('Sample Labels') xlabel('log2(Expression Value)') title('Before Normalization')
Различие интенсивности между выборками на одном чипе и выборки на разных чипах не кажется слишком большим. Первый BeadChip, содержа выборки T2, T1, T6, T4, T8 и N11, кажется, немного больше переменной, чем другие.
Использование графиков MA и XY для парного сравнения массивов на BeadChip может быть информативным. На графике MA среднее значение (A) уровней выражения двух массивов строится на оси x и различие (M) в измерении на оси y. График XY является графиком поля точек одного массива против другого. Этот пример использует функцию helper maxyplot построение графиков MAXY для парного сравнения трех массивов на первом чипе, гибридизованных с тератозооспермическими выборками (T2, T1 и T6).
Примечание: Вы также можете использовать mairplot функция для создания графиков MA или IR (Intensity/Ratio) для сравнения конкретных массивов.
inspectIdx = 1:3;
maxyplot(rawMatrix, inspectIdx)
sgtitle('Before Normalization')
На графике MAXY графики MA для всех попарных сравнений находятся в правом верхнем углу. В левом нижнем углу расположены графики XY сравнений. График MAXY показывает, что два массива, T1 и T2, довольно похожи, в то время как другие массивы, T6.
Графики выражения box и графики MAXY показывают, что существуют различия в уровнях экспрессии внутри чипов и между чипами; следовательно, данные требуют нормализации. Используйте quantilenorm функция для применения квантовой нормализации к необработанным данным.
Примечание: Вы также можете попробовать инвариантную нормализацию набора с помощью mainvarsetnorm функция.
normExprs = rawMatrix; normExprs(:, :) = quantilenorm(rawMatrix.(':')(':')); log2NormExprs = log2(normExprs);
Отобразите и проверьте нормированные уровни выражения на прямоугольном графике.
figure; maboxplot(log2NormExprs, 'ORIENTATION', 'horizontal') ylabel('Sample Labels') xlabel('log2(Expression Value)') title('After Quantile Normalization')
Отображение и просмотр графика MAXY трех массивов (T2, T1 и T6) на первом чипе после нормализации.
maxyplot(normExprs, inspectIdx)
sgtitle('After Quantile Normalization')
Многие гены в этом исследовании не экспрессируются или имеют только небольшую изменчивость по выборкам.
Во-первых, можно удалить гены с очень низкими абсолютными значениями экспрессии при помощи genelowvalfilter.
[mask, log2NormExprs] = genelowvalfilter(log2NormExprs); detectionConf = detectionConf(mask, :);
Во-вторых, отфильтровывайте гены с небольшим отклонением между выборками, используя genevarfilter.
[mask, log2NormExprs] = genevarfilter(log2NormExprs); detectionConf = detectionConf(mask, :);
Производители микромассивов обычно предоставляют аннотации набора функций для каждого типа чипа. Файлы аннотации чипа содержат такие метаданные, как имя гена, символ, номер присоединения NCBI, расположение хромосомы и информация о пути. Перед сборкой ExpressionSet объект для эксперимента получает аннотации о функциях или зондах на BeadChip. Вы можете скачать Human_WG-6.csv файл аннотации для Sentrix Human 6 (WG6) BeadChips со страницы поддержки на веб-сайте Illumina и сохраните файл локально. Считайте файл аннотации в MATLAB как dataset массив. Установите переменную annotPath в путь и директорию, в который вы загрузили файл аннотации.
annotPath = fullfile('C:\Examples\illuminagedemo\Annotation'); WG6Annot = dataset('xlsfile', fullfile(annotPath, 'Human_WG-6.csv'));
Просмотр свойств этого dataset массив.
get(WG6Annot)
Description: ''
VarDescription: {}
Units: {}
DimNames: {'Observations' 'Variables'}
UserData: []
ObsNames: {}
VarNames: {1×13 cell}
Получите имена переменных, описывающих функции чипа.
fDataVariables = get(WG6Annot, 'VarNames')
fDataVariables =
1×13 cell array
Columns 1 through 5
{'Search_key'} {'Target'} {'ProbeId'} {'Gid'} {'Transcript'}
Columns 6 through 10
{'Accession'} {'Symbol'} {'Type'} {'Start'} {'Probe_Sequence'}
Columns 11 through 13
{'Definition'} {'Ontology'} {'Synonym'}
Проверьте количество целевых идентификаторов зондов в файле аннотации.
numel(WG6Annot.Target)
ans =
47296
Поскольку данные выражения в этом примере являются только небольшим набором полных значений выражения, вы будете работать только с функциями в DataMatrix log2NormExprs объекта. Найдите соответствующие функции в log2NormExprs и WG6Annot.Target.
[commTargets, fI, WGI] = intersect(rownames(log2NormExprs), WG6Annot.Target);
Можно хранить предварительно обработанные значения выражения и пределы обнаружения аннотированных целей зонда как ExptData объект.
fNames = rownames(log2NormExprs);
TNExptData = bioma.data.ExptData({log2NormExprs(fI, :), 'ExprsValues'},...
{detectionConf(fI, :), 'DetectionConfidences'})
TNExptData = Experiment Data: 42313 features, 13 samples 2 elements Element names: ExprsValues, DetectionConfidences
Выборочные данные в Header.Samples поле TNGEOData структура может быть подавляющей и трудно ориентироваться. Из данных в Header.Samples поле, можно собрать необходимую информацию о выборках, таких как заголовки образцов, номера присоединения к GEO и т.д., и сохранить выборочные данные как MetaData объект.
Найдите описание выборочных характеристик.
sampleChars = cellfun(@(x) char(regexp(x, '\w*tile', 'match')),... TNGEOData.Header.Samples.characteristics_ch1, 'UniformOutput', false)
sampleChars =
1×13 cell array
Columns 1 through 4
{'Infertile'} {'Infertile'} {'Infertile'} {'Infertile'}
Columns 5 through 8
{'Infertile'} {'Fertile'} {'Infertile'} {'Infertile'}
Columns 9 through 13
{'Infertile'} {'Fertile'} {'Fertile'} {'Fertile'} {'Fertile'}
Создайте dataset массив для хранения только что извлеченных выборочных данных.
sampleDS = dataset({TNGEOData.Header.Samples.geo_accession', 'GSM'},...
{strtok(TNGEOData.Header.Samples.title)', 'Type'},...
{sampleChars', 'Characteristics'}, 'obsnames', sampleLabels')
sampleDS =
GSM Type Characteristics
T2 {'GSM160620'} {'Teratozoospermic'} {'Infertile'}
T1 {'GSM160621'} {'Teratozoospermic'} {'Infertile'}
T6 {'GSM160622'} {'Teratozoospermic'} {'Infertile'}
T4 {'GSM160623'} {'Teratozoospermic'} {'Infertile'}
T8 {'GSM160624'} {'Teratozoospermic'} {'Infertile'}
N11 {'GSM160625'} {'Normospermic' } {'Fertile' }
T3 {'GSM160626'} {'Teratozoospermic'} {'Infertile'}
T7 {'GSM160627'} {'Teratozoospermic'} {'Infertile'}
T5 {'GSM160628'} {'Teratozoospermic'} {'Infertile'}
N6 {'GSM160629'} {'Normospermic' } {'Fertile' }
N12 {'GSM160630'} {'Normospermic' } {'Fertile' }
N5 {'GSM160631'} {'Normospermic' } {'Fertile' }
N1 {'GSM160632'} {'Normospermic' } {'Fertile' }
Сохраните выборочные метаданные как объект MetaData класс, включая краткое описание для каждой переменной.
TNSData = bioma.data.MetaData(sampleDS,... {'Sample GEO accession number',... 'Spermic type',... 'Fertility characteristics'})
TNSData =
Sample Names:
T2, T1, ...,N1 (13 total)
Variable Names and Meta Information:
VariableDescription
GSM {'Sample GEO accession number'}
Type {'Spermic type' }
Characteristics {'Fertility characteristics' }
Набор метаданных функций для Sentrix Human 6 BeadChips велик и разнообразен. Выберите информацию о функциях, уникальных для эксперимента, и сохраните информацию как MetaData объект. Извлеките интересующие аннотации, например Accession и Symbol.
fIdx = ismember(fDataVariables, {'Accession', 'Symbol'});
featureAnnot = WG6Annot(WGI, fDataVariables(fIdx));
featureAnnot = set(featureAnnot, 'ObsNames', WG6Annot.Target(WGI));
Создайте MetaData объект для информации аннотации функции с кратким описанием двух переменных метаданных.
WG6FData = bioma.data.MetaData(featureAnnot, ... {'Accession number of probe target', 'Gene Symbol of probe target'})
WG6FData =
Sample Names:
GI_10047089-S, GI_10047091-S, ...,hmm9988-S (42313 total)
Variable Names and Meta Information:
VariableDescription
Accession {'Accession number of probe target'}
Symbol {'Gene Symbol of probe target' }
Большинство описаний экспериментов в Header.Series поле TNGEOData структура может быть реорганизована и сохранена как MIAME объект, который вы будете использовать, чтобы собрать ExpressionSet объект для эксперимента.
TNExptInfo = bioma.data.MIAME(TNGEOData)
TNExptInfo =
Experiment Description:
Author name: Adrian,E,Platts
David,J,Dix
Hector,E,Chemes
Kary,E,Thompson
Robert,,Goodrich
John,C,Rockett
Vanesa,Y,Rawe
Silvina,,Quintana
Michael,P,Diamond
Lillian,F,Strader
Stephen,A,Krawetz
Laboratory: Wayne State University
Contact information: Stephen,A,Krawetz
URL: http://compbio.med.wayne.edu
PubMedIDs: 17327269
Abstract: A 82 word abstract is available. Use the Abstract property.
Experiment Design: A 61 word summary is available. Use the ExptDesign property.
Other notes:
{'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE6nnn/GSE6967/suppl/GSE6967_RAW.tar'}
Теперь, когда вы создали все компоненты, можно создать объект ExpressionSet класс для хранения значений выражения, информации о выборке, аннотаций функций микросхемы и описательной информации об этом эксперименте.
TNExprSet = bioma.ExpressionSet(TNExptData, 'sData', TNSData,... 'fData', WG6FData,... 'eInfo', TNExptInfo)
TNExprSet =
ExpressionSet
Experiment Data: 42313 features, 13 samples
Element names: Expressions, DetectionConfidences
Sample Data:
Sample names: T2, T1, ...,N1 (13 total)
Sample variable names and meta information:
GSM: Sample GEO accession number
Type: Spermic type
Characteristics: Fertility characteristics
Feature Data:
Feature names: GI_10047089-S, GI_10047091-S, ...,hmm9988-S (42313 total)
Feature variable names and meta information:
Accession: Accession number of probe target
Symbol: Gene Symbol of probe target
Experiment Information: use 'exptInfo(obj)'
Примечание: The ExprsValues элемент в ExptData объект, TNExptData, переименован в Expressions в TNGeneExprSet.
Можно сохранить объект ExpressionSet класс как MAT файл для дальнейшего анализа данных.
save TNGeneExprSet TNExprSet
Загрузите данные эксперимента, сохраненные из предыдущего раздела. Вы будете использовать эти данные, чтобы найти дифференциально экспрессированные гены между тератозооспермией и нормальными выборками.
load TNGeneExprSet
Идентифицировать дифференциальные изменения уровней транскриптов в нормоспермических Ns и тератозооспермическая Tz выборки, сравните значения экспрессии генов между двумя группами данных: Tz и Ns.
TNSamples = sampleNames(TNExprSet); Tz = strncmp(TNSamples, 'T', 1); Ns = strncmp(TNSamples, 'N', 1); nTz = sum(Tz) nNs = sum(Ns) TNExprs = expressions(TNExprSet); TzData = TNExprs(:,Tz); NsData = TNExprs(:,Ns); meanTzData = mean(TzData,2); meanNsData = mean(NsData,2); groupLabels = [TNSamples(Tz), TNSamples(Ns)];
nTz =
8
nNs =
5
Выполните сочетание t-тест с помощью mattest функция для транспозиции столбцов матрицы данных экспрессии генов TzData и NsData. Примечание: В зависимости от размера выборки, может быть нецелесообразно учитывать все возможные сочетания. Обычно случайное подмножество сочетаний рассматривается в случае большого размера выборки.
Используйте nchoosek функция в Statistics and Machine Learning Toolbox™, чтобы определить количество всех возможных сочетаний выборок в этом примере.
perms = nchoosek(1:nTz+nNs, nTz); nPerms = size(perms,1)
nPerms =
1287
Используйте PERMUTE опция mattest функция для вычисления значений p всех сочетаний.
pValues = mattest(TzData, NsData, 'Permute', nPerms);
Вы также можете вычислить дифференциальный счет от значений p с помощью следующей анонимной функции [1].
diffscore = @(p, avgTest, avgRef) -10*sign(avgTest - avgRef).*log10(p);
Дифференциальный счет 13 соответствует p-значению 0,05, дифференциальный счет 20 соответствует p-значению 0,01 и дифференциальный счет 30 соответствует p-значению 0,001. Положительный дифференциальный счет представляет повышающую регуляцию, а отрицательная счет понижающую регуляцию генов.
diffScores = diffscore(pValues, meanTzData, meanNsData);
Определите количество генов, которые, как считается, имеют дифференциальный счет, превышающую 20. Примечание: Вы можете получить разное количество генов из-за результата теста сочетания.
up = sum(diffScores > 20) down = sum(diffScores < -20)
up =
3741
down =
3033
В нескольких проверки гипотез, где мы одновременно тестируем нулевую гипотезу тысяч генов, каждый тест имеет специфическую ложноположительную частоту или частоту ложных открытий (FDR) [2]. Оцените FDR и q-значения для каждого теста с помощью mafdr функция.
figure;
[pFDR, qValues] = mafdr(pValues, 'showplot', true);
diffScoresFDRQ = diffscore(qValues, meanTzData, meanNsData);
Определите количество генов с абсолютным дифференциальным счетом более 20. Примечание: Вы можете получить разное количество генов из-за теста сочетания и результатов bootstrap.
sum(abs(diffScoresFDRQ)>=20)
ans =
3122
Постройте график -log10 значений p от изменений складок на графике вулкана.
diffStruct = mavolcanoplot(TzData, NsData, qValues,... 'pcutoff', 0.01, 'foldchange', 5);
Примечание: Из графика вулкана UI, вы можете в интерактивном режиме изменить предел отсечения и складывания значения p и экспортировать дифференциально экспрессированные гены.
Определите количество дифференциально экспрессированных генов.
nDiffGenes = numel(diffStruct.GeneLabels)
nDiffGenes = 451
Получите список регулируемых генов для Tz выборки по сравнению с Ns выборки.
up_genes = diffStruct.GeneLabels(diffStruct.FoldChanges > 0); nUpGenes = length(up_genes)
nUpGenes = 223
Получите список регулируемых вниз генов для Tz выборки по сравнению с Ns выборки.
down_genes = diffStruct.GeneLabels(diffStruct.FoldChanges < 0); nDownGenes = length(down_genes)
nDownGenes = 228
Извлеките список дифференциально экспрессированных генов.
diff_geneidx = zeros(nDiffGenes, 1); for i = 1:nDiffGenes diff_geneidx(i) = find(strncmpi(TNExprSet.featureNames, ... diffStruct.GeneLabels{i}, length(diffStruct.GeneLabels{i})), 1); end
Вы можете получить подмножество данных эксперимента, содержащих только дифференциально экспрессированные гены.
TNDiffExprSet = TNExprSet(diff_geneidx, groupLabels);
Анализ основного компонента (PCA) на дифференциально экспрессируемых генах показывает линейную разделяемость Tz выборки из Ns выборки.
PCAScore = pca(TNDiffExprSet.expressions);
Отобразите коэффициенты первого и шестого главных компонентов.
figure; plot(PCAScore(:,1), PCAScore(:,6), 's', 'MarkerSize',10, 'MarkerFaceColor','g'); hold on text(PCAScore(:,1)+0.02, PCAScore(:,6), TNDiffExprSet.sampleNames) plot([0,0], [-0.5 0.5], '--r') ax = gca; ax.XTick = []; ax.YTick = []; ax.YTickLabel = []; title('PCA Mapping') xlabel('Principal Component 1') ylabel('Principal Component 6')
Можно также использовать интерактивный инструмент, созданный mapcaplot функция для выполнения анализа основного компонента.
mapcaplot((TNDiffExprSet.expressions)')
Выполните неконтролируемую иерархическую кластеризацию значимых профилей генов из Tz и Ns группы, использующие корреляцию в качестве метрики расстояния для кластеризации выборок.
sampleDist = pdist(TNDiffExprSet.expressions','correlation'); sampleLink = linkage(sampleDist); figure; dendrogram(sampleLink, 'labels', TNDiffExprSet.sampleNames,'ColorThreshold',0.5) ax = gca; ax.YTick = []; ax.Box = 'on'; title('Hierarchical Sample Clustering')
Используйте clustergram функция для создания иерархической кластеризации дифференциально экспрессированных генов и применения палитры redbluecmap в кластерграмму.
cmap = redbluecmap(9); cg = clustergram(TNDiffExprSet.expressions,'Colormap',cmap,'Standardize',2); addTitle(cg,'Hierarchical Gene Clustering')
Кластеризация наиболее дифференциально распространенных транскриптов явно разделяет тератозооспермические (Tz) и нормоспермический (Ns) сперматозоальные РНК.
[1] Platts, A.E., et al., «Success and отказа in human spermatogenesis as showed by teratozoospermic RNAs», Human Molecular Genetics, 16 (7): 763-73, 2007.
[2] Storey, J.D. and Tibshirani, R., «Statistical valuance for genomewide studies», PNAS, 100 (16): 9440-5, 2003.