В этом примере показано, как анализировать сводные данные экспрессии генов Illumina BeadChip с использованием функций MATLAB ® и Bioinformatics Toolbox™.
В этом примере показано, как импортировать и анализировать данные экспрессии генов с платформы микрочипов Illumina BeadChip. Данные, приведенные в этом примере, взяты из исследования профилей экспрессии генов сперматогенеза человека Platts et al. [1]. Уровни экспрессии измеряли на Illumina Sentrix Human 6 (WG6) BeadChips.
Данные большинства экспериментов по экспрессии генов микрочипов, как правило, состоят из четырех компонентов: значения данных эксперимента, информация о пробе, пояснения к особенностям и информация об эксперименте. Вы будете работать с данными микрочипов, создавать каждый из четырех компонентов, собирать их в ExpressionSet объект и найти дифференциально экспрессируемые гены. Дополнительные примеры о ExpressionSet см. раздел Работа с объектами для данных эксперимента микрочипов.
Образцы гибридизовали на трех Illumina Sentrix Human 6 (WG6) BeadChips. Извлечение GSE6967 данных серии GEO с помощью getgeodata функция.
TNGEOData = getgeodata('GSE6967')
TNGEOData =
struct with fields:
Header: [1×1 struct]
Data: [47293×13 bioma.data.DataMatrix]
TNGEOData структура содержит Header и Data поля. Header поле имеет два поля, Series и Samples, содержащий описание эксперимента и образцов соответственно. 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 из файлов сводных данных. Значения экспрессии генов идентифицируют с помощью ID-мишеней зонда 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 BeadChip в 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 - это график рассеяния одного массива относительно другого. В этом примере используется вспомогательная функция maxyplot построить графики MAXY для попарного сравнения трех массивов на первом чипе, гибридизированном с тератозооспермическими образцами (T2, T1 и Т6).
Примечание: Вы также можете использовать mairplot создание графиков MA или IR (интенсивность/отношение) для сравнения конкретных массивов.
inspectIdx = 1:3;
maxyplot(rawMatrix, inspectIdx)
sgtitle('Before Normalization')
В графике MAXY графики MA для всех попарных сравнений находятся в правом верхнем углу. В левом нижнем углу находятся графики XY сравнений. График MAXY показывает, что два массива, T1 и T2, очень похожи, в то время как отличаются от другого массива, T6.
Графики поля экспрессии и графики 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)'
Примечание: 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 в 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. Примечание: Вы можете получить другое количество генов из-за теста перестановки и результатов начальной загрузки.
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) сперматозоидные РНК.
Platts, A.E., et al., «Успех и неудача в сперматогенезе человека, выявленные тератозооспермическими РНК», Human Molecular Genetics, 16 (7): 763-73, 2007.
[2] Стори, Дж. Д. и Тибширани, Р., «Статистическая значимость для исследований по всему геному», PNAS, 100 (16): 9440-5, 2003.