Этот пример показывает, как анализировать сводные данные Illumina BeadChip экспрессии гена с помощью функций Bioinformatics Toolbox™ и MATLAB®.
Этот пример показывает, как импортировать и анализировать данные об экспрессии гена из платформы Illumina BeadChip микромассивов. Набор данных в примере от исследования профилей экспрессии гена человеческого сперматогенеза Platts и др. [1]. Уровни экспрессии были измерены на Человеке Illumina Sentrix 6 BeadChips (WG6).
Данные из большинства экспериментов экспрессии гена микромассивов обычно состоят из четырех компонентов: экспериментируйте значения данных, демонстрационная информация, аннотации функции и информация об эксперименте. Вы будете работать с микроданными массива, создать каждый из этих четырех компонентов, собрать их в объект ExpressionSet и найти дифференцированно выраженные гены. Для большего количества примеров о классе ExpressionSet смотрите Работу с Объектами для Данных об Эксперименте Микромассивов.
Выборки были гибридизированы на трех Человеке Illumina Sentrix 6 BeadChips (WG6). Получите ГЕО Серийные данные GSE6967 с помощью функции 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 из файлов сводных данных. Значения экспрессии гена идентифицированы с целевыми идентификаторами зонда 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) на третьем чипе.
Используйте коробчатую диаграмму, чтобы просмотреть необработанные уровни экспрессии каждой выборки в эксперименте.
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 для попарного сравнения этих трех массивов на первом чипе, гибридизированном с teratozoospermic выборками (T2, T1 и T6).
Примечание: можно также использовать функцию mairplot, чтобы создать MA или IR (Интенсивность/Отношение) графики для сравнения определенных массивов.
inspectIdx = 1:3;
maxyplot(rawMatrix, inspectIdx)
suptitle('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)
suptitle('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 для Человеческих 6 BeadChips (WG6) Sentrix от Страницы поддержки на веб-сайте 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 можно собрать важную информацию о выборках, таких как демонстрационные заголовки, ГЕО демонстрационные инвентарные номера, и т.д., и сохранить выборочные данные как объект 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'
Набор метаданных функции для Человеческих 6 BeadChips Sentrix является большим и разнообразным. Выберите информацию о функциях, которые уникальны для эксперимента и сохраняют информацию как объект 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
Загрузите данные об эксперименте, сохраненные от предыдущего раздела. Вы будете использовать эти данные, чтобы найти дифференцированно выраженные гены между teratozoospermia и нормальными выборками.
load TNGeneExprSet
Чтобы идентифицировать дифференциальные изменения на уровнях расшифровок стенограммы в normospermic Ns и teratozoospermic выборках 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]. Оцените ФРГ и q-значения для каждого теста с помощью функции mafdr.
figure;
[pFDR, qValues] = mafdr(pValues, 'showplot', true);
diffScoresFDRQ = diffscore(qValues, meanTzData, meanNsData);
Определите количество генов с абсолютным дифференциальным счетом, больше, чем 20. Примечание: можно получить различное количество генов из-за теста перестановки и результатов начальной загрузки.
sum(abs(diffScoresFDRQ)>=20)
ans =
3121
Постройте-log10 p-значений против изменений сгиба в графике вулкана.
diffStruct = mavolcanoplot(TzData, NsData, qValues,... 'pcutoff', 0.01, 'foldchange', 5);
Примечание: Из графика вулкана пользовательский интерфейс, можно в интерактивном режиме изменить сокращение 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 к clustergram.
cmap = redbluecmap(9); cg = clustergram(TNDiffExprSet.expressions,'Colormap',cmap,'Standardize',2); addTitle(cg,'Hierarchical Gene Clustering')
Кластеризация наиболее дифференцированно богатых расшифровок стенограммы ясно разделы teratozoospermic (Tz) и normospermic (Ns) spermatozoal РНК.
[1] Platts, A.E., и др., "Успех и провал в человеческом сперматогенезе, как показано teratozoospermic РНК", Человеческая Молекулярная Генетика, 16 (7):763-73, 2007.
[2] Ярус, J.D. и Tibshirani, R., "Статистическое значение для genomewide исследований", PNAS, 100 (16):9440-5, 2003.