В этом примере показано, как анализировать сводные данные 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)
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 файл аннотации для Человеческих 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 =
3122
Постройте-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 к кластерграмме.
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.