Кратность прочтения при 30% ошибке секвенирования

sunny82

Возник праздный вопрос по тестируемому сейчас в мире новому типу секвенаторов MiniIon
http://www.nanoporetech.com/technology/the-minion-device-a-...
Ситуация - если секвенатор читает кусочки ДНК длиной скажем 10 000 пар оснований, и допускает 30% ошибок, т.е. 30% оснований в куске будет прочитано неправильно (вместо 1 конкретного основания будет прочитано любое из 3-х оставшихся, для простоты с равной вероятностью). Вопрос - сколько раз должен быть независимо прочитан каждый кусочек (с учетом того, что он может быть гетерозиготным, т.е. что нуклеотидов в одной координате может быть 2 разных, и это нормально чтобы скажем на 99% быть уверенным в каждом прочитанном нуклеотиде этого куска?
Для любопытствующих - кратность прочтения определяется количеством индивидуальных нанопор. Сейчас их 2000, будет 8000. Т.е. это максимальная кратность прочтения в теории.

sunny82

Поясню на конкретнм примере. Есть вот такой кусок
aggtctgagc
aggcctgagc
Он гетерозиготен по 4 нуклеотиду
Был прочитан 5 раз с 30% ошибкой (примерно и результат приведен ниже
tggtcagtgc
agctctcatc
aaggctgact
gagtcttagc
cgacctgagc

Martika1

Я прикинул: получается, достаточно шести прочтений. Более того, шести прочтений будет достаточно, даже если вероятность ошибки — 40%. Однако это в предположении равных вероятностей ошибок по каждому из оснований (то есть что вероятности прочтения A таковы: A - 70%, а G, T, C - по 10%). Неравные вероятности неверных прочтений — а они существенно неравны — сильнее влияют на необходимое количество секвенирований.
В любом случае, если вероятность ошибки 30%, то 99%-ная точность достигается на не более чем 30 прочтениях. Это в случае, когда все 30% приходятся на одно и то же основание.
Добавлю, что на практике поступают иначе. Так как вероятности прочтения (base calling) скошены по парам "верное-видимое" и зависят от соседних оснований, то применяются всякие представления типа HMM и EM-алгоритмы типа Баума-Велша.

BSCurt

Тебя интересует на 99% процентов быть уверенным в правильности каждого нуклеотида куска или на 99% в правильности всего куска?

Martika1

да, уместный вопрос. Моя оценка — на 99% правильности каждого основания.
Если нужна 99%-процентная вероятность того, что весь кусок из 10000 пар оснований считан целиком верно, то для одного основания требуется вероятность ошибки, не превышающая
1 - (1 - 99%)^(1/10000) = 0.00046,
что достигается при 22 прочтениях фрагмента (в предположении вероятностей правильного:неправильных прочтений 70:10:10:10).
Да, предыдущая оценка была с ошибкой из-за бага в программе: точность 99% у каждого отдельного основания достигается при 13 прочтениях.
В худшем случае неравных вероятностей ошибок (например, A читается как A в 70% случаев и как G — в 30%) потребуется 63, 65, 67, 68, 69, 70, ... прочтений (чётное число прочтений дают худшую точность, чем следующее за ним нечётное число: лучше прочесть 63 раза, чем 66).
Это если на сей раз я правильно считаю мультиномиальные коэффициенты.

BSCurt

Я для худшего случая насчитал 31 прочтение.
Но я правда думал исключительно про гомозиготный случай
и понятно в этом случае формула для правильного прочтения нуклеотида
P=Sum_{k=0}^{(n-1)/2}C^n_k p^(n-k) q^(k)
Про гетерозиготный случай надо что-то более умное думать.

Martika1

У меня для биномиального случая, конечно, та же формула.
Но опять ошибка: если нужна 99% вероятность идеального прочтения фрагмента, то отдельное основание должно читаться правильно с вероятностью (1 - 1e-6 а не 0.9995.
Но после исправления этой ошибки получилось, что требуется ещё больше прочтений: 129
Альфа подтверждает:
sum[binomial(129,n) * 0.7^(129-n) * 0.3^n,{n, 0, 128/2}] = 0.99999899
0.99999899^10000 = 0.98997 ~ 99%

elenakozl

Одним словом, не секвенаторы, а говно!

Martika1

Возвращаясь к полиномиальному случаю.
Нужная точность считывания одного нуклеотида: 1 - 1e-6.
Вероятности неверных считываний: 0.1:0.1:0.1.
Тогда нужно 39 считываний.
Сразу предупреждаю, что где-то скопипащенная функция для подсчёта мультиномиальных коэффициентов понемногу накапливает ошибку. В частном биномиальном случае для 20 считываний отличий от Вольфрам-Альфы ещё не было, но для шестидесяти считываний ошибка уже была в пятом или шестом знаке — который нас и интересовал.
Надеюсь, что для 39 считываний в мультиномиальном случае ошибка не столь большая, и на самом деле примерно столько считываний и требуется. В общем, около сорока.
Вот код программы:

void Main
{
double p = 0.7;
double q = (1.0 - p) / 3;
//double q1 = 1.0 - p; // Для биномиального распределения
double q1 = q; // Для равных вероятностей разных ошибок
double q2 = (1.0 - p - q1) / 2;
double q3 = q2;

for (int n = 2; n < 80; n++)
{
double s = 0.0;

for (uint i = 0; i <= n / 2 + 1; ++i)
for (uint i1 = 0; i1 <= n - i; ++i1)
for (uint i2 = 0; i2 <= n - i - i1; ++i2)
{
uint i3 = (uintn - i - i1 - i2);
if (Math.Max(i1, Math.Max(i2, i3 >= i)
s += Multinomial(new uint[] {i, i1, i2, i3}) * Math.Pow(p, i) * Math.Pow(q1, i1) * Math.Pow(q2, i2) * Math.Pow(q3, i3);
}

Console.WriteLine("{0}\t{1}", n, s);
}
}


/**
* Где-то стыренный код. Мною не проверялся, но для небольших значений n в случае бинома работает верно.
*/
public static ulong Multinomial(params uint[] numbers)
{
uint numbersSum = 0;
foreach (var number in numbers)
numbersSum += number;
uint maxNumber = numbers.Max;
var denomFactorPowers = new uint[numbers.Max + 1];

foreach (var number in numbers)
for (int i = 2; i <= number; i++)
denomFactorPowers[i]++;

for (int i = 2; i < denomFactorPowers.Length; i++)
denomFactorPowers[i]--; // reduce with nominator;

uint currentFactor = 2;
uint currentPower = 1;
double result = 1;

for (uint i = maxNumber + 1; i <= numbersSum; i++)
{
uint tempDenom = 1;
while (tempDenom < result && currentFactor < denomFactorPowers.Length)
{
if (currentPower > denomFactorPowers[currentFactor])
{
currentFactor++;
currentPower = 1;
}
else
{
tempDenom *= currentFactor;
currentPower++;
}
}

result = result / tempDenom * i;
}

return (ulong)result;
}

sunny82

Ты не прав - если порядок нужного каверджа (повторности прочтения) лежит хотя бы в 100 повторных прочтениях одного фрагмента, чтобы быть на 99% уверенным в прочитанном куске, то даже 2000 пор запросто хватит для прочтения таргетного нескольких коротких и парочки длинных генов, к примеру, для верификации диагноза наследственной болячки или для проверки носительства при анамнезе и семейной истории в любой ПЦР лаборатории в регионе. Или же для прочтения де ново небольших геномов и хорошей их последующей сборки - число собранных контигов будет единичным, а это круто. Плюс такие длинные прочтения критичны для проверок транслокаций и амплификаций с точным местом встраивания, некоторых CNV и всех инверсий - другими методами это проверять зачастую проблематично.
Оставить комментарий
Имя или ник:
Комментарий: