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

Данные в этом примере содержат данные подсчета синтетических генов для 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)

Identity строка представляет модель Пуассона, где отклонение идентично среднему значению, описанному в [3]. Обратите внимание, что данные, по-видимому, чрезмерно разбиты (то есть большинство точек выше Identity линия). Constant строка представляет отрицательную биномиальную модель, где дисперсией является сумма члена (среднего) шумов выстрела и константы, умноженной на среднее в квадрате, как описано в [2]. 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)