В этом примере показано, как анализировать сводные данные 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.