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

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

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

load('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

Из 5000 генов 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)

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

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

plotChiSquaredFit(tLocal)

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

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

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

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

ans = 2×2

         493           7
           5        4495

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

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

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

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

plotChiSquaredFit(tPoisson)