Отвечаем на вопросы про биоинформатику
“
Нужно ли удалять (отмечать) дупликаты перед определением дифференциально экспрессированных генов, если библиотека Paired-end? На Biostars большинство говорит, что не нужно, хотя кое-где сомневаются. Если бы библиотека была SE, ответ очевиден — дупликаты удалять нельзя, т.к. многие высоко экспрессированные транскрипты сгенерируют много идентичных ридов просто за счет большого покрытия и тогда при дедупликации мы убьем их экспрессию. Но если библиотека PE и каждый фрагмент у нас прочитался 2 ридами по 100, между которыми еще и вставка 100-200? Казалось бы, в этом случае вероятность, что случайно получатся 2 одинаковых фрагмента очень маленькая и дупликаты имеет смысл удалять. Особенно, если я вижу их разное количество в ряду образцов.
Павел Мазин
Cотрудник Сколковского института науки и технологий и ИППИ РАН. Преподаватель ФББ МГУ и НИУ ВШЭ
Отвечает
Для дубликатов в данных РНК-Сек есть два возможных определения:
- Риды с полностью идентичными последовательностями.
- Риды, которые картируются в одно и то же место.
Здесь и далее и одиночные, и парные риды обозначаются словом «рид». Парные риды одинаковы, если оба соответствующих конца одинаковы.
Дубликаты могут появляться по следующим причинам:
- Высокая концентрация конкретной РНК в образце. Это не проблема, так как просто отражает уровень экспрессии генов.
- Наличие нескольких РНК с одинаковым участком (геномные дупликации). В целом тоже не проблема, так как такие риды будут картироваться в несколько мест генома и в большинстве случаев либо исключаются из анализа, либо «размазываются» по всем локациям.
- Артефакты ПЦР: при слишком большом числе циклов ПЦР фрагменты с большей склонностью копироваться постепенно вытесняют другие фрагменты. Такие артефакты могут сильно искажать оценки экспрессии и, несомненно, создают проблему.
Схлопывание всех дубликатов в один рид решит проблему с артефактами, но одновременно приведёт к уменьшению динамического диапазона определения экспрессии генов. Максимальное количество схлопнутых ридов, которые могут картироваться на транскрипт, равно его длине.
На самом деле число позиций, на которые может картироваться рид, равно L-l+1, где L и l — длины транскрипта и рида соответственно. Но мы этим пренебрегаем, так как а) обычно l<<L и б) при определении RPKM тоже следует делить не на L, а на L-l+1.
В случае парных ридов число уникальных позиций, на которые может картироваться рид, увеличивается на возможную вариабельность в длине вставки. Примем эту вариабельность за 50 нт.
Речь идёт о вариабельности в длине вставки, а не о самой длине вставки. Для наших рассуждений мы считаем вставку распределенной равномерно на интервале некоторой ширины (равной вариабельности вставки, 50 нт). В реальности размер вставки имеет примерно нормальное распределение и адекватной оценкой вариабельности будет полуширина этого распределения. Отметим, что максимальное количество схлопнутых парных ридов, которые могут картироваться на транскрипт, уменьшается с длиной вставки, так как рид должен целиком войти в транскрипт.
Проблемы с динамическим диапазоном начнутся у транскриптов, для которых риды покрывают хотя бы половину возможных позиций: у таких транскриптов не удастся получить увеличение экспрессии более, чем в два раза (обычный порог при анализе дифференциальной экспрессии). Давайте рассчитаем, когда это случится.
При случайном распределении ридов по транскрипту распределение позиций по покрытию является Пуассоновским. Нас интересует случай, когда вероятность того, что позиция не покрыта, равна 0.5. Соответственно,
poison(0, λ) = e-λ = 0.5,
где λ — вероятность рида картироваться на данную позицию, которая вычисляется так:
λ = rpkm*L*N/1000/L = rpkm*N/1000,
где rpkm — экспрессия транскрипта (ридов на тысячу нуклеотидов на миллион откартированных ридов),
L — длина транскрипта (нт),
N — размер библиотеки (млн картируемых ридов).
В случае парных ридов λ будет меньше в 50 раз (вариабельность вставки). Соответственно, пороговый RPKM для библиотеки в 20 млн ридов составит -log(0.5)*1000/N ≈ 35 для одиночных ридов и в 50 раз больше (1750) для парных ридов.
poison(0, λ) = e-λ = 0.5,
где λ — вероятность рида картироваться на данную позицию, которая вычисляется так:
λ = rpkm*L*N/1000/L = rpkm*N/1000,
где rpkm — экспрессия транскрипта (ридов на тысячу нуклеотидов на миллион откартированных ридов),
L — длина транскрипта (нт),
N — размер библиотеки (млн картируемых ридов).
В случае парных ридов λ будет меньше в 50 раз (вариабельность вставки). Соответственно, пороговый RPKM для библиотеки в 20 млн ридов составит -log(0.5)*1000/N ≈ 35 для одиночных ридов и в 50 раз больше (1750) для парных ридов.
В обычном полиА-образце тканей млекопитающих гены с такой экспрессией составляют более 25% белок-кодирующих генов для непарных библиотек и более 0.1% для парных. Для больших библиотек ситуация будет хуже. Видно, что схлопывание одиночных ридов противопоказано при хоть сколько-нибудь больших библиотеках. В случае парных ридов существенные эффекты будут проявляться только при больших размерах библиотеки (2% генов при размере библиотеки 100 млн ридов). Если библиотека не очень большая, то с некоторой осторожностью схлопывать можно.
Как описано выше, ожидаемое количество дублированных из-за экспрессии ридов в образце можно оценить, исходя из равномерного распределения ридов по транскриптам (зная длины, ожидаемые уровни экспрессии и размер библиотеки). По этим оценкам получается, что в библиотеке размером в 30 млн одиночных ридов около 50% ридов будут дубликатами, а для парных — около 4-7% (в зависимости от вариабельности вставки). Эти оценки приемлемо согласуются с реальными данными (рис 1). Таким образом, даже в парных библиотеках следует ожидать заметного числа дубликатов просто из-за высокой экспрессии отдельных генов.
Рисунок 1. Зависимость доли дубликатов (ось y) от размера библиотеки (ось х).
Красным показаны образцы с одиночными ридами, синим и зелёным — с парными ридами при вариабельности вставки в 50 и 100 нт соответственно.
Теоретические оценки показаны линиями, результаты для реальных данных — точками.
Красным показаны образцы с одиночными ридами, синим и зелёным — с парными ридами при вариабельности вставки в 50 и 100 нт соответственно.
Теоретические оценки показаны линиями, результаты для реальных данных — точками.
Артефакты ПЦР приводят не просто к перепредставленности отдельных РНК: перепредставлены оказываются отдельные фрагменты транскриптов. На практике это выглядит как крайне неровное покрытие транскрипта (экзона) ридами, с характерными «стопками» ридов.
Рисунок 2. Первый ряд — нормальное покрытие, второй — покрытие курильщика (со «стопками»), третий — структура генома на данном участке.
Оценить такую неровность можно по отношению дисперсии к среднему у распределения числа ридов, которые начинаются с каждой из позиций транскрипта (экзона). В идеале распределение должно быть Пуассоновским и отношение должно быть равно 1. В хороших данных для большинства экзонов это отношение оказывается не более, чем 2-3.
Рисунок 3. Распределение внутренних экзонов по отношению дисперсии к среднему (в log шкале) для четырех датасетов. Зеленый и синий в целом нормальные, красный и фиолетовый — не очень.
Обнаружение «стопок» мало поможет при анализе текущих данных — проблемы со схлопыванием остаются в силе. Зато вы сможете дать рекомендации для последующих экспериментов — например, уменьшить число циклов ПЦР.
Итак, схлопывание дубликатов имеет смысл использовать при наличии «стопок» в парных библиотеках не очень большого размера. В остальных случаях схлопывание противопоказано, особенно для одиночных ридов.
R-код для оценки доли дубликатов
calcDuplicateFreqInTransc = function(r,L,N,insert=1){
l = r*N/1000/insert
z = r*N*L/1000
c(total=z,dupl=z - (1-exp(-l))*L*insert)
}
#' @param rs vector of gene RPKMs
#' @param ls vector of gene lengths (longest transcript was used)
#' @param N lib size (in millions)
#' @param insert insert variability
calcDuplicateFreqInLib = function(rs,ls,N,insert=1){
rs = rs / sum(rs*ls/1e9) #RPKM was estimated using different gene models (only constant exons were considered) so, RPKM need to be adjusted to fit library size
z=apply(sapply(1:length(rs),function(i)calcDuplicateFreqInTransc(rs[i],ls[i],N,insert)),1,sum)
z[2]=z[2]/z[1]
z
}
#mean.exp = real human RPKM of protein coding genes
#gene.len = lengths of longest transcripts of the same genes as in mean.exp
calcDuplicateFreqInLib(mean.exp,gene.len,20,1)
Наши контакты
Телефон: +7 916 088 13 07
E-mail: hello@ksivalue.com
ООО «Ксивелью»
ИНН 7702424959, ОГРН 5177746030831
Почтовый и фактический адрес: 119049, Москва, Ленинский проспект, 30А. Схема проезда
E-mail: hello@ksivalue.com
ООО «Ксивелью»
ИНН 7702424959, ОГРН 5177746030831
Почтовый и фактический адрес: 119049, Москва, Ленинский проспект, 30А. Схема проезда