Этот пример показывает, как выполнить непарный тест гипотезы для синтетических данных о количестве короткого чтения из двух различных биологических условий.

Данные в этом примере содержат синтетические генные данные о количестве для 5 000 генов, представляя два различных биологических условия, такие как больные и нормальные ячейки. Для каждого условия существует пять выборок. Только 10% генов (500 генов) дифференцированно выражаются. А именно, половина из них (250 генов) точно 3-кратная сверхвыраженный. Другие 250 генов являются 3-кратным underexpressed. Остальная часть данных об экспрессии гена сгенерирована от того же отрицательного биномиального распределения для обоих условий. Каждая выборка также имеет различный фактор размера (то есть, покрытие или глубина выборки).

Загрузите данные.

clear all
load(fullfile(matlabroot,'examples','bioinfo','nbintest_data.mat'),'K','H0');

Переменная K содержит генные данные о количестве. Строки представляют гены, и столбцы представляют выборки. В этом случае первые пять столбцов представляют выборки от первого условия. Другие пять столбцов представляют выборки от второго условия. Отобразите первые несколько строк K.

K(1:5,:)

ans = 5×10

       13683       14140        8281       14309       12208        8045        9446       11317       14597       14592
       16028       16805        9813       16486       14076        9901       10927       13348       16999       17036
         814         862         492         910         758         521         573         753         870         936
       15870       16453        9857       16454       14267        9671       10997       13624       17151       17205
        9422        9393        5734        9598        8174        5381        6315        7752        9869        9795

В этом примере нулевая гипотеза верна, когда гена дифференцированно не выражают. Переменная H0 содержит булевы индикаторы, которые указывают, для которых генов нулевая гипотеза верна (отмеченный как 1). Другими словами, H0 содержит известные метки, которые вы будете использовать позже, чтобы соответствовать предсказанным результатам.

sum(H0)

ans = 4500

Из 5 000 генов, 4500 дифференцированно не выражаются в этих синтетических данных.

Запустите непарный тест гипотезы для выборок от двух условий с помощью nbintest. Предположение - то, что данные прибыли из отрицательного биномиального распределения, где отклонение соединяется со средним значением через локально регрессировавшую сглаженную функцию среднего значения, как описано в [1] установкой 'VarianceLink' к 'LocalRegression'.

tLocal = nbintest(K(:,1:5),K(:,6:10),'VarianceLink','LocalRegression');

Используйте plotVarianceLink, чтобы построить график рассеивания для каждого экспериментального условия (для условий X и Y) с демонстрационным отклонением в общей шкале по сравнению с оценкой зависимого условием среднего значения. Используйте линейную шкалу для обеих осей. Включайте кривые для всех других опций связи установкой 'Compare' к true.

plotVarianceLink(tLocal,'Scale','linear','Compare',true)

Строка Identity представляет модель Poisson, где отклонение идентично среднему значению, как описано в [3]. Заметьте, что данные, кажется, сверхрассеиваются (то есть, большинство точек выше строки Identity). Строка Constant представляет отрицательную биномиальную модель, где отклонение является суммой термина шума выстрела (среднее значение) и константа, умноженная на среднее значение в квадрате, как описано в [2]. Local Regression и опции связи Constant, кажется, соответствуют лучше сверхрассеянным данным.

Используйте plotChiSquaredFit, чтобы оценить качество подгонки для регрессии отклонения. Это строит эмпирический CDF (ecdf) вероятностей в квадрате хи. Вероятности являются отношением между наблюдаемым и предполагаемым отклонением, расслоенным уровнями количества короткого чтения в пять интервалов равного размера.

plotChiSquaredFit(tLocal)

Каждые данные показывают пять кривых ecdf. Каждая кривая представляет один из пяти уровней количества короткого чтения. Например, синяя строка представляет кривую ecdf для низких количеств короткого чтения между 0 и 1264. Красная линия представляет высокие количества (больше чем 11 438).

Один способ интерпретировать кривые состоит в том, чтобы проверять, ли кривые ecdf выше диагональной строки. Если они выше строки, то отклонение завышено. Если они ниже строки, то отклонение недооценено. В обеих фигурах отклонение, кажется, правильно оценивается для более высоких количеств (то есть, красная линия следует за диагональной строкой), но немного завышенный для более низких уровней количества.

Чтобы оценить производительность теста гипотезы, создайте матрицу беспорядка использование известных меток и предсказанных p-значений.

confusionmat(H0,(tLocal.pValue > .001))

ans = 2×2

         493           7
           5        4495

Из 500 дифференцированно выраженных генов, 493 правильно предсказаны (истинные положительные стороны), и 7 из них неправильно предсказаны как не дифференцированно выраженные гены (ложные отрицательные стороны). Из 4 500 генов, которых дифференцированно не выражают, 4495, правильно предсказаны (истинные отрицательные стороны), и 5 из них неправильно предсказаны как дифференцированно выраженные гены (ложные положительные стороны).

Для сравнения, запущенного тест гипотезы, снова принимающий, что количества моделируются распределением Пуассона, где отклонение идентично среднему значению.

tPoisson = nbintest(K(:,1:5),K(:,6:10),'VarianceLink','Identity');

Постройте кривые ecdf. Заметьте, что все кривые ниже диагональной строки, подразумевая, что отклонение недооценено. Поэтому отрицательная биномиальная модель соответствует данным лучше.

plotChiSquaredFit(tPoisson)