KR102408080B1 - 전사체 결정 방법 - Google Patents

Abstract

본 발명은 다음 단계를 포함하는, 전사체 존재비를 추정하는 방법에 관계한다:
a) 관심되는 유전자 좌위의 전사체의 잠재적 혼합체로부터 전사체 단편 염기서열결정 데이터를 획득하는 단계,
b) 상기 단편 염기서열결정 데이터를 관심되는 상기 좌위의 유전자 좌표에 배정하여, 단편 유전자 좌표 커버리지의 데이터 세트를 획득하되, 각 유전자 좌표에 대한 상기 커버리지가 통합되어 커버리지 포락선을 형성하는 단계,
c) 상기 혼합체의 전사체의 숫자를 세팅하는 단계,
d) 각 전사체 i에 대한 모형화된 유전자 커버리지의 확률 분포 함수를 미리 세팅하되, i가 전사체에 대한 수치 식별자를 나타내고, 상기 확률 분포 함수가 상기 전사체 i의 가중 인자 α_i의 수학적 결과물 및 최소한 2개의 확률 부함수 j의 합계로 구성되고, j가 확률 부함수에 대한 수치 식별자를 나타내고, 각 확률 부함수 j가 가중 인자 β_i,j에 의해 독립적으로 가중되는 단계,
e) 각 전사체에 대한 확률 분포 함수를 더하여 합계 함수를 획득하는 단계,
f) 합계 함수를 커버리지 포락선에 적합시키고, 따라서 α_i와 β_i,j에 대한 값을 최적화하여 적합을 증가시키는 단계,
g) 미리 세팅된 수렴 규준 실현될 때까지 단계 e)와 f)를 반복하고, 따라서 수렴 규준이 실현된 후에 최적화될 때 가중 인자 α_i에 의해 제공된 혼합체의 각 전사체에 대한 추정된 전사체 존재비를 획득하는 단계.

Description

전사체 결정 방법{TRANSCRIPT DETERMINATION METHOD}

발명의 분야

본 발명은 차세대 염기서열결정 (NGS) 리드 (read)에 기초하여 전사체 (예를 들면, mRNA) 존재비의 정보를 제공하는 것에 관계한다.

배경

차세대 염기서열결정 기술은 핵산 표본을 염기서열결정할 때, 많은 양의 짧은 리드를 생산한다. 차세대 염기서열결정에서 필수적인 단계는 라이브러리 제조 또는 요약하면 라이브러리 프렙 (prep)이다. 이러한 과정은 mRNA 또는 cDNA를 입력으로서 취하고, 그리고 mRNA 분자의 섹션에 각각 상응하는 짧은 cDNA 단편의 라이브러리를 생산한다. 이들 단편은 이후, 통상적으로 전체가 아닌 그들의 시작에서 및/또는 그들의 끝에서 부분적으로 NGS 서열분석기에 의해 염기서열화된다. 이것은 리드로 불리고, 그리고 NGS 서열분석기에 의해, 유전자 코드의 핵염기를 나타내는 4개의 ASCII 문자, 예를 들면, A, C, G, T 또는 0, 1, 2, 3의 군의 서열로서 가장 흔히 저장되는 뉴클레오티드의 짧은 서열을 유발한다. 어떤 mRNA 분자가 본래 표본에서 존재하였는지를 추론하기 위해, 이들 리드는 참고 유전체 위에 매핑된다.

차세대 염기서열결정은 예로서, 매핑된 유전체를 이용하여 서열 리드를 일정한 생물체 변이체 (WO 2009/085412 A1)에 연관시킴으로써, 다양한 유전체 매핑 절차 (US 2013/110410 A1) 또는 DNA 확인 방법에서 이용되었다.

WO 2009/091798 A1은 생물체의 전사체의 프로필을 획득하기 위한 방법을 설명하고, 상기 방법은 다음을 포함한다: 하나 또는 그 이상의 cDNA 분자를 염기서열결정하여 서열 리드를 획득하고; 각 서열 리드를 참고 서열에 맞춰 정렬한다.

하지만, 이전 방법에서 알려지지 않은, 짧은 서열 리드를 이용한 전사체 분석의 근원적인 주요 문제점은 전사체 변이체, 예를 들면, 서열 편차에서 다른 동종형, 예컨대, 예를 들면 유전자의 계통 차이, 점 돌연변이, 또는 한 단백질의 스플라이스 변이체의 경우에 정렬 단계이다. 짧은 서열 리드를 하나의 전사체 변이체에 맞춰 정확하게 정렬하는 것은 통상적으로 어렵다.

서열 리드에 기초하여 전사체 염기서열결정 데이터를 집합시키는 가장 흔한 방법은 "Cufflinks" 방법 (Trapnell et al. 2010)이다. Cufflinks는 RNA-Seq 실험에서 관찰된 리드를 "설명하는" 전사체의 축약형 세트를 작제한다. 이것은 비교 어셈블리 문제를, 이분 그래프에서 최고 정합에서 문제로 환원시킴으로써 그러한다. 본질적으로, Cufflinks는 리드 정렬에서 커버링 관계를 구축하고, 그리고 이러한 관계에 대한 비순환 방향 그래프에서 최소 경로 커버를 발견함으로써, 딜워스의 정리의 건설적인 증명을 실행한다. 이러한 통계학적 방법을 이용하여, Cufflinks는 공지된 참고 주해를 이용하여, 또는 단지 참고 유전체만을 이용하여 전사체의 순이론적 어셈블리 후, 표본 내에 존재하는 전사체 동종형의 존재비를 추정할 수 있다. Cufflinks는 한 세트의 단편을 고려할 때 한 세트의 전사체의 존재비에 대한 우도를 도출하기 위해, 대합된-단부 염기서열결정 실험의 통계 모형을 이용한다. 이러한 우도 함수는 독특한 최대값을 갖는 것으로 보일 수 있는데, 이것은 수치 최적화 알고리즘을 이용하여 Cufflinks에 의해 발견된다. 상기 프로그램은 이후, 이들 확률을 곱하여, 전사체에서 제안된 존재비를 고려할 때 실험에서 단편이 관찰될 전체 우도를 연산한다. Cufflinks의 통계 모형이 선형이기 때문에, 우도 함수는 독특한 최대값을 갖는데, 이것은 수치 최적화 알고리즘으로 Cufflinks에 의해 발견된다.

Roberts et al.(2011)은 단편 편향을 교정함으로써, RNA-Seq 발현 추정값을 향상시키는 방법에 관계한다. Wen-Ping et al.(2007)은 자연 개체군에서 전사체 존재비 부류의 혼합체 모형화를 설명한다.

이전 방법은 전사체 변이체를 정확하게 구별하고, 그리고 다른 전사체에 관하여 정확한 전사체 양 또는 존재비를 획득하는데 실패하였다. 본원에서 비교에 의해 보여진 바와 같이, Cufflinks 방법조차도 여러 실험에서, 정확한 전사체 존재비 정보에 도달하는데 실패하였다.

본 발명의 목적은 전사체 존재비의 더욱 정확한 사정을 허용하는 향상된 방법을 제공하는 것이다.

발명의 요약

본 발명은 다음 단계를 포함하는, 전사체 존재비를 추정하는 방법을 제공한다: a) 관심되는 유전자 좌위 내에서 전사체의 잠재적 혼합체로부터 전사체 단편 염기서열결정 데이터를 획득하는 단계,

b) 상기 단편 염기서열결정 데이터를 관심되는 상기 좌위의 유전자 좌표에 배정하여, 단편 유전자 좌표 커버리지의 데이터 세트를 획득하되, 각 유전자 좌표에 대한 상기 커버리지가 통합되어 커버리지 포락선 (또한, 전체 커버리지 히스토그램 또는 전체의 히스토그램으로서 지칭됨)을 형성하는 단계,

c) 상기 혼합체의 전사체의 숫자를 세팅하는 단계,

d) 각 전사체 i에 대한 모형화된 유전자 커버리지의 확률 분포 함수를 미리 세팅하되, i가 전사체에 대한 수치 식별자를 나타내고, 상기 확률 분포 함수가 상기 전사체 i의 가중 인자 α_i의 수학적 결과물 및 최소한 2개의 확률 부함수 j의 합계로 구성되고, j가 확률 부함수에 대한 수치 식별자를 나타내고, 각 확률 부함수 j가 가중 인자 β_i,j에 의해 독립적으로 가중되는 단계,

e) 각 전사체에 대한 확률 분포 함수를 더하여 합계 함수를 획득하는 단계,

f) 합계 함수를 커버리지 포락선에 적합시키고, 따라서 α_i와 β_i,j에 대한 값을 최적화하여 적합을 증가시키는 단계,

g) 미리 세팅된 수렴 규준이 실현될 때까지 단계 e)와 f)를 반복하고, 따라서 수렴 규준이 실현된 후에 최적화될 때 가중 인자 α_i에 의해 제공된 혼합체의 각 전사체에 대한 추정된 전사체 존재비를 획득하는 단계.

본 발명은 예로서, 컴퓨터에서 상기 방법과 단계를 수행하거나 또는 보조하기 위한 기계 코드를 내포하는, 상기 방법을 이용하는 컴퓨터 프로그램 제품을 더욱 제공한다. 컴퓨터 프로그램 제품은 임의의 종류의 기억 장치에 담겨 제공될 수 있다. 본 발명의 방법의 단계를 수행하는데 보조하도록 프로그램된 시스템, 예를 들면, 컴퓨터 장치 역시 제공된다. 계산 단계는 통상적으로, 오퍼레이터의 도움 없이 수행된다. 입력과 세팅 단계는 예로서, 단계 d)에서 확률 부함수의 숫자와 유형에 대한 옵션 제안을 제안함으로써 프로그램 또는 시스템에 의해 보조될 수 있다. 당연히, 프로그램 또는 시스템은 또한, 오퍼레이터로부터 추가 입력 없이 디폴트 파라미터에서 수행될 수 있다.

다음의 상세한 설명과 바람직한 구체예는 본 발명의 모든 양상에 적용되고, 그리고 명시적으로 지시된 경우를 제외하고, 제한 없이 서로 결합될 수 있다. 바람직한 구체예와 양상은 청구항에서 규정된다.

발명의 상세한 설명

본 발명은 전사체 단편 서열의 표본으로부터 전사체 존재비 정보를 획득하기 위해 수치 방법을 이용한다.

상기 방법은 일반적으로 전사체 단편 서열로서 지칭되는 리드를 참고 서열, 예를 들면, 참고 유전체에 맞춰 정렬하여 (NGS), 유전자 커버리지 정보를 획득한다 (단계 b). 이런 목적으로 이용된 이전 통계학적 도구는 종종, 관찰된 데이터의 성격에 관해 비현실적인 가정을 하고, 이런 이유로 전사체 농도의 추정값이 부정확하다. 가장 폭넓게 이용되는 도구 중에서 일부, 예를 들면, Cufflinks는 전사체를 따라서 리드의 분포가 균일하다고 가정하는데, 이것은 현재 mRNA-Seq 프로토콜과 상충된다. 본 발명은 전사체를 따라서 리드 분포의 편향 및 전사체 존재비를 동시에 학습할 수 있는 통계 모형을 제공한다. 이런 목적으로, 전사체의 리드 또는 단편 분포는 적합 단계에서 전사체 존재비와 함께 훈련되는 함수의 혼합체에 의해 모형화된다. 적합 단계에서 이용된 방법은 기대값 최대화 알고리즘으로 이전 최대화 또는 최소화 절차, 예를 들면, 최고 우도 프레임워크로부터 추론될 수 있다. 본 발명의 모형에서 리드의 전체 확률 분포가 혼합체의 혼합체이기 때문에, 이러한 모형은 Mix² (리드: Mixquare) 모형으로 불린다. 아래에서, Mix² 모형은 매우 다능하고, 그리고 임의선택적 파라미터 묶기에 의해 데이터에 내재하는 상이한 구조에 맞게 조정될 수 있는 것으로 나타난다. 특히, 전사체 존재비를 획득하는데 이용된 이들 방법은 전사체 관련된 확률 분포에 적합할 수 있다. 실험에서, Mix² 모형은 Cufflinks 프로그램에서 이용된 통계 모형보다 전사체 존재비에 대한 훨씬 우수한 추정값을 달성하는 것으로 나타난다. 심지어 부정확한 전사체 주해로부터 시작할 때에도, Mix² 모형은 데이터로부터 정확한 주해를 학습하고, 그리고 선행 기술의 것들보다 훨씬 뛰어난 존재비 추정값을 산출할 수 있다. 적합 단계 동안 우수한 학습 능력으로 인해, 예로서 배정 단계 a) 동안 또는 단계 d)에서 (예를 들면, 무작위) 확률 분포 함수의 선택 동안 선택된 시작 파라미터는 필수적이지 않다. 심지어 전사체의 가정된 숫자가 상이할 수 있다. 잘못된 전사체 주해 또는 전사체 숫자 가정은 예로서, 하나 또는 그 이상의 전사체에 대한 확률 분포 함수를 제로의 존재비에 적합시킴으로써 교정될 것이다. 확률로서 모형화될 수 있는 가중 인자 알파는 수렴 후 전사체의 존재비를 나타낼 것이다.

본원에서 이용된 바와 같이, "관심되는 유전자 좌위"는 상기 방법을 염색체 상에서 유전자의 하나의 연속적 서열 스트레치에 한정하지 않는다. 이것은 일반적으로, 유전자 서열의 하나 또는 그 이상의 섹션을 지칭한다. 유전자 좌위는 유전체 좌표와 연관되고, 서열 리드 (전사체 단편 서열)에 대한 위치 정보를 제공할 수 있다. 유전자 "위치" 또는 "좌표"는 참고 서열의 시작으로부터 일정한 거리를 갖는 참고 서열 상에서 수치적으로 확인된 뉴클레오티드를 지칭하기 위해 본원에서 이용된다. 유전자 좌표는 유전자 좌표가 전사체와 양립하면 유전체 좌표 또는 전사체 좌표에서 표현될 수 있는데, 이것은 유전자 좌표가 전사체 엑손 중에서 한 가지 내에 위치되면 그러하다. 양립성 유전체 좌표는 전사체 시작으로부터 상대적 거리를 계산하고 유전체 좌표를 선행하는 인트론의 길이를 감산함으로써 전사체 좌표로 전환된다. 전사체 좌표는 이러한 과정을 되돌림으로써 유전체 좌표로 전환될 수 있다. 한 전사체의 좌표 상에서 규정된 확률 분포는 확률 부함수의 수축, 신장과 이동에 의해, 다른 전사체의 전사체 좌표 상에서 확률 분포로 전환될 수 있다.

핵염기 유형 (예를 들면, A, T/U, G, C)은 유전체 좌표 또는 위치와 연관되거나 또는 연관되지 않을 수도 있다. 통상적으로, 한 가지 유형의 핵염기가 각 전사체에 대한 각 유전자 좌표와 연관된다. 하지만, 각 유전자 좌표에 대해 상이한 중복 전사체는 핵염기 구성에서 상이할 수 있는데, 그 이유는 본 발명이 점 돌연변이를 식별하는데 이용될 수 있기 때문이다, 다시 말하면, 특히 표본이 상이한 생물체의 또는 상이한 대립형질의 핵산 분자를 내포하면, 하나 또는 그 이상의 핵염기가 상이한 전사체 사이에서 상이할 수 있기 때문이다. 심지어 이들 전사체, 이들의 단편 및 관심되는 좌위 사이에 가능한 부정합의 경우에도, 통상적인 유전자 좌표가 할당될 수 있다. 단계 b), 단편 서열을 유전자 좌표에 배정하는 것은 이들 서열 사이에 서열 비교 또는 정렬을 포함할 수 있는데, 이것은 유전자 좌표를 제공한다. 서열 비교는 당분야에서 널리 공지되고 (예를 들면, 언급된 개방 소스 cufflinks 방법에 의해) 뉴클레오티드의 비교를 포함한다. 본 발명의 방법은 관심되는 소정의 유전자 좌위 상에서 일탈하는 서열을 갖는 상이한 균주 또는 생물체의 데이터에서 이용될 수 있거나, 또는 상이한 스플라이스 변이체, 다시 말하면, 엑손 서열의 상이한 조합에 의해 식별될 수 있는 상이한 전사체를 구분하는데 이용될 수 있다. 따라서 배정 단계 a) 동안 부정합은 허용되거나 또는 불허될 수 있다. 이런 부정합은 바람직하게는, 100개 염기당 기껏해야 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 또는 11개 부정합이다.

본 발명은 전사체의 혼합체에서 개별 전사체의 존재비를 모형화할 수 있다 (예를 들면, 확률에 상응하는 상대적 양으로서). 전사체 단편 염기서열결정 데이터는 최소한 2개, 바람직하게는 3, 4, 5, 6, 7, 8, 9, 10개, 또는 최소한 또는 기껏해야 15개, 최소한 또는 기껏해야 20개, 최소한 25개, 최소한 또는 기껏해야 30개, 최소한 또는 기껏해야 40개의 전사체 서열을 포함할 수 있다. 전사체의 숫자는 예로서, 핵산을 증폭하거나 또는 제거함으로써, 예로서 핵산 표본에서 수행된 선택 단계에서 선택될 수 있다. 이런 제거된 또는 증폭된 핵산은 통상적인 서열 스트레치, 예를 들면, 선택되는 올리고뉴클레오티드 프로브, 앵커 서열 또는 프라이머 서열에 상응하는 서열을 내포할 수 있다. 하나 또는 그 이상의 유전체 좌위가 이러한 방식으로 선택될 수 있다. 본 발명의 특정 바람직한 구체예에서, 한 유전자의 전사체 변이체가 조사된다. 그럼에도 불구하고 일반적으로, "전사체"는 본원에서 이용된 바와 같이, 임의의 유전자 또는 유전자 조합의 임의의 핵산 또는 이의 서열 및 이의 임의의 변이체, 특히 이의 mRNA 또는 cDNA 서열 변이체를 지칭할 수 있다.

바람직하게는, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 5 내지 800개 뉴클레오티드, 바람직하게는 6 내지 600개 뉴클레오티드, 7 내지 400개 뉴클레오티드, 8 내지 200개 뉴클레오티드 또는 9 내지 150개 뉴클레오티드, 훨씬 바람직하게는 10 내지 100개 뉴클레오티드, 특히 바람직하게는 12 내지 70개 뉴클레오티드의 길이를 갖는다.

전사체 단편 서열의 숫자는 각 전사체에 대해 최소한 또는 기껏해야 100개, 최소한 또는 기껏해야 500개, 최소한 또는 기껏해야 1000개, 최소한 또는 기껏해야 5000개, 최소한 또는 기껏해야 10000개, 바람직하게는 최소한 또는 기껏해야 10개, 또는 최소한 또는 기껏해야 20개, 최소한 또는 기껏해야 50개일 수 있다. 조합으로 또는 대안으로서, 전사체 단편 서열의 숫자는 기껏해야 400000, 기껏해야 300000, 기껏해야 200000, 기껏해야 100000 또는 기껏해야 50000일 수 있다.

최소한 하나 또는 그 이상의, 예를 들면, 모든 전사체의 전사체 길이는 예로서, 100 내지 1000000개 뉴클레오티드, 바람직하게는 1000 내지 100000개 뉴클레오티드 또는 2000 내지 10000개 뉴클레오티드일 수 있다.

바람직한 구체예에서, 관심되는 유전자 좌위는 예로서, 하나 또는 그 이상의 유전자 또는 유전자 원소의 전사체 서열을 인코딩하는 하나 또는 그 이상의 동종형을 포함하고, 바람직하게는 한 유전자 또는 유전자 원소의 최소한 2, 3, 4개 또는 그 이상의 스플라이스 변이체를 포함한다. 이것은 다른 유전자 또는 유전자 원소의 하나 또는 그 이상의 다른 스플라이스 변이체를 포함할 수 있다. 스플라이스 변이체에 더하여 또는 이의 대안으로, 이것은 상이한 대립형질을 포함할 수 있다. 바람직한 구체예에서, 유전자 또는 유전자 원소는 단백질 (예를 들면, mRNA)을 인코딩하고, 또한 비단백질-코딩 전사체, 예를 들면, 마이크로RNA, snoRNA 또는 rRNA를 비롯한 조절성 또는 촉매성 RNA뿐만 아니라 그들의 전구체, 특히 프리-마이크로RNA 또는 프리-rRNA를 포함한다.

본원에서 이용된 바와 같이, "유전자" 및 "유전자 원소"는 전사되면 하나 또는 그 이상의 전사체를 형성하는 서열을 갖는 유전자 뉴클레오티드에 관계한다.

본원에서 이용된 바와 같이, "동종형"은 전사체의 특정 변이체에 관계하는데 이용된다. 한 "유전자" 또는 "유전자 원소"의 전사체는 예로서, 스플라이스 변이체의 경우에 상이하고, 따라서 상이한 동종형을 발생시킬 수 있다. 다른 동종형 변이는 예로서, 상이한 미생물 또는 균주의 혼합체에서 전사체 물질의 상이한 대립형질 또는 상이한 공급원에 의해 유발될 수 있다.

전사체의 숫자를 세팅하는 단계는 관심되는 유전자 좌위로부터 미리 주해된 서열 데이터를 획득하고, 그리고 전사체의 숫자를 최소한, 상기 관심되는 유전자 좌위로부터 예상되는, 상이한 전사체 서열로서 계수되는 스플라이스 변이체를 포함하는 상이한 전사체 서열의 숫자로 세팅하는 것을 포함할 수 있다. 앞서 말한 바와 같이, 전사체의 세트 숫자는 전사체의 실제 숫자 (이것은 정확하게 공지되거나 또는 공지되지 않을 수도 있다)를 초과할 수 있는데, 그 이유는 부정확한 전사체가 합계 함수의 적합 동안 제거되고, 하나 또는 그 이상의 가중 인자 알파가 제로로 수렴하도록 유발될 수 있기 때문이다. 통상적으로, 관심되는 유전자 좌위는 예로서, 핵산 선택 단계 중에서 임의의 한 가지의 cDNA 산출 단계로부터 알려진다. 주해된 유전자 데이터를 이용하여, 전사체의 양, i_max에 대한 시작 숫자에 도달하는 것이 가능하다. 각 전사체에 대해, 존재비는 단계 d)에서 확률 분포 함수를 세팅함으로써 모형화될 것이다. 이러한 함수는 가중 인자 알파를 내포하는데, 이것은 반복적인 과정인 합계와 확률 분포 적합의 수렴 후 혼합체에서 전사체의 존재비에 상응한다. 전사체의 혼합체는 본 발명의 모형의 첫 번째 혼합체이다. 각 알파는 적합 과정 동안 별개로 변형될 수 있다. 따라서, 전사체의 숫자를 세팅하는 것은 대략 추정값일 수 있다. 예를 들면, 이것은 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 20, 30, 40 또는 그 이상 중에서 최소한 한 가지일 수 있다.

두 번째 혼합체는 각 확률 분포 함수에서 수학적으로 내포된다. 확률 분포 함수는 2개 또는 그 이상 (예를 들면, 2, 3, 4, 5, 6, 7, 8, 9, 10개 또는 그 이상)의 확률 부함수 또는 "블록"으로 구성된다. 전사체의 길이 (전사체 좌표 범위)에 걸쳐 비대칭 리드 분포를 모형화할 수 있기 위해, 최소한 2개의 확률 부함수를 갖는 것이 필수적이다. 전사체당 확률 부함수의 총량은 j_max이다. J_max는 통상적으로 대략 4-8이다. 너무 많은 양의 확률 부함수는 과훈련 및 감소된 연산 효율을 야기할 수 있다.

확률 분포 함수는 확률 부함수의 합계이다. 주목할 점은 인자 알파가 전체 확률 분포 함수, 그리고 따라서, 각 기본적 확률 부함수를 동등하게 가중한다는 것이다. 확률 분포 함수 내에서, 확률 부함수는 가중 인자 β_i,j. (베타_i,j)에 의해 별개로 가중된다. 중요하게는, 각 전사체, 다시 말하면, 각 확률 분포 함수에 대해, 베타는 독립 또는 종속 인자 ("묶인")일 수 있다. 독립성의 경우에, 각 베타는 적합 단계 동안 개별적으로 변형될 수 있다. 상이한 확률 분포 함수 사이에 베타, 예를 들면, 첫 번째 베타를 묶는 것이 또한 가능하거나, 또는 일반적으로 임의의 소정의 베타_j가 각 확률 분포 함수에 대해 동일할 것이다 (예를 들면, 베타_1,j = 베타_2,j). 이것은 모형화된 가중 인자 베타의 숫자를 i_max x j_max 대신에 단지 j_max로 감소시키는데, 이것은 적합 과정을 단순화하고 연산적 자원 이용을 감소시킨다. 이러한 단순화는 각 전사체에 대해 유사한 단편 커버리지 분포가 예상될 수 있으면, 예를 들면, 다른 부분에서보다 (2개 또는 그 이상, 예를 들면, 전체 전사체 ("묶인" 전사체)에 대해) 소정의 서열 부분에서 항상 더욱 큰 존재비가 예상되면, 최고로 작동한다. 이러한 가정은 통상적으로, 대략 동일한 또는 유사한 길이의 전사체, 예를 들면, 다른 묶인 전사체의 0.3 내지 3.2, 바람직하게는 0.5 내지 2.1의 길이 (핵염기에서)를 갖는 전사체에 대해 사실이다. 당연히, 아래에서 더욱 명확해지는 바와 같이, 확률 부함수의 다른 파라미터 역시 상이한 전사체 사이에, 또는 또한, 소정의 전사체의 확률 부함수 사이에 유사한 방식으로 묶일 수 있다 (예를 들면, 전사체에 대한 전체 확률 분포 함수를 유전자 좌표 방향에서 횡단적으로 이동할 때). 또한, 단계 d)에서 전사체의 확률 부함수를 서로에 대해 대략 동일한 최고 높이 또는 최대값 (통상적으로, 확률 분포 함수의 부분 내지 전체의 확률 또는 확률 기여로서 모형화됨)에 정규화하는 것이 가능하다.

파라미터의 묶기를 위한 추가 옵션에 따라, 단지 하나 또는 그 이상의 파라미터, 예를 들면, 하나 또는 그 이상의 베타를 일군의 전사체, 다시 말하면, 각 확률 분포 함수 내에서 묶는 것이 또한 가능하다 - 하지만 모든 전사체, 다시 말하면, 각 확률 분포 함수를 통해서는 그렇지 않다. 이런 군은 확률 분포 함수의 합계 함수의 예상된 유사한 모양에 의해 규정될 수 있다. 예를 들면, 군은 유사한 길이 및/또는 유사한 GC-함량에 의해 규정될 수 있다. 유사한 길이는 예로서, 군의 모든 구성원이 모든 구성원의 평균 크기의 +/- 50%, 바람직하게는 +/- 35%, 더욱 바람직하게는 +/- 30% 이내에 길이를 가질 때이다. 유사한 GC-함량은 예로서, 모든 구성원의 평균 GC-함량의 +/- 10%, 바람직하게는 +/- 5% 이내일 수 있다.

"약"은 본원에서 이용된 바와 같이, 동일한 값 또는 소정의 값의 +/- 10% 차이나는 값을 지칭할 수 있다.

전사체 i에 대한 확률 부함수의 합계는 소정의 가중 인자와 함께, 전사체 i에 대한 확률 분포 함수를 구성한다. 확률 부함수로서 임의의 수학적 함수는 이들 확률 부함수의 합계가 확률 함수를 형성할 수 있다는 요건에서 선택될 수 있다. 확률 부함수가 적합되고 최적화될 수 있는 연산 모형에 대한 기초를 형성하기 때문에, 소수의 파라미터를 갖는 단순한 함수가 선호된다. 확률 부함수는 예로서, 유전자 좌표의 의존에서 함수의 모양을 결정하는 1, 2, 3, 4, 5, 6, 7, 8, 9 또는 10개의 함수 파라미터를 포함하거나 또는 이들로 구성된다. 확률 부함수 j는 바람직하게는, 각 유전자 좌표에 대해 양성 값으로 구성된다. 바람직하게는, 확률 부함수는 비주기 함수이고 및/또는 특히 바람직하게는, 확률 부함수는 밀도 함수 또는 확률 함수이지만, 이것을 확률 부함수의 합계 (이것은 확률 분포 함수 (소정의 전사체에 대한)로서 지칭된다)로부터 식별하기 위해, 본원에서 확률 부함수 또는 "블록"으로서 지칭된다. 이것은 통상적으로, 일정한 유전자 좌표에서 최대값을 내포하고, 그리고 상기 최대값으로부터 떨어진 양성과 음성 유전자 좌표의 경우에 꾸준히 감소한다. 바람직하게는, 확률 부함수는 단지 하나의 최대값을 갖는다. 실례 확률 부함수 j는 가우스 함수, 사각형 함수, 삼각형 함수, 바람직하게는 가우스 함수에서 선택된다.

유전자 좌표 또는 이의 범위는 바람직하게는, 각 전사체에 대해 동일하다. 이를 위해 전사체에 대한 확률 분포 함수는 유전자 좌표 방향에서 수축되거나 또는 신장되거나, 또는 이동될 수 있다. 당연히, 이런 변환은 참고 서열, 예를 들면, 유전체에서 정확한 공급원 좌표에 다시 한 번 도달하기 위해 역전될 수 있다. 하지만 이런 변환은 존재비 모형화 알고리즘에 대해 유익할 수 있다.

이러한 수축, 신장과 이동은 정확한 (또는 가능성이 더욱 높은) 전사체 길이에 관한 지식을 획득하거나 또는 향상시키는데 이용될 수 있다. 등용에서, 전사체 길이 추정값이 이러한 단계에 의해 획득되는데, 이것은 단편의 숫자 (예를 들면, 양 또는 농도)를 결정하는데 이용된다. 단편 숫자는 상관관계: FPKM = 단편의 숫자/전사체 길이 (x 가중 인자)에 기초하여, FPKM (백만 리드당 킬로 염기쌍마다 단편) 값으로부터 추정될 수 있다. 따라서, 단편의 숫자는 FPKM x 전사체 길이와 상관한다.

전사체 길이 추정값은 단계 e)의 수렴 후 합계 함수의 임의선택적으로 수축된 또는 신장된 및/또는 이동된 확률 분포 함수에 따라 전사체 시작과 끝 위치를 세팅함으로써 획득될 수 있다. 규정된 시작과 끝 위치를 가질 수 없고, 단지 이에 대한 확률 값만을 갖는 확률 분포 함수로 구성되는 합계 함수의 성격으로 인해, 시작과 끝은 예로서, 상기 유전자 좌표까지 (시작의 경우에) 또는 상기 유전자 좌표로부터 시작하여 (끝의 경우에) 부분적인 구역이 곡선 아래 전체 면적 (합계 함수 또는 첫 번째 또는 마지막 확률 분포 함수의)의 분율인 유전자 좌표로서 추정될 수 있다. 이러한 분율은 1% 내지 10% 사이에 값일 수 있다. 숙련자는 이용된 확률 분포 함수의 모양에 따라 이런 구역 컷오프를 결정하기 위해 적합한 값을 쉽게 시험할 수 있다. 이것은 공지된 시작과 끝 위치를 갖는 모형 핵산을 이용하고, 다른 시작과 끝 위치를 갖는 본 발명의 함수를 설정하고, 그리고 본 발명의 알고리즘이 설명된 바와 같이 향상된 시작과 끝 위치를 훈련하게 함으로써 쉽게 시험될 수 있다. 컷오프 값에 대해 이용된 실례 값은 0.5%, 1%, 2%, 3%, 4%, 5%, 6%, 7%, 8%, 9% 또는 10%이다.

추가 옵션에 따라, 단편 길이 분포 또는 평균 단편 길이는 단계 e)의 수렴 후 본 발명의 모형 함수를 이용하여 계산된다.

특히, 하나 또는 그 이상, 바람직하게는 각 전사체 (및 이의 확률 분포 함수)에 대한 유전자 좌표는 관심되지 않는 유전자 영역이 제외되도록 임의선택적으로 변환된, 유전체 내에 뉴클레오티드 위치에 상응할 수 있는데, 여기서 바람직하게는, 관심되지 않는 상기 유전자 영역은 상기 전사체 단편 염기서열결정 데이터에 의한 커버리지를 내포하지 않는다 - 예를 들면, 어떤 리드 또는 서열 단편도 정렬하지 않는 인트론. 따라서 연속적 확률 분포 함수는 간헐적 제로 분포 없이 모형화될 수 있다. 관심되지 않는 이런 영역은 개별 확률 분포 함수로부터, 그리고 부가적으로 또는 대안으로 또한, 커버리지 포락선 (이에 대해, 다양한 확률 분포 함수의 합계 함수가 각 반복 단계 g) 동안 적합된다)으로부터 제거될 수 있다. 포락선은 통합된 서열 데이터 커버리지의 히스토그램을 나타낸다.

따라서 바람직하게는, 본 발명의 방법은 인트론을 갖는 유전자 좌표 위치를 상기 커버리지 포락선으로부터 제거하는 것을 포함하는 단계 b2)를 더욱 포함한다. 이것은 변형된 커버리지 포락선을 획득하는 것을 유발하는데, 이것은 단계 f) 및 바람직하게는, 모든 다른 단계에서 - 당연히, 확률 분포 함수가 본래 단편 염기서열결정 데이터 또는 관심되는 좌위의 유전자 또는 유전체 좌표에 다시 재조회되는 경우를 제외하고, 커버리지 포락선의 이용을 대체한다.

유전자 좌표 커버리지는 존재비를 모형화하기 위해, 리드의 전체 서열 정보를 이용하는 것이 필요하지 않다. 또한, 리드 서열의 단지 제한된 숫자의 뉴클레오티드, 예를 들면, 단지 첫 번째 하나의 또는 시작 부위 정보만을 이용하는 것으로 충분하다. 단편 유전자 좌표 커버리지는 유전자 좌표에 배정된 각 단편 서열에 대한 최소한 하나의 뉴클레오티드의 수치를 내포할 수 있고, 바람직하게는 여기서 최소한 하나의 뉴클레오티드는 단편 시작 부위 또는 전체 단편 서열을 포함한다. "포함한다"는 본원에서 이용된 바와 같이, 추가 구성원이 내포되는 것을 허용하는 열린 정의로서 이해될 것이다. "구성되는"은 다른 한편, 구성되는 정의 특질의 추가 원소 없이 닫힌 정의로서 고려된다. 따라서 "포함하는"은 더욱 넓은 정의이고 "구성되는" 정의를 내포한다. 본원에서 "포함하는" 언어를 이용하는 임의의 정의는 또한, 본 발명의 특수한 구체예에서 구성 제한으로 읽혀질 수도 있다.

확률 분포 함수의 성분인 2개 또는 그 이상의 확률 부함수는 바람직하게는, 상기 확률 분포 함수 (전사체의)의 유전체 좌표 범위 사이에서 균등하게 분포된다. 일반적으로, 바람직하게는 전사체에 대한 확률 부함수는 각각 상이한 유전자 좌표에서 최대값을 포함한다. 본 발명의 특정 바람직한 구체예에서, 각 전사체에 대한 확률 부함수의 최대값은 임의선택적으로, 앞서 설명된 바와 같은 변환 후 (이것은 예를 들면, 인트론의 제거 후, 변형된 유전자 커버리지를 제공할 수 있다), 단계 a)에서 전사체에 배정된 유전자 좌표의 첫 번째와 마지막 뉴클레오티드 사이에 최소한 약 1/j_max 배 차이로 분리된다. 따라서 확률 부함수는 전사체에 의해 커버되는 유전자 좌표 (및 그러므로, 커버리지를 모형화하는 확률 밀도 함수)를 따라서 균등하게 분포될 수 있다. 당연히, 최적 균등한 분포로부터 일탈하는 것이 또한 가능하다. 따라서 유전자 좌표 상에서 최고 위치는 예로서, 각 확률 부함수에 대한 균등한 분포 최대값으로부터, 전사체 길이 또는 j_max에 의해 나눠진 전사체 길이의 0%-50%, 바람직하게는 10%-40% 일탈할 수 있다.

따라서, 단계 d)에서 전사체에 대한 확률 부함수는 양성 값을 갖는 전사체의 전체 길이를 커버하도록 유전자 좌표에서 배치되거나 또는 이동될 수 있다.

바람직한 구체예에서, 상기 방법은 상기 전사체 단편 염기서열결정 데이터를 제공하기 위해, 최소한 하나의 전사체, 바람직하게는 mRNA 또는 cDNA의 서열 리드를 결정하는 것을 포함하고, 여기서 상기 리드는 상기 전사체의 단편의 서열을 포함한다. 결정 단계는 당분야에서 공지된 임의의 방법, 예를 들면, PCR 염기서열결정에 의해 수행될 수 있다. 이런 방법은 Maxam-Gilbert 염기서열결정, 사슬-종결 방법, 샷건 염기서열결정, 가교 PCR, 대량 평행 서명 염기서열결정 (MPSS), 폴로니 염기서열결정, 파이로시퀀싱, Illumina (Solexa) 염기서열결정, SOLiD 염기서열결정, 이온 반도체 염기서열결정, DNA 나노볼 염기서열결정, 헬리스코프 단일 분자 염기서열결정, 단일 분자 실시간 (SMRT) 염기서열결정, 나노포어 DNA 염기서열결정, 혼성화에 의한 염기서열결정, 질량 분광분석법을 동반한 염기서열결정, 미소유체 생어 염기서열결정, 현미경검사-기초된 기술, RNAP 염기서열결정, 시험관내 바이러스 고처리량 염기서열결정을 포함한다.

대안으로 또는 부가적으로, 본 발명의 방법은 예로서, 미리 준비된 서열 데이터를 이용하는 경우에 전사체 단편 염기서열결정 데이터를 본 발명의 방법을 수행할 수 있는 연산 장치에 공급하는 단계를 포함한다. 상기 연산 장치는 이후, 본 발명의 방법을 수행하여, 공급된 염기서열결정 데이터로 존재비 정보에 도달할 것이다.

존재비 정보는 예로서, 아래에 더욱 설명된 바와 같이 출력 장치, 예를 들면, 비디오 스크린, 프린터 또는 컴퓨터 판독가능 매체에서 표시될 수 있다.

단계 f)에서 적합시키는 것은 기대값 최대화 알고리즘에 의해 수행될 수 있다. 기대값 최대화 알고리즘은 Dempster et al., 1977에 의해 설명되었다. 당연히 임의의 다른 최대화 또는 최소화 알고리즘이 본 발명의 합계 함수를 포락선에 적합시키기 위해 조정될 수 있다. 본질적으로 가중치 알파 i와 베타 i,j - 및 임의선택적으로, 확률 부함수를 규정하는 임의의 추가 파라미터 -는 예로서, 확률 부함수의 파라미터를 변화시킴으로써 합계 함수 (모형) 및 커버리지 포락선 (단편 서열에 기초된 실제 데이터) 사이에 차이를 최소화하도록 변형된다. 특히, 계산 단계 e)와 f)를 반복함으로써 여러 반복 후 최소값을 찾기 위해, 이런 최소화 또는 변동을 만드는 많은 알고리즘이 존재한다.

확률 부함수의 이런 추가 파라미터는 예로서, 반치전폭 값 또는 특히, 가우스 함수의 경우에, 표준 편차 값 (시그마 i,j)이다. 이들은 함수 폭의 파라미터이다. 바람직하게는, 전사체 i의 각 확률 부함수에 대해 확률 부함수의 최대값으로부터 떨어진 임의의 유전자 좌표에서 함수 폭 (예를 들면, 바람직하게는 반치전폭 값 또는 표준 편차에서 전체 폭)은 대략 동일하고, 바람직하게는 동일하다. 대안으로, 폭 대 최대값의 비율은 전사체 i의 각 확률 부함수에 대해 대략 동일하고, 바람직하게는 동일하다. 이것은 전사체의 상이한 확률 분포 함수 사이에 베타에 대해 앞서 설명된 바와 유사하게 함수 폭을 "묶는다".

일반적으로, 바람직한 구체예에서, 전사체의 확률 부함수의 군 내에서 및/또는 별개의 전사체에 대해 동일한 식별자 j를 갖는 확률 부함수의 군 내에서 확률 부함수의 최소한 1, 2, 3, 4, 5, 또는 6개 파라미터가 적합 단계 f) 동안 묶인다, 예를 들면, 서로 결합되어 변형된다.

마지막 단계는 적합의 수렴이다. 수렴 규준은 오퍼레이터에 의해 세팅될 수 있다. 이것은 파라미터 α_i 또는 β_i,j 중에서 임의의 한 가지 또는 이들의 조합, 예를 들면, 모든 α들 또는 모든 β들, 또는 모든 α들과 β들의 최고 조정 또는 예로서, 실시예에서 방정식 (15)에서처럼 확률 분포 함수에서 최고 조정일 수 있다. 예를 들면, 수렴 규준은 방정식 (21)에 따른 로그 확률 증가가 0.5 미만인지 여부일 수 있다. 임의의 수렴 규준은 원하는 품질을 달성하기 위해 오퍼레이터에 의해 선택될 수 있다.

등용에서, 본 발명의 방법은 다른 편향 감소 방법과 더욱 결합된다. 예를 들면, 단계 e)의 획득된 추정된 전사체 존재비는 편향 계수에 의해 더욱 가중되고, 따라서 가중된 추정된 전사체 존재비를 획득한다. 이러한 편향 계수는 일정한 전사체가 습성-화학적 전사체 또는 단편 산출 동안 증가되거나 또는 감소된다는 것을 고려할 수 있다. 일부 단편, 예를 들면, C로 시작하는 단편은 지나치게 많을 수 있고, 그리고 바람직하게는 1보다 낮은 편향 계수에 의해 가중되고, 따라서 단계 e)의 추정된 전사체 존재비와 비교하여, 가중된 추정된 전사체 존재비를 감소시킨다. C로 시작하는 단편이 지나치게 많은 화학적 이유는 프라이머 동역학에 있다: 프라이머는 G 또는 T로 시작하는 주형과 비교하여 C로 시작하는 주형과 더욱 잘 어닐링한다. 이런 이유로, G 및/또는 T로 시작하는 단편은 C로 시작하는 단편에 대한 편향 계수보다 큰 편향 계수에 의해 가중될 수 있다. 일반적으로, 숙련자는 체계적인 편향이 존재하는 임의의 파라미터에 대해 이런 편향 계수를 획득할 수 있는데, 이들은 임의의 다중핵산 성질에 의존할 수 있고, 그리고 따라서, 단계 3)의 결과를 가중하는 합치하는 편향 계수를 포함한다.
본 발명은 본 발명에 따른 방법을 컴퓨터에서 수행하기 위한 또는 본 발명에 따른 방법을 컴퓨터에 의해 뒷받침하기 위한 컴퓨터 프로그램 제품을 포함하는 컴퓨터 판독가능한 기억 장치에 더욱 관계하고, 특히 단계 a), b), b2), c), d), e), f)와 g) 중에서 임의의 한 가지가 컴퓨터에서 수행될 수 있다. 임의의 창의적인 방법 또는 단계가 컴퓨터-실행된 방법으로서 수행될 수 있다. 심지어, 서열 리드를 결정하는 통상적으로 습성-화학 단계는 예로서, 자동화된 또는 반자동화된 서열 리더로부터 데이터를 제어하고 획득하기 위해 컴퓨터에 의해 보조될 수 있다. 컴퓨터 프로그램 제품 또는 기억 장치는 또한, 표본으로부터 짧은 리드를 획득하는 리드 산출 성분, 예를 들면, 서열분석기, 바람직하게는 컴퓨터 성분을 포함하는 서열분석기가 제공될 수 있다. 예를 들면, 컴퓨터 판독가능 매체는 자성 저장 장치 (예를 들면, 하드 디스크, 플로피 디스크, 자성 스트립, ...), 광학 디스크 (예를 들면, 콤팩트 디스크 (CD), 디지털 다용도 디스크 (DVD), ...), 스마트 카드, 그리고 플래시 메모리 장치 (예를 들면, 카드, 스틱, 키 드라이브, ...)를 포함할 수 있지만 이들에 한정되지 않는다.

삭제

본 발명은 본 발명의 이들 구체예에 한정됨 없이, 다음의 도면과 실례에 의해 더욱 예시된다 - 각 원소는 본 발명의 임의의 다른 구체예와 병합가능하다. 특히, 아래에 제공된 공식 중에서 임의의 한 가지는 본 발명의 방법에서 단계를 설명하거나 또는 규정하기 위해 별개로, 개별적으로 또는 조합으로 이용될 수 있다.

도면:
도면 1: cDNA 산출, 단편화, 염기서열결정, 서열 리드를 참고 서열에 매핑, 그리고 단편 서열의 매핑된 데이터 세트의 분석의 단계 이후에 NGS 작업 흐름;
도면 2: 2개 전사체의 커버리지 히스토그램의 혼합체;
도면 3, 참고 도면: Cufflinks 모형에 의해 가정된 단편 시작 부위의 분포;
도면 4, 참고 도면: Cufflinks 모형에 의해 가정된 커버리지;
도면 5: Uqcrq 유전자 상에서 전사체 리드의 커버리지;
도면 6: 유전자 좌표 변환; 평균 확률 p_(avg)(r)의 환산되고 이동된 이형에 의한 전체 확률 p_(total)(r)의 분해; 전체 존재비 포락 (위쪽), 전사체 1에 대한 존재비 (중앙), 전사체 1에 대한 변형된 유전자 커버리지를 예증하기 위해 접합부 (이동)를 제거하고 환산함으로써 전사체 1에 대한 초기 유전자 커버리지의 변환 (아래쪽);
도면 7: 유전자 좌표 상에서 동등한 존재비 (각각 전체의 1/3)를 갖는 3개의 모형 전사체의 엑손 구조와 접합부;
도면 8: 접합부 (이동)의 제거 후, 도면 7에서 전사체에 대한 전사체 좌표에서 시작 부위 확률 분포 p(s(r)|t=i);
도면 9: 도면 7에서 전사체에 대한 유전체 좌표에서 시작 부위 확률 분포 p(s(r)|t=i);
도면 10: i=1,2,3에 대해 동등한 존재비 α_i = 1/3을 갖는, 도면 7과 8에서 전사체에 대한 시작 부위의 전체 확률 분포;
도면 11: 도면 7에서 전사체의 신장된 시작과 끝 부위. 모형은 환산과 이동 파라미터 람다 i와 뉴 i로 정확한 시작과 끝 부위 분포를 추정한다;
도면 12: 시그마의 묶기: 확률 부함수에서 임의선택적인 추가 표준 편차 (폭) 파라미터인 시그마는 전사체 i에 대한 각 확률 부함수 내에서 일정하게 유지되는데, 합계 함수 (아래쪽)와 함께 2개의 전사체 1 (중앙)과 2 (위쪽)가 도시된다; 파선은 전사체의 유전자 커버러지 범위 전체에 대하여, 만약 그렇지 않으면 균등하게 분포된 확률 부함수 (여기서는 가우스 함수)의 합계 함수를 모형화하기 위한 확률 부함수를 결여한다;
도면 13: 도면 10에서 전체 단편 시작 부위 분포로부터 표본추출된 시작 부위의 히스토그램;
도면 14: 표본추출된 단편 길이의 히스토그램;
도면 15: 가우스의 혼합체의 성분 (확률 부함수) (점선); 실선: 전사체에 대한 확률 분포 함수;
도면 16: 1tie-Mix² 모형에 대한 가중 인자 알파의 수렴;
도면 17: 1tie-Mix² 모형에 대한 가중 인자 베타의 수렴;
도면 18: 1tie-Mix² 모형에 대한 최종 가중 인자 베타;
도면 19: 1tie-Mix² 모형에 대한 전사체 2 [p(s(r)|t=2)]의 최종 확률 분포 함수 (pdf);
도면 20: 1tie-Mix² 모형에 대한 합계 함수 [p_(total)(s(r))]의 수렴;
도면 21: 수렴 후, 알파₍₃₎=0.2에서 1tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 22: 수렴 후, 알파₍₃₎=0.4에서 1tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 23: 수렴 후, 알파₍₃₎=0.6에서 1tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 24: 수렴 후, 알파₍₃₎=0.8에서 1tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 25: 정확한 및 부정확한 전사체 주해의 엑손 구조;
도면 26: 5tie-Mix² 모형에 대한 알파_i의 수렴;
도면 27: 5tie-Mix² 모형에 대한 베타_j의 수렴;
도면 28: 5tie-Mix² 모형에 대한 뮤_j (확률 부함수 j의 이동에 대한 파라미터)의 수렴;
도면 29: 5tie-Mix² 모형에 대한 시그마_j (확률 부함수 j에 대한 폭 파라미터)의 수렴;
도면 30: 5tie-Mix² 모형에 대한 뉴_i (전사체 i에 대한 확률 분포 함수에 대한 이동 또는 번역 인자)의 수렴;
도면 31: 5tie-Mix² 모형에 대한 람다_i (전사체 i에 대한 확률 분포 함수에 대한 환산 인자)의 수렴;
도면 32: 전사체 1과 5tie-Mix² 모형에 대한 확률 분포 함수 p(r|t=i)의 수렴; 수직의 실선은 전사체의 정확한 시작과 끝 부위를 지시한다;
도면 33: 전사체 2와 5tie-Mix² 모형에 대한 확률 분포 함수 p(r|t=i)의 수렴; 수직의 실선은 전사체의 정확한 시작과 끝 부위를 지시한다;
도면 34: 전사체 3과 5tie-Mix² 모형에 대한 확률 분포 함수 p(r|t=i)의 수렴; 수직의 실선은 전사체의 정확한 시작과 끝 부위를 지시한다;
도면 35: 5tie-Mix² 모형에 대한 합계 함수 p_(total)(r)의 수렴;
도면 36: 수렴 후, 알파₍₃₎=0.2에서 5tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 37: 수렴 후, 알파₍₃₎=0.4에서 5tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 38: 수렴 후, 알파₍₃₎=0.6에서 5tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 39: 수렴 후, 알파₍₃₎=0.8에서 5tie-Mix² 모형, Cufflinks 모형에 대한 추정된 알파_i와 참 알파_i (존재비, 전체의 비율);
도면 40: 수렴 후, 알파₍₃₎=0.2에서 5tie-Mix² 모형, 1tie-Mix² 모형에 대한 추정된 알파_i와 참 알파_i;
도면 41: 수렴 후, 알파₍₃₎=0.4에서 5tie-Mix² 모형, 1tie-Mix² 모형에 대한 추정된 알파_i와 참 알파_i;
도면 42: 수렴 후, 알파₍₃₎=0.6에서 5tie-Mix² 모형, 1tie-Mix² 모형에 대한 추정된 알파_i와 참 알파_i;
도면 43: 수렴 후, 알파₍₃₎=0.8에서 5tie-Mix² 모형, 1tie-Mix² 모형에 대한 추정된 알파_i와 참 알파_i.

실시예 :

실시예 1: NGS 방법의 소개

어떤 mRNA 분자가 본래 표본에서 존재하였는지를 추론하기 위해, NGS 리드는 공지된 방법, 예를 들면, 버로우즈-휠러 변환으로 참고 유전체 위에 매핑된다. 각 리드에 대해 이것은 스플라이스 부위에 관한 정보를 잠재적으로 포함하는 한 세트의 유전자 좌표를 제공한다. 매핑 과정은 도면 1에서 가시화된다. 여기에서 서열분석기에 의해 생산된 짧은 리드의 위치가 참고 유전체 내에서 확인된다. 이러한 과정은 서열분석기에 의해 산출된 모든 리드에 대해 반복되는데, 이것은 도면 1에서 채워진 흑색 곡선 아래에 짧은 직선에 의해 지시된 바와 같이, 유전자 축 상에서 다수의 짧은 서열을 유발한다. 매핑된 리드의 통합된 통계는 유전자 축 상에서 상이한 유형의 히스토그램 (즉, 상이한 유형의 커버리지 포락선)을 야기한다. 예를 들면, 도면 1에서 채워진 흑색 곡선이 이러한 커버리지 (포락)을 묘사한다. 유전자 축 상에서 소정의 위치에서 곡선의 값은 상기 위치를 커버하는 리드의 숫자이다. 다른 히스토그램 역시 다른 유형의 커버리지 포락선인 단편 시작 부위 히스토그램과 유사하게 연구된다. 각 유전자 위치에서 이러한 히스토그램의 값은 이러한 위치에서 시작하는 리드의 숫자이다. 본 발명의 방법은 특수한 유형의 히스토그램에 의존하지 않고, 유전자 축 및 유전자 좌위 내에 단편의 세트 둘 모두에서 임의의 유형에 적용가능하다.

유전자 좌위, 예를 들면, 유전자는 다수의 상이한, 아마도 중복 영역을 내포할 수 있는데, 이들은 각각, 그들의 자체 mRNA를 산출한다. 이런 영역은 전사체로 불린다. 유전자 좌위에서 히스토그램 h_(locus)(r)은 이런 이유로, 좌위 내에 각 개별 전사체와 연관된 히스토그램의 혼합체, 다시 말하면,

이고,

여기서 h(r|t=i)은 전사체 t=i와 연관된 히스토그램이고, 그리고 i는 전사체의 명칭, 예를 들면, Apoe이다. 도면 2는 도면 1에서 히스토그램의 2개 전사체의 히스토그램으로의 가능한 분해를 보여준다. 전사체 1은 유전자 좌위의 전체 길이에 걸쳐 신장하고 전사체 2의 왼쪽으로 시작한다. 이에 더하여, 전사체 1은 2개 엑손 및 1개 접합부를 갖는다. 전사체 2는 다른 한편, 전사체 1의 첫 번째 엑손의 단부에서 끝나는 단일 엑손을 갖는다. 도면 2에서 채워진 흑색 곡선은 전사체 1의 히스토그램이고, 그리고 전사체 1의 히스토그램을 제외하고 이의 위쪽에 곡선은 전사체 2의 히스토그램 (이것은 미가공 단편 염기서열결정 데이터로부터 알려져 있지 않다)이다. 전사체 t=i의 가중 α_i는 다음에 의해 제공된다:

이것은 전사체 t=i에 속하는 전체 히스토그램 h_locus(r)에서 수치의 비율이다. h_locus(r)이 예로서, 단편 시작 부위의 히스토그램이면, α_i는 전사체 t=i에 의해 산출되었던 단편의 백분율이다. 이러한 경우에, α_i는 전사체 t=i의 존재비로 불리고, 그리고 mRNA 표본 내에서 전사체 t=i의 농도에 직접적으로 상관된다. 도면 2에서 전사체의 가중은 히스토그램 h(r|t=i) 아래 면적 대 h_locus(r) 아래 전체 면적의 비율이다. 양쪽 전사체에 대한 히스토그램의 면적이 거의 동일하기 때문에, 각 전사체에 대한 가중은 약 0.5이다. 따라서 각 전사체는 유전자 좌위 내에 위치를 커버하는 약 50% 기회를 갖는다.

방정식 (1)에서 히스토그램 h_locus(r)의 분해는 통상적으로 알려져 있지 않고, 그리고 α_i의 도출은 정교한 수학적 기구를 필요로 한다. 분해 (1)를 발견하는 문제를 수학적으로 다루기 쉽게 만들기 위해, (1)은 확률적 프레임워크, 다시 말하면,

에서 재공식화되고,

여기서 p_locus(r)은 좌위에서 단편 r을 관찰하는 전체 확률이고, p(r|t=i)는 전사체 t=i를 고려할 때 좌위에서 단편 r을 관찰하는 확률이고, 그리고 α_i는 전사체 t=i가 좌위에서 관찰되는 확률이다. 확률 분포 p_locus(r)와 p(r|t=i)가 이들을 1까지 합계에 정규화함으로써 히스토그램으로부터 도출되면, (3)은 (1)과 (2)의 직접적인 결과이다. 아래에서 개발된 방법이 일반적인 확률적 세팅에서 적용가능하고 단지 차세대 염기서열결정의 분야에만 적용가능한 것이 아니라는 사실을 강조하기 위해, 아래첨자 "locus (좌위)"는 "total (전체)"에 의해 대체될 것이고, 따라서

p_(locus)(r) = p_(total)(r)(4)이다.

전사체 확률 α_i의 추정에서 전사체 분포의 모양은 근본적으로 중요하다. 현재 방법, 예를 들면, Cufflinks (Trapnell et al., 2010)는 하지만, 이들 분포를 정확하게 모형화하지 못한다. 그 대신에, 단편 시작 부위 분포에 대해 Cufflinks는 단지 단편 길이의 분포에만 의존하는 모형을 이용한다. Cufflinks에서 단편 길이에 대한 디폴트 모형은 평균 200 및 80의 표준 편차를 갖는 가우스이다. Cufflinks에서 다른 가정과 함께, 이것은 200 bp의 리드 길이 및 3000 bp의 전사체에 대해 전사체의 단편 시작 부위 분포와 커버리지가 도면 3과 4에서 도시된 함수에 의해 모형화된다는 것을 암시한다.

대조적으로, 도면 5는 사각형 라이브러리 프렙의 커버리지를 보여주는데, 이것은 도면 4에서 Cufflinks에 의해 가정된 분포로부터 유의미하게 일탈한다. 전사체의 5' 및/또는 3' 단부를 향한 커버리지에서 편향은 현재 NGS 라이브러리 사이에서 공통 특질이고, 그리고 Cufflinks의 가정은 이런 이유로 무효하다. 결과로서, Cufflinks에 의해 추정된 전사체 존재비와 농도는 그들의 참값과 실제적으로 상이할 수 있다. 본원에서 설명된 방법은 전사체 분포와 그들의 확률을 동시에 추정하고, 그리고 이런 이유로, 전사체 존재비의 추정에서 훨씬 정확하다.

유전자 위치에서 관찰된 분포는 종종, 어떤 전사체가 유전자 위치에서 내포되는 지에 관한 시각 단서 및 그들의 확률을 표시한다. 도면 2에서 실례를 다시 한 번 고려한다. 전사체의 평균 분포가 도면 6에서 것과 유사한 모양을 갖고, 그리고 전사체 확률 p(r|t=i)이 환산과 이동에 의해 p_avg(r)로부터 대략적으로 도출되는 것으로 알려져 있으면, 도면 2에서 중첩의 성분을 발견하는 과제는 평균 분포 p_avg(r)의 모양과 유사한 전체 확률 분포 p_total(r)에서 모양을 찾음으로써 해결될 수 있다. 이러한 과정은 도면 6에서 가시화되는데, 이것은 p_avg(r)의 환산되고 이동된 이형의 중첩에 의한 p_total(r)에 최고 적합이 도면 2에서 중첩과 유사한 중첩, 따라서 전사체 확률 α_i의 정확한 추정을 야기한다는 것을 지시한다.

p(r|t=i)가 p_avg(r)의 환산되고 이동된 이형이기 때문에, 이들은 다음에 의해 제공된다:

여기서 λ_i와 ν_i는 전사체 t=i의 환산과 이동 파라미터이다. 이것은 전사체 t=i가 전사체에 독특한 단지 2개의 파라미터 λ_i와 ν_i를 갖는다는 것을 의미한다. 모든 다른 파라미터, 다시 말하면, p_avg(r)의 파라미터는 상이한 p(r|t=i) 사이에 묶인다. 이것은 상이한 전사체에 공통적인 구조, 예를 들면, p_avg(r)이 p(r|t=i)의 파라미터를 적절하게 묶음으로써 추정될 수 있다는, 아래에서 설명된 방법의 중심 개념을 강조한다. 아래에서 p(r|t=i)는 함수의 혼합체에 의해 모형화될 것이고, 그리고 (3)에서 p_total(r)의 분해는 이런 이유로, 기능의 혼합체의 혼합체이다. 가독성을 증가시키기 위해, 이러한 모형은 본 발명의 모형으로서 지칭될 것이다. 다음의 섹션은 Mix² 모형으로 또한 불리는 본 발명의 모형에 대한 일반적인 소개를 제공하고, 그리고 이의 여러 변이체를 논의한다. 이에 더하여, 실험은 본 발명의 모형, Mix² 모형이 확률 분포 p(r|t=i)를 확실하게 추정하고, 그리고 Cufflinks 모형의 것들보다 훨씬 정확한 α_i에 대한 추정치를 산출할 수 있다는 것을 보여준다.

실시예 2: 좌표 변환

2.1. 유전체와 전사체 좌표에서 위치

유전자 축은 생물체에 대해 염기서열화된 염기쌍의 서열인데, 이것은 통상적으로, 제로 또는 1에서 시작하고, 그리고 생물체의 복잡성에 따라, 수억 개 염기쌍 길이까지 이를 수 있다. 이에 더하여, 유전자 축은 통상적으로, 염색체 또는 콘틱으로 세분된다. 유전자 축은 도면 5의 위쪽에 가시화되는데, 이것은 이러한 그래픽이 대략 염기쌍 53,242,500과 53,244,200 사이에 염색체 11 상에서 유전체의 선택을 나타낸다는 것을 지시한다. 전사체는 통상적으로, 유전자 축 상에서 엑손의 서열 (엑손₁,...,엑손_N)으로서 규정되고, 여기서 i-번째 엑손은 유전자 축 상에서 간격 [s(엑손_i),...,e(엑손_i)]인데, 이것은 s(엑손_i)에서 시작하고 e(엑손_i)에서 끝난다. 2개 연속 엑손 사이에 갭 [e(엑손_i)+1,s(엑손_i+1)-1]은 인트론으로 불리고, 그리고 인트론을 선행하는 마지막 뉴클레오티드로부터 인트론 이후에 첫 번째 뉴클레오티드까지 결합은 접합부로 불린다. 전사체의 3가지 실례가 도면 7에서 도시된다. 이러한 도면에서 x 축은 유전자 좌표를 제공하는데, 이것은 1000 내지 5500 범위에서 변하고, 반면 y 축은 전사체 id를 제공한다. 따라서, 전사체 1은 위치 1000에서 시작하고 위치 2500에서 끝나는 단일 엑손으로 구성된다. 전사체 2는 다른 한편, 엑손 서열 ([1500,3200],[4000,5000])에 의해 규정되고, 반면 전사체 3은 엑손 서열 ([2700,3200],[4000,5500])에 의해 규정된다. 따라서, 전사체 2와 3은 도면 7에서 파선 화살표에 의해 표시되는 동일한 접합부를 갖는다. 아래에서, 전사체의 길이는 l(t)에 의해 표시될 것이다. 따라서 l(전사체 1)=1501, l(전사체 2)=2702 및 l(전사체 3)=2002. 도면 7에서 전사체의 엑손 서열과 길이는 또한, 표 1에서 요약된다.

도면 7에서 전사체의 주해

	엑손 1	엑손 2	길이
전사체 1	[1000,2500]	해당 없음	1501
전사체 2	[1500,3200]	[4000,5000]	2702
전사체 3	[2700,3200]	[4000,5500]	2002

전사체의 엑손 범위 안에 있는 유전자 축 상에서 지점 역시 이러한 전사체의 좌표계 내에서 참조될 수 있다. 전사체 좌표는 선행하는 인트론의 길이를 제외하고 전사체의 시작으로부터 거리이고, 그리고 이런 이유로, 1과 l(t) 사이에 숫자이다. 따라서 전사체의 엑손 e_i에서 위치 P에 대해 전사체와 유전체 좌표 P_trans는 아래와 같이 전환될 수 있다:

그리고

예를 들면, 도면 7에서 유전자 좌표에서 위치 4500을 고려한다. 이러한 위치는 두 번째와 세 번째 전사체 둘 모두의 엑손 범위 안에 있다. 전사체 2의 경우에 이러한 위치의 전사체 좌표는 4500-1500+1-801=2200이고, 반면 전사체 3의 경우에 전사체 좌표는 4500-2700+1-801=1000이다. 위치 4500이 전사체 1의 엑손 범위 안에 있지 않기 때문에, 이러한 위치를 전사체 1에 대한 전사체 좌표로 전환하는 것이 가능하지 않다. 전사체의 엑손 범위 안에 있는 유전자 축 상에서 위치는 전사체와 양립성인 것으로 불릴 것이다. 전사체와 양립하는 유전자 (예를 들면, 유전체) 좌표에서 위치의 전사체의 좌표로의 변환은 T에 의해 표시될 것이다. T가 연관되는 전사체는 전후관계로부터 명확해질 것이다. 전사체 좌표에서 확률 분포 p_trans(P_trans|t=i)는 아래와 같이 유전자 또는 유전체 좌표에서 확률 분포로 전환될 수 있다:

표기법을 단순화하기 위해, "genome (유전체)"와 "트랜스 (trans)" 아래첨자는 아래에서 편의할 때는 언제든지 누락될 것이다. 도면 8과 9는 전사체 좌표와 유전체 좌표에서 확률 분포 사이에 관계를 보여준다. 도면 8은 도면 7에서 3개의 전사체 각각에 대한 환산된 확률 분포를 보여준다. 혼합체의 성분으로서 이들의 용도를 기대하면서, 각 확률 분포는 1/3의 인자로 곱셈되었다. 실선은 3개 전사체 중에서 가장 짧고 1501개 염기쌍의 길이를 갖는 전사체 1에 속한다. 전사체 2와 3은 다른 한편, 2002개와 2702개 염기쌍 길이를 갖고, 그리고 전사체 좌표에서 이들의 확률 분포가 이런 이유로, 각각 1 내지 2002와 2702 범위에서 변한다.

도면 9는 유전체 좌표에서 전사체 확률 분포를 보여준다. 여기에서 전사체의 상이한 시작 부위와 접합부가 고려되었다.

유전체 좌표에서 전사체 확률의 합계를 취하는 것은 관심되는 전체 유전자 좌위에서 위치에 대한 전체 확률 분포 p_total(r)를 산출하는데, 이것은 도면 10에서 도시된다. 도면 10에서 곡선은 NGS 데이터로부터 유전자 좌위에서 도출되는 분포의 평탄화된 이형과 유사하다. 추후 섹션에서, 도면 10에서 것과 같은 분포가 단편 시작 부위로서 해석되는 유전자 좌위에서 표본을 산출하는데 이용될 것이다. 이들 표본은 이후, 본 발명의 모형의 확률 밀도 함수의 파라미터를 훈련하는데 이용될 것이고, 그리고 결과의 p(r|t=i)와 p_total(r)은 도면 9와 10에서 곡선의 정확한 추정값인 것으로 나타날 것인데, 이들은 차후에, 존재비를 가중 인자 α_i로서 산출한다.

2.2 유전체와 전사체 좌표에서 단편

단편은 전사체 내에서 연속 서열이다. 전사체와 유사하게, 단편 r은 이런 이유로, 유전자 축 상에서 간격의 서열, r=(rint₁,...,rint_K)로 구성되고, 여기서 rint_i=[s(rint_i),e(rint_i)]는 시작 s(rint_i) 및 끝 e(rint_i)를 갖는 i-번째 간격이다. 단편 r은 이의 시작과 끝이 전사체의 엑손 내에 위치되면 전사체와 양립하고, 그리고 인접한 간격 사이에 갭 [e(rint_i)+1,s(rint_i+1)-1]은 전사체의 인트론이다, 다시 말하면,

s(int₁) ∈ 엑손_i (9)

e(int_K) ∈ 엑손_k (10)

일부 i ≤ k, 그리고

단편 r이 전사체와 양립하면, 이것은 전사체 좌표에서 간격 [T(s(rint₁)),T(e(rint_K))]으로 전환될 수 있다. 단편의 간격 사이에 갭은 이런 이유로, 전사체 좌표로의 변환 동안 제거된다. 이전 섹션에서처럼, 단편 r의 전사체 좌표로의 변환은 T에 의해 표시될 것이다. 예를 들면, 단편 ([2000,3000],[4000,4500]) 및 도면 7에서 3가지 전사체를 고려한다. 이러한 단편의 시작과 끝은 각각, 전사체 2의 첫 번째와 두 번째 엑손에 있다. 동시에, 상기 단편의 갭 [3001,3999]은 전사체 2의 인트론과 일치한다. 따라서 상기 단편은 전사체 2와 양립하고 전사체 좌표에서 간격 [501,2200]으로 변환된다. 다른 한편, 전사체 1과 전사체 3의 경우에 단편의 시작 또는 끝 중에서 어느 한쪽은 이의 엑손 중에서 한 가지 범위 안에 있지 않고, 이런 이유로 상기 단편은 이들 전사체 중에서 어느 한쪽과 양립하지 않는다.

전사체 좌표에서 단편 상에서 확률 분포 p_trans(r_trans|t=i)는 이런 이유로, 아래와 같이 유전체 좌표에서 단편 상에서 확률 분포로 전환될 수 있다.

앞서와 같이, "genome (유전체)"와 "trans (트랜스)" 아래첨자는 편의할 때는 언제든지 누락될 것이다. 이의 시작과 끝을 이용하는 대신에, 간격은 또한, 이의 시작 s(r)과 길이 l(r)에 의해 대표될 수 있다. 이것은 전사체 좌표에서 확률 분포의 하기 인수분해를 허용한다.

그리고

이들 인수분해는 NGS 데이터에 편의한데, 그 이유는 단편 길이의 전역 분포가 통상적으로, 라이브러리의 바이오분석기 추적으로부터 추론될 수 있기 때문이다. (13)에서 인수분해가 Cufflinks에 의해 이용되는데, 이것은 소정의 길이의 단편이 전사체 상에서 동등한 확률로 어디든지 배치될 수 있다는 것을 부가적으로 가정한다. 대조적으로, 섹션 4에서 실험에서 Mix² 모형은 (14)에서 인수분해를 이용하는데, 이것은 단편 시작 부위의 전사체 특정한 가변성의 더욱 효율적인 추정을 허용한다.

실시예 3: 전사체 확률과 전사체 특정한 확률 분포의 추정

아래에서 설명되는 모형은 함수의 혼합체의 혼합체를 이용하고, 그리고 이런 이유로, Mix² 모형으로 불릴 것이다.

3.1 Mix² 모형의 수학적 기초

아래에서, r은 단편과 위치 둘 모두를 나타낼 수 있다. 하지만, 편의를 위해, r은 항상 단편으로서 지칭될 것이다. 유전자 좌위 p_total(r)에서 특정 단편 r을 관찰하는 확률은 전사체가 단편을 산출하는 확률에 의해 가중된 전사체에 대한 단편을 관찰하는 확률의 합계이다. 따라서 p_total(r)은 확률 분포의 하기 혼합체에 의해 제공된다.

섹션 2에서 설명된 바와 같이, r이 t=i와 양립하면,

이고, 그리고 만약 그렇지 않으면

이다. 상기 방법은 확률 분포 p _trans (r|t=i)가 혼합체, 다시 말하면,

이라고 가정한다.

여기서 M_i은 혼합체 성분의 숫자이고, 그리고 p_trans(r|t=i,b=j)는 확률 분포이다. β_i,j ≥ 0은 1까지 합계하는 가중, 다시 말하면,

이다.

앞서와 같이, 변수 t는 전사체를 나타내는 것으로 의미되고, 반면 새로 도입된 변수 b는 일반적으로 확률 부함수로서 또한 지칭되는 "빌딩 블록"을 나타내는 것으로 의미된다. 이러한 모형에서 전사체 t=i 및 블록 b=j가 단편 r을 고려할 때 관찰되었던 사후 확률을 계산하는 것이 가능하다. 이러한 사후 확률은 다음에 의해 제공된다:

주의할 점은 단편 r이 전사체 t=i와 양립하지 않으면,

및 따라서,

이라는 것이다. t와 b에 대한 사후 확률은 다음에 의해 제공된다:

그리고

상기 모형은 p_total(r)의 파라미터 중에서 일부가 상이한 전사체 (확률 분포 함수)와 블록 (확률 부함수) 사이에 묶인다는 것을 더욱 가정한다. 이것은 혼합체 가중치 β_i,j 뿐만 아니라, 확률 분포 p_trans(r|t=i,b=j), i=1,...,N, j=1,...,M의 확률 부함수의 세트를 결정하는 임의의 파라미터 Θ 둘 모두를 포함할 수 있다. 이러한 묶기는 통합된 또는 전사체와 블록의 군 내에서 모든 전사체와 블록에 적용될지도 모른다. 확률 분포 p_trans(r|t=i)의 파라미터를 묶는 것이 p_trans(r|t=i) 사이에 일정한 유사성을 암시하기 때문에, 파라미터 묶기는 파라미터 묶기의 유형에 의해 암시된 유사성의 종류를 표시하는 전사체에만 적용되어야 한다. 예로서, 파라미터 묶기가 p_trans(r|t=i)가 도입부에서 제안되었던 바와 같이, 환산과 이동에 의해 단일 p_avg(r)로부터 도출된다는 것을 암시하면, 파라미터 묶기는 p_trans(r|t=i)가 서로의 환산되고 이동된 이형인 전사체에만 적용되어야 한다. 일부 라이브러리 프렙의 경우에 후자는 일정한 범위에서 길이를 갖는 전사체에만 유효하다. 길이가 상이한 범위 안에 있는 전사체는 이런 이유로, p(r|t=i)이 서로의 환산되고 이동된 이형일 것을 필요로 하는 파라미터 묶기에 적용되지 않아야 한다.

유전자 영역에서 단편 r의 데이터 세트 R에 대해, Mix² 모형의 전사체 특정한 분포 p(r|t=i)와 전사체 확률 α_i은 모형 p_total(r) 하에 데이터 세트 R의 우도를 최대화함으로써 추정된다. 파라미터 α_i,β_i,j,Θ의 우도 함수는 이런 이유로, 다음에 의해 제공된다:

(21)의 최대화는 근사 해답이 다수의 상이한 최적화 방법으로 발견될 수 있는 제약된 비선형 최적화 문제이다. 아래에서 이용된 최적화 방법은 기대값 최대화 (EM) 알고리즘이다. EM 알고리즘은 우도 함수의 국부 최대값을 찾는 반복적인 절차이다. 전역 최적값에 가까운 국부 최대값을 획득하기 위해, 모형 파라미터에 대한 상이한 초기화로 EM 알고리즘을 여러 번 시도하는 것이 필요할 수도 있다. 가장 높은 우도를 산출하는 결과는 이후, 최적화 문제에 대한 해답으로서 선택된다. EM 알고리즘은 p_trans(r|t=i) 또는 이들의 파라미터의 묶기의 형태에 관한 어떤 가정도 없이 α_i를 추정하는데 이용될 수 있다. Cufflinks에 의해 또한 이용되는, α_i에 대한 EM 갱신 공식은 다음에 의해 제공된다:

여기서 |R|은 데이터 세트 R에서 단편의 숫자이고, α_i ⁽ⁿ⁺¹⁾은 n+1-번째 반복 후 α_i에 대한 추정값이고, 그리고 p⁽ⁿ⁾(t=i|r)은 n-번째 반복으로부터 파라미터에서 추정된, 단편 r을 고려할 때 전사체 t=i를 관찰하는 사후 확률이다. 이러한 경우에 최적화 문제가 오목하기 때문에, 상이한 초기화로 EM 알고리즘을 반복하는 것이 필요하지 않다. 오목 문제에 대해 EM 알고리즘이 항상 수렴하는 단일 전역 최적값이 있다. 이들 모형 파라미터는 모형 p_total(r)의 관찰된 데이터 R로의 적합을 최적화함으로써 추정된다. 하지만, 궁극적 목적은 p_trans(r|t=i)에 우수한 적합, 그리고 이런 이유로, 전사체 확률 α_i의 우수한 추정값을 찾는 것이다. p_total(r)의 최적화는 그 자체로, α_i의 우수한 추정값을 암시하지 않는다. p_trans(r|t=i)의 파라미터가 적절하게 묶이는 경우에만 p_total(r)의 최적화가 전사체 확률 α_i의 우수한 추정값을 산출할 것이다. 다음 섹션은 상이한 묶기 전략을 이용하는 Mix² 모형의 여러 변이체를 소개한다. 이들 변이체 중에서 2가지가 섹션 4에서 연구될 것인데, 여기서 이들은 환산되고 이동된 전사체 확률 분포 p(r|t=i)의 경우에, Cufflinks 모형보다 α_i에 대한 훨씬 정확한 추정값을 산출하는 것으로 나타날 것이다.

3.2 Mix² 모형의 변이체

마지막 섹션에서 Mix² 모형의 수학적 기초는 상당히 일반적이다. 본 섹션은 Mix² 모형의 다수의 구체적인 실현을 논의하고, 그리고 이들의 다양한 이점과 한계를 강조한다.

3.2.1 파라미터의 단일 군 내에서 묶기

본 섹션에서 논의된 가장 단순한 Mix² 모형은 상이한 전사체 사이에 i=1,...,N에 대한 가중치 β_i,j만을 묶는다. 따라서 이러한 모형의 파라미터의 세트는 {α_i,β_j: i=1,...,N, j=1,...,M}이고, 그리고 p_trans(r|t=i)는 다음에 의해 제공된다:

이러한 Mix² 모형은 파라미터의 단일 군 내에서 파라미터, 다시 말하면 β_j만을 묶기 때문에, 이러한 모형은 아래에서 1tie-Mix² 모형으로 불릴 것이다. 1tie-Mix² 모형 하에 단편 r을 관찰하는 확률은 다음에 의해 제공된다:

방정식 (24)에서 합계의 순서는 교체될 수 있고, 이것은 다음을 산출한다:

앞서와 같이, (26)에서 피가수는 단편 r이 전사체 t=i와 양립하지 않으면 제로이다. p(r|b=j)가 확률 분포이기 때문에, 방정식 (25)는 α_i를 고려할 때, β_j를 추정하는 문제가 α_i의 추정과 개념적으로 동일하다는 것을 보여준다. 따라서 EM 알고리즘이 적용될 수 있고 β_j에 대한 하기 갱신 공식을 제공한다:

이것은 1tie-Mix² 모형의 α_i와 β_j 둘 모두 공식 (22)와 (27)의 반복적인 적용에 의해 동시에 훈련될 수 있다는 것을 암시한다. 원칙적으로, 1tie-Mix² 모형에 대한 EM 알고리즘은 α_i와 β_j의 초기화에 따라 상이한 해답에 수렴할 수 있는데, 그 이유는 Cufflinks 모형과 비교하여, 1tie-Mix² 모형의 우도 함수가 오목하지 않기 때문이다. 하지만, 섹션 4에서 실험에서, α_i와 β_j의 단일 초기화는 만족스러운 결과에 수렴하기 위해 충분하였다. 1tie-Mix² 모형에 대한 확률 분포 p_trans(r|t=i,b=j)는 데이터를 적합시키도록 선택되어야 한다. 따라서 데이터의 구조에 관한 이전 지식이 필요하다. 관찰된 p_trans(r|t=i)가 상이한 t=i에 대한 서로의 환산된 이형이면, p_trans(r|t=i,b=j) 역시 상이한 t=i에 대한 환산된 이형이어야 한다. 이에 더하여, 동일한 b=j, 하지만 상이한 t=i에 대해 p_trans(r|t=i,b=j)는 동일한 환산 인자 β_j에 의해 조절되는 전사체 좌표계의 영역 안에 있어야 한다. 후자는 전사체의 정확한 시작과 길이가 알려져 있으면 달성될 수 있다. 이러한 경우에 p_trans(r|t=i,b=j)는 t=i의 전사체 좌표계를 따라서 등거리 위치에서 배치될 수 있다. 정확한 시작과 길이가 알려져 있지 않으면, p_trans(r|t=i,b=j)를 어떻게 위치시키는 지를 즉각적으로 알 수는 없다. 다음 섹션에서 모형은 데이터로부터 전사체 시작과 길이를 학습함으로써 이러한 제한을 방지한다.

3.2.2 파라미터의 5가지 군 내에서 묶기

마지막 섹션에서 모형은 전사체의 정확한 시작 부위와 길이에 관한 지식에 의존한다. 본 섹션은 p _trans (r|t=i,b=j)를 자동적으로 위치시키고 환산하는 모형을 논의한다. 이러한 모형이 파라미터의 5가지 군 내에서 파라미터를 묶기 때문에, 이것은 아래에서 5tie-Mix² 모형으로 불릴 것이다. 5tie-Mix² 모형은 분포 p _trans (r|t=i,b=j)에 대한 가우스를 이용하는데, 이의 내부 파라미터, 다시 말하면 이들의 평균 μ _i,j 와 표준 편차 σ _i.j 는 파라미터 θ의 세트 =

: i = 1 , . . . , N, j = 1, . . . , M}로부터 도출된다. 특히 μ _i,j 와 σ _i.j 는 다음에 의해 제공된다:

와

는 상이한 전사체 사이에 묶이고, 따라서

= μ _j 및

= σ _j 이고, 반면

와

는 상이한 블록 사이에 묶이고, 따라서

= λ _i 및

= ν _i 이다. 결과적으로 방정식 (28)과 (29)는 다음과 같이 환산된다:

전사체 t = i 및 블록 b = j에 대한 pdf는 이런 이유로, 다음에 의해 제공된다:

앞서와 같이, β _i.j 는 전사체 i = 1 , . . . , N 사이에 블록 j에 대해 묶이고, 따라서 β _i.j = β _j 이다. 확률 분포

(여기서

은 평균 μ _j 와 표준 편차 σ _j 를 갖는 가우스이다)는 전사체 t = i에 대한 pdf p(r|t=i)를 획득하기 위해 ν _i 에 의해 이동되고 λ _i 에 의해 환산되는 좌위에서 전사체의 평균 pdf P _avg (r)로서 해석될 수 있다. 상기 모형의 힘을 제한하지 않으면서, 선택되는 전사체에 대해 λ _i = 1 및 ν _i = 0을 세팅하는 것이 가능하다. 방정식 (32)에서 pdf는 전사체 t = i의 좌표계에서 확률 분포가 아니지만 실제 연속적 공간에 있다. 따라서 결과의 pdf의 p(r|t=i)와 P _total (r)은 또한 연속적 pdf이다. 이러한 점을 강조하기 위해, 이들 확률 분포는 아래에서

,

및

에 의해 표시될 것이고, 그리고

의 가중치는

에 의해 표시될 것이다. 이들 연속 확률 분포로부터, 유전체와 전사체 좌표에서 P _total (r)과 p _trans (r|t=i)은

을 전사체 t = i에 대한 가능한 단편 F(t)의 세트에 제한하고 이들의 값을 1까지 합계에 정규화함으로써 계산될 수 있다, 다시 말하면,

여기서 (33)에서 분모의 합계에서 병합은 좌위 내에 모든 전사체에 걸쳐 확장된다. α_i는 이후, 아래와 같이

로부터 도출될 수 있다:

방정식 (35)은 연속 모형 가중치 알파_(i,R)의 재환산에 상당한다. 재환산의 효과는 천천히 변하는 확률 분포 p_trans(r|t=i) 및 모든 전사체가 합리적인 길이를 갖는 유전자 좌위의 경우에 별로 확연하지 않다. 따라서 추후 설명된 실험에서 재환산 공식의 효과는 미미하다. 좌위와 p_trans(r|t=i)가 이들 요건에 합치하는 상황에서 알파_(i,R)은 이런 이유로, 전사체 확률 α_i로서 직접적으로 이용될 수 있는데, 이것은 섹션 4에서 행위될 것이다. λ_i는 전사체의 상이한 길이를 설명하는 환산 파라미터이고, 반면 ν_i는 전사체 주해로부터 전사체 시작 부위의 상이한 상대적 오프셋을 설명한다. 이들 2가지 파라미터가 전사체 주해에서 부정확한 시작 부위와 길이를 교정할 수 있기 때문에, 시작과 끝의 정의를 모두 함께, 전사체 주해로부터 제거하는 것이 가능하다. 예를 들면, 모든 전사체의 시작 부위가 동일한 값에 세팅될 수 있는데, 이것은 좌위 내에 전사체의 모든 시작 부위보다 작다. 유사하게, 모든 전사체의 끝 부위가 동일한 값에 세팅될 수 있는데, 이것은 좌위 내에 전사체의 모든 끝 부위보다 크다. 전사체 주해의 이러한 확장은 도면 11에서 예시된다. 전사체 주해를 확장하는 것은 전사체 시작 또는 끝과의 비적합성으로 인해 전사체에 대해 이전에 무효하였던 단편이 이러한 전사체에 대해 이제부터 유효해지는 효과를 갖는다. 도면 11에서, 예로서, 전사체 2에 대해 유효한 모든 단편은 또한, 전사체 3에 대해 유효하다. 이에 더하여, 인트론에 걸치거나 들어가지 않는 모든 단편은 3가지 전사체 모두에 대해 유효하다.

본 섹션에서 모형의 파라미터는 EM 알고리즘으로 효율적으로 훈련될 수 있다. 평균, 분산과 오프셋 파라미터에 대한 EM 갱신 공식은 아래에 제공된다.

환산 파라미터 λi에 대하여 EM 알고리즘의 보조 함수의 유도체를 계산하는 것은 하기 이차 방정식을 야기한다.

대부분의 실질적인 상황에서 방정식 (39)는 1개의 양성과 1개의 음성 해답을 갖는다. λ_i이 환산 파라미터이기 때문에, 이것은 양성이어야만 하고, 그리고 이런 이유로, (39)의 단지 양성 해답만 관심된다. 섹션 4에서 실험은 5tie-Mix² 모형이 부정확한 전사체 주해를 실제로 교정할 수 있고, 그리고 이러한 경우에, 정확한 전사체 주해에 의존하는 모형보다 α_i에 대한 훨씬 정확한 추정값을 산출한다는 것을 보여준다.

3.2.3 파라미터의 2가지 군 내에서 묶기

이전 모형에서와 같이, 본 섹션에서 모형은 블록 j에 대한 전사체 전체에 대하여 β_i.j 를 묶고, 그리고 이런 이유로 β_i.j = β_j 이다. 이에 더하여 본 섹션에서 모형은 5tie-Mix² 모형으로부터

를 묶지만, 소정의 블록에 대한 전사체 전체에 대해서, 더욱 일반적으로 (i, j)의 쌍의 세트, i = 1, . . . , N, j = 1, . . . , M 전체에 대해서는 그렇지 않다. 예를 들면, N ≥ 2, M ≥ 3이면, 가능한 세트는 L = {(1, 1), (1, 2), (2, 3)}이고, 따라서

이다. 본 섹션에서 모형이 2가지 군 내에서 파라미터만을 묶기 때문에, 이것은 2tie-Mix² 모형으로 불릴 것이다. 따라서 2tie-Mix² 모형의 경우에

이다.

2tie-Mix² 모형에 대한 동기는 다음과 같다. 도면 12에서 묘사된 상황을 고려한다. 이러한 도면에서 임의의 접합부가 없는 2개의 전사체는 1000 bp에서 시작하고 4000 bp에서 끝나는 영역에 있다. 두 번째 전사체는 완전한 영역에 걸치고, 반면 첫 번째 전사체는 위치 1500에서 시작하고 위치 3500에서 끝난다. 각 전사체의 pdf p(r|t=i)는 도면 12에서 가우스 위에 2개 곡선의 모양을 갖는 것으로 추측된다. 하지만, 도면 12의 아래쪽에서 묘사되는 전체 pdf P _total (r)은 이들 두 전사체에 대한 p(r|t=i)의 정확한 합계가 아니다. 그 대신에, P _total (r)은 2250 bp와 2750 bp 사이에 영역에서 딥 (dip)을 갖는다. 이것은 NGS 염기서열결정의 공통 특질인데, 여기서 일정한 영역 내에서 뉴클레오티드 서열이 이러한 영역을 교차하는 모든 전사체에 대한 단편의 산출에 영향을 줄 수 있다. 이것은 염기서열결정을 개시하는 염기서열결정 프라이머가 일정한 서열에 결합하는 능력 없음의 결과일 수 있다. 도면 12에서 실례에서, 예로서, 딥의 영역 내에서 서열은 염기서열결정 프라이머에 의해 덜 선호될지도 모르고, 그리고 이러한 영역에서 전체 pdf P _total (r)은 이런 이유로 현저하게 감소된다. 이러한 문제에 대한 해답은 딥의 영역을 교차하지 않는 다른 전사체에 의해 공유될지도 모르는 β _j 의 묶기가 아니고, 상기 영역과 교차점을 갖는 블록 p(r|t=i,b=j)에 속하는

의 묶기이다. 도면 12에서 이들 블록은 파선에 의해 표시된다. 이들 모든 블록의

가 매우 큰 값에 세팅되면, 가우스 p(r|t=i,b=j)는 전체 연속 x 축에 걸쳐 확산하고, 그리고 유전체 좌위에서 전체 pdf P _total (r)에 거의 기여하지 않는다. P _total (r)에서 딥이 일어날 수 있는 영역이 사전에 알려져 있지 않으면, 서로에 가까운 블록에 속하는 모든

가 묶일 수 있다. 본 실례에서 전사체 1에 대한 2개의 인접한 블록의

는 그들 위에 놓여 있는 전사체 2의 블록의

, 예를 들면, L = {(1, 1), (1, 2), (2, 3)}과 묶일지도 모르고, 그리고 따라서, 본 섹션의 시작에서처럼,

이다.

에 대한 EM 갱신 공식은 방정식 (37)로부터 도출될 수 있고, 그리고 다음에 의해 제공된다:

상기 모형의 다른 파라미터에 대한 EM 갱신 공식은 이전 섹션에서와 동일하다. 이러한 상황에서 도면 12에서 β _j 를 갱신하는 것이 합리적이다. 다른 파라미터를 갱신하는 것은

에 대한 p(r|t=i,b=j)의 공간적 연관을 파괴할지도 모르고, 그리고 이런 이유로, 잠재적으로 유해하다.

3.2.4 파라미터의 6가지 군 내에서 묶기

5tie-Mix² 모형으로 단순한 확장은

를 환산하는 다른 파라미터

의 도입에 의해 제공된다. 특히

따라서, 5tie-Mix² 모형과 비교하여, 본 섹션에서 모형은

와

를 독립적으로 환산하고, 그리고 파라미터가 파라미터의 6가지 군 내에서 묶이기 때문에, 이것은 6tie-Mix² 모형으로 불린다. 6tie-Mix² 모형의 경우에

는 상이한 블록 사이에 동일한 전사체에 대해 묶이고, 따라서

= κ _i 이다. 6tie-Mix² 모형에서 λ _i 와 κ _i 는 평균 μ _j 와 표준 편차 σ _j 와 유사한 역할을 수행하고 하기 공식으로 훈련될 수 있는데, 이것은 방정식 (36)과 (37)과 유사한 구조를 갖는다:

5tie-Mix² 모형에 비하여 6tie-Mix² 모형의 이점은 이것이 상이한 전사체 사이에 평활도에서 차이를 모형화할 수 있다는 것이다. 하지만, 파라미터의 증가된 숫자는 또한, 나쁜 모형 수렴을 야기할 수 있고, 그리고 이런 이유로, NGS 데이터가 6tie-Mix² 모형에 의해 암시된 평활도에서 분산을 표시하는 경우에만 이용되어야 한다.

3.2.5 고차 Mix² 모형

이전 섹션에서 가우스 p_trans(r|t=i,b=j)의 평균 μ_i,j와 표준 편차 σ_i,j는 방정식 (30)과 (31)에 의해 제공된 λ의 전사체 특정한 아핀 선형 변환을 통해 μ_j와 σ_j로부터 도출되었다. 본 섹션은 이러한 개념을 일반화하고, 그리고 μ_i,j와 σ_i,j가 λ에서 다항식에 의해 산출되도록 허용한다, 다시 말하면,

에 대한 EM 갱신 공식은 하기 표현에 의해 제공된다.

나머지 모형 파라미터에 대한 EM 갱신 공식은 비선형 방정식을 야기하는데, 이들은 단지 그들 자체로, 반복적인 절차에 의해 해결될 수 있다.

실시예 4: 환산되고 이동된 전사체 특정한 확률 분포에 대한 실험

본 섹션에서 실험에서, 확률 p _trans (r|t=i,b=j)는 아래와 같이 인수분해된다:

여기서 s(r)과 l(r)은 단편 r의 시작과 길이이다. 이에 더하여, 확률

은 s(r) 및 전사체 l(r)의 길이에 의존하는 것으로 가정된다. 이런 이유로 (48)은

로 환산된다.

결과로서, 혼합체 모형 (16)은 단편 길이에 대한 확률 분포에 의해 곱셈된 시작 부위 분포 p _trans (s(r)|t=i)에 대한 혼합체 모형이다, 다시 말하면,

따라서, 상기 모형의 정확한 p _trans (r|t=i)로의 수렴은 분포

의 정확한 시작 부위 분포로의 수렴 및 α _i 의 정확한 존재비로의 수렴을 점검함으로써 평가될 수 있다.

2가지 유형의 실험이 본 섹션에서 논의된다. 첫 번째 유형은 10000개 단편의 세트에 대한 Mix² 모형의 수렴의 상술된 분석인데, 여기서 단편 시작 부위는 가중치 α₁=0.28, α₂=0.32와 α₃=0.4에서 도면 9에서 p(r|t=i)의 중첩으로부터 무작위로 뽑혔다. 이에 더하여, 각 단편 시작 부위 s(r)의 경우에 단편 길이는 s(r)에서 시작하는 단편의 가능한 길이, 다시 말하면, 1,...,l(t)-s(r)+1에 정규화된, 200의 평균 및 80의 표준 편차를 갖는 가우스로부터 무작위로 뽑혔다. 표본추출된 시작 부위와 단편 길이의 히스토그램은 도면 13과 14에서 도시된다. 목격될 수 있는 바와 같이, 시작 부위 히스토그램은 도면 10에서 분포의 모양을 거의 추종하고, 반면 단편 길이 히스토그램은 가우스의 모양을 추종한다.

본 섹션에서 두 번째 유형의 실험은 가중치 α₁,α₂,α₃의 60개 상이한 세트에 대해 첫 번째 유형을 반복하고, 그리고 Mix² 모형에 의해 추정된 α₁,α₂,α₃을 참 가중치와 비교한다. 전자 실험은 Mix² 모형 및 이들의 파라미터의 수렴에 대한 상세한 통찰력을 제공하고, 반면 후자 실험은 Mix² 모형에 의해 추정된 가중치 α₁,α₂,α₃의 전반적인 정확도에 주력한다. 60개 가중치의 세트는 α₃을 값 0.2, 0.4, 0.6, 0.8 각각에 세팅함으로써 선택되었다. α₃의 각 값에 대해, α₁의 15개 상이한 값이 0 및 1-α₃ 사이에 등거리 간격에서 선택되었고, 그리고 α₂는 α₂=1-α₁-α₃에 세팅되었다. 가중치의 각 세트에 대해 쿨백 라이블러 (KL) 발산이 참 가중치와 추정된 가중치 사이에서 계산되었고, 그리고 가중치의 모든 세트에 걸쳐 KL 발산의 평균이 Cufflinks와 Mix² 모형의 정확도를 정량하는 척도로서 이용되었다.

4.1 1tie-Mix² 모형

4.1.1 α₁=0.28, α₂=0.32, α₃=0.4에 대한 1tie-Mix² 모형의 수렴

앞서 언급된 바와 같이, 본 섹션에서 모형을 추정하는데 이용된 데이터는 가중치 α₁=0.28, α₂=0.32와 α₃=0.4에서 도면 9에서 p(r|t=i)의 중첩으로부터 단편 시작 부위를 표본추출함으로써 산출되었고, 반면 단편 길이는 200의 평균 및 80의 표준 편차를 갖는 재정규화된 가우스로부터 표본추출되었다. 이들 표본추출된 데이터 세트의 히스토그램은 도면 13과 도면 14에서 도시된다. 본 섹션에서 1tie-Mix² 모형은 8개의 빌딩 블록

을 이용하는데, 이들은 가우스가 되도록 선택되었다. 이들 가우스의 평균은

의 위치에서 시작하고

의 단계에서 진행하는 전사체의 길이에 걸쳐 균등하게 분포되는데, 여기서 양쪽 숫자는 가장 가까운 정수로 반올림되었다. 따라서 2702 bp의 길이를 전사체 2의 경우에, 예로서, 이것은 169, 507, 845, 1183, 1521, 1859, 2197, 2535에서 평균을 유발한다. 각 가우스의 표준 편차는 첫 번째 평균에 동등하도록 세팅되었다. 따라서 전사체 2의 경우에 각 가우스의 표준 편차는 169이었다. 상기 모형에서 초기 β _j 는 1/8로 세팅되었다. 가중된 블록

와 함께 전사체 2에 대한 p(r|t=i)의 결과적인 초기화는 도면 15에서 묘사된다. 첫 번째와 마지막 블록은 전사체 내에서 것들보다 약간 높다. 이것은 이들 가우스가 그들의 꼬리의 부분을 결여하고, 그리고 그들의 정규화 상수가 이런 이유로, 다른 가우스의 것보다 높다는 사실에 기인한다.

상기 모형에서 초기 존재비는 α₁=α₂=α₃=1/3에 세팅되었다. EM 알고리즘은 차후 반복 사이에 α_i와 β_j 사이에 차이가 0.001 미만이거나 또는 차후 반복 사이에 로그 우도에서 증가가 0.5 미만일 때까지 적용되었다. 이들 종결 조건에서 EM 알고리즘은 20회 반복 후 수렴하였다.

도면 16은 α_i의 수렴을 보여준다. 이러한 도면의 x 축은 EM 알고리즘의 반복을 나타내고, 반면 y 축은 α_i의 값을 나타낸다. 도면 16에서 3개 곡선 모두 1/3의 초기 값으로부터 기원하고 정확한 해답에 가까운 값에 수렴하는데, 이들은 도면 16에서 3개의 수평선에 의해 표시된다. 최종 값은 α₁=0.29, α₂=0.3, α₃=0.41이었다. 대조적으로, Cufflinks 모형으로 도출된 추정값은 α₁=0.37, α₂=0.12, α₃=0.51이다. 이들 숫자는 α_i의 참값과 추정된 값 사이에 쿨백 라이블러 (KL) 발산과 함께 표 2에서 요약된다. KL 발산은 동일한 분포에 대해 0이고, 그리고 분포 사이에 차이가 증가함에 따라서 증가한다. 표 2는 Cufflinks 모형에 대한 참 α_i와 추정된 α_i 사이에 KL 발산이 1tie-Mix² 모형의 KL 발산을 2 크기 자릿수로 초과한다는 것을 보여준다.

도면 17은 8개 블록에 대한 8개 베타의 수렴을 보여준다. 앞서와 같이, x 축은 EM 알고리즘의 반복을 나타내고, 그리고 y 축은 β_j의 값을 나타낸다. β_j 모두 1/8로 초기화되었기 때문에, 이들 곡선은 동일한 위치로부터 기원한다. α_i의 경우에서처럼 β_j에 대한 그라운드 트루스가 없기 때문에, 이들의 품질은 단지, 결과의 분포 p(r|t=i)를 점검함으로써 평가될 수 있다. 최종 β_j는 도면 18에서 도시되는데, 여기서 이들은 그들의 블록 id에 대하여 플롯팅된다. 이것은 β_j가 p(r|t=i)의 모양을 거의 추종한다는 것을 보여주고, 이것은 차례로, 도면 19로부터 목격될 수 있는 바와 같이, 정확한 p(r|t=i)의 우수한 근사를 야기한다.

최종적으로, 도면 20은 전체 확률 분포 p_total(r)의 수렴을 보여준다. 이것은 전체 우도의 면에서 EM 알고리즘이 첫 번째 반복에서 정확한 해답의 방향에서 가장 큰 단계를 차지한다는 것을 지시한다. 차후 반복은 전체 우도의 증가에 대한 상대적으로 중간 효과를 갖는다. 하지만, 초기 반복에서 α_i는 여전히, 그들의 정확한 해답으로부터 거리가 멀다. 이것은 전체 우도에서 중간 변화 역시 파라미터의 의미있는 변화를 야기할 수 있다는 것을 지시한다. 결과적으로, EM 알고리즘의 종결 조건은 전체 우도의 작은 증가에 대하여 관용적이어야 한다.

4.1.2 가중치의 60개 세트에서 1tie-Mix²와 Cufflinks 모형 사이에 비교

본 섹션은 가중치 α₁,α₂,α₃의 60개 상이한 세트로 실험을 논의하는데, 이들은 섹션 4의 시작 부분에서 설명된 절차에 따라 선택되었다. 따라서 α₃은 값 0.2, 0.4, 0.6, 0.8을 취하고, 그리고 α₁과 α₂는 0과 1-α₃ 사이에 등거리 간격에서 분포된다. 가중치의 각 세트에 대해 p(r|t=i)은 도면 15에서처럼 전사체 2에 대해 초기화되었다. 다른 2개 전사체에 대한 분포 p(r|t=i)는 그에 맞춰 초기화되었다. EM 알고리즘은 수렴 규준 중에서 한 가지가 부합될 때까지 실행되었다. 최종 반복에서 획득된 α_i가 1tie-Mix² 모형의 결과로서 선택되었다. 유사하게 EM 알고리즘은 수렴 규준 중에서 한 가지가 부합될 때까지 Cufflinks 모형에 대해 수행되었고, 그리고 최종 반복으로부터 α_i가 Cufflinks 모형의 결과로서 선택되었다. 도면 21은 α₃=0.2에 대한 이들 실험의 결과를 보여준다. 이러한 그래프의 x 축은 α_i의 참값을 제공하고, 반면 y 축은 참 α_i 및 Cufflinks와 1tie-Mix² 모형에 의해 추정된 α_i를 제공한다. 도면 21에서 쇄선은 참 α_i를 지칭하고, 반면 점선과 파선은 각각, Cufflinks와 1tie-Mix² 추정값을 지칭한다. 도면 21은 α₃=0.2의 경우에 1tie-Mix² 추정값은 α_i의 참값에 매우 잘 정합하고, 반면 Cufflinks 추정값은 상당히 불량하다는 것을 보여준다. α₂=0의 경우에만 Cufflinks와 1tie-Mix² 모형의 추정값이 일치한다. 이러한 상황에서 도면 7로부터 목격될 수 있는 바와 같이, 완전하게 별개인 전사체 1과 3만 존재한다. 이러한 경우에 이런 이유로, EM 알고리즘은 전사체 1과 3에 배정되는 단편의 숫자를 단순히 계수하고, 그리고 p(r|t=i)의 모양과는 관계가 없다. 도면 22는 α₃=0.4에 대한 Cufflinks와 1tie-Mix² 모형의 추정값을 보여준다. α₃=0.2에서처럼, 1tie-Mix² 모형의 추정값은 Cufflinks 모형과 비교하여 매우 정확하다. α₂=0의 경우에만 양쪽 추정값이 일치한다. α₃=0.6 및 α₃=0.8에 대해 유사한 모습이 나타난다.

표 3으로부터 목격될 수 있는 바와 같이, 이것은 전반적으로, Cufflinks 모형의 경우에 0.12368의 평균 KL 발산 및 1tie-Mix² 모형의 경우에 3.6369e-04의 평균 KL 발산을 야기한다. 따라서, 평균 KL 발산의 면에서 1tie-Mix² 모형의 정확도는 Cufflinks 모형의 정확도를 3 크기 자릿수로 초과한다.

4.2 5tie-Mix² 모형

본 섹션에서 모형은 전사체 특정한 오프셋 및 환산 파라미터 ν_i와 λ_i를 추정하고, 그리고 이런 이유로, 정확한 전사체 주해에 의존하지 않는다. 따라서, 이러한 모형의 잠재력을 증명하기 위해, 이것은 3개의 부정확한 전사체 주해로 훈련된다. 이들 주해는 도면 25에서 목격될 수 있다. 이러한 도면에서 실선은 도면 7 및 표 1에서와 동일한 정확한 전사체 주해를 지시한다. 도면 25에서 점선은 본 섹션에서 모형을 훈련하는데 이용되었던 부정확한 전사체 주해를 보여준다. 부정확한 전사체 주해에서 엑손의 정확한 시작과 끝 위치 및 정확한 주해로부터 이들의 차이는 표 4에서 제공된다. 본 섹션에서 5tie-Mix² 모형은 확장된 전사체 주해를 이용하는데, 이들은 위치 1에서 시작하고 위치 10000에서 끝난다. 따라서 전사체 2와 3의 확장된 주해는 동일하다. 5tie-Mix² 모형의 ν_i와 λ_i는 부정확한 주해에 합치하도록 초기화되고 표 5에서 제공된다. 전사체 2의 λ_i는 1이 되도록 선택되었고, 그리고 전사체 1과 3의 λ_i는 이런 이유로, 1275/2900=0.4397 및 2300/2900=0.7931에 의해 제공된다. 5tie-Mix² 모형의 파라미터의 복합적인 상호작용으로 인해, 이러한 모형의 우도 표면은 준최적 국부 최대값을 갖는다. 이들 최대값 중에서 한 가지에서 포획되는 것을 방지하고, 그리고 이런 이유로, α_i에 대한 준최적 추정을 획득하기 위해, 모형 파라미터를 적절하게 초기화하는 것이 필요하다. 여기에서 전략은 부정확한 주해에 대한 실시예 4.1.1에서 1tie-Mix² 모형에 대해 획득된 초기 값을 이용하고 이들을 일정한 양으로 무작위로 변화시키는 것이었다. 이러한 방식으로 200개의 상이한 초기 파라미터 세트가 산출되었는데, 여기에서 EM 알고리즘은 수렴 규준 중에서 한 가지가 부합될 때까지, 다시 말하면, 로그 우도에서 차이가 0.5 미만이거나 또는 α_i와 β_j 사이에 차이가 차후 반복 사이에 0.001 미만일 때까지 수행되었다. 200개의 획득된 결과로부터 가장 높은 우도를 갖는 것들이 추정값으로서 선택되었다.

부정확한 주해를 교정할 수 없는 본 섹션에서 연구된 다른 모형, 다시 말하면, Cufflinks와 1tie-Mix² 모형의 경우에, 추정된 α₂가 너무 높다는 것을 예상할 수 있는데, 그 이유는 전사체 2가 이의 이동된 시작으로 인해 상실하는 것보다 이의 증가된 길이에 의해 더욱 많은 단편을 획득하기 때문이다. 유사하게, α₁은 α₃보다 더욱 강하게 과소평가되는 것으로 예상될 수 있는데, 그 이유는 전사체 3이 이의 시작에서 상실하는 것보다 전사체 1이 이의 끝에서 더욱 많은 단편을 상실하기 때문이다.

부정확한 전사체 주해

	엑손 1	엑손 2	길이
전사체 1	[900,2174]	해당 없음	1275
전사체 2	[1600,3200]	[4000,5298]	2900
전사체 3	[2800,3200]	[4000,5898]	2300

정확한 대 부정확한 전사체 주해. 상대적 차이 및 초기 ν_i와 λ_i

	시작 차이	길이 차이	초기 νi	초기 λi
전사체 1	-100	-226	900	0.4397
전사체 2	+100	+198	1600	1
전사체 3	+198	+298	2800	0.7931

4.2.1 α₁=0.28, α₂=0.32, α₃=0.4에 대한 5tie-Mix² 모형의 수렴

도면 26은 가중치 α₁=0.28, α₂=0.32, α₃=0.4에 대한 EM 알고리즘의 수렴 후 가장 높은 우도를 갖는 초기 파라미터 세트에 대한 α_i의 수렴을 보여준다. 이러한 초기 파라미터 세트에 대해 EM 알고리즘은 149회 반복 후 수렴하였다. 도면 26에서 x 축은 EM 알고리즘의 반복을 제공하고, 그리고 y 축은 상응하는 α를 제공한다. 도면 26에서 점선은 EM 알고리즘의 코스 동안 α_i를 보여주고, 반면 수평 쇄선은 참 α_i를 보여준다. 도면 26은 α_i이 참값에 매우 가까운 값에 수렴한다는 것을 지시한다. 이것은 또한, 표 6에서 반영되는데, 이것은 EM 알고리즘의 최종 반복 후 정확한 가중치로부터 그들의 KL 발산과 함께 추정된 α_i를 제공한다. EM 알고리즘 동안 다른 모형 파라미터의 값은 도면 27, 28, 29, 30과 31에서 도시된다. 앞서와 같이, β_j, μ_j, σ_j의 획득된 값의 품질은 단지, 결과의 p(r|t=i)를 조사함으로써 평가될 수 있다. 대조적으로, 전사체 특정한 이동과 환산 파라미터 ν_i와 λ_i는 최종 모형이 부정확한 초기 가정을 교정하는 지에 관한 지시를 제공한다. ν_i에서 증가는 오른쪽으로 이동을 지시하고, 반면 λ_i에서 증가는 길이에서 증가를 지시한다. 도면 30과 31은 이런 이유로, 5tie-Mix² 모형이 전사체 1의 시작을 오른쪽으로 이동시키고 이의 길이를 증가시키고, 반면 전사체 2와 3의 시작을 왼쪽으로 이동시키고 이들의 길이를 줄인다는 것을 지시한다. 이것은 표 5에 나타나 있는 바와 같이, 정확한 주해로부터 부정확한 주해의 일탈과 일치한다.

EM 알고리즘 동안 부정확한 전사체 주해의 교정은 또한, 도면 32, 33과 34에서 목격될 수 있는데, 이들은

의 수렴을 보여준다. 이들 도면에서 수직선은 문제되는 전사체의 전사체 좌표에서 정확한 시작과 끝을 지시한다. 도면 32는 전사체 1에 대한

의 수렴을 보여준다. 이러한 도면에서 쇄선은 초기 전사체 특이적

를 제공하고, 반면 실선은 149회 반복 후 EM 알고리즘의 결과인

를 보여준다. 도면 32는 도면 30과 31에 따라,

가 EM 알고리즘 동안 약간 오른쪽으로 움직이고 길이에서 증가한다는 것을 보여준다.

에서 전사체 시작과 전사체 길이의 결과의 값은 거의 완벽하다. 도면 32는 또한, 정확한 주해를 향한 주요 단계가 EM 알고리즘의 첫 몇 회 반복에서 발생한다는 것을 보여준다. 이것은 도면 30과 31과 일치하는데, 이들은

과

이 20회 EM 반복 후 실질적으로 일정하게 머물러 있다는 것을 보여준다. 도면 33과 34는 전사체 2와 3에 대한

의 수렴을 보여준다. 이들 도면은 전사체 1에서처럼, EM 알고리즘이 양쪽 전사체를 왼쪽으로 이동시키고 이들의 길이를 줄임으로써 초기 부정확한 전사체 주해를 교정한다는 것을 보여준다. 최종적으로, 도면 35는

의 수렴을 보여준다. 도면 35에서 실선에 의해 목격될 수 있는 바와 같이, 최종

은 도면 13에서 히스토그램의 모양에 우수한 근사이다.

4.2.2 가중치의 60개 세트에서 5tie-Mix², 1-tie-Mix²와 Cufflinks 모형 사이에 비교

본 섹션은 가중치 α₁, α₂, α₃의 60개 상이한 세트 및 5tie-Mix², 1tie-Mix²와 Cufflinks 모형으로 실험을 논의한다. 도면 25에서 부정확한 전사체 주해가 1tie-Mix²와 Cufflinks 모형에 대한 참고의 고정된 프레임으로서 5tie-Mix² 모형의 파라미터를 초기화하는데 이용되었다. 도면 36은 α₃=0.2의 참값에 대한 5tie-Mix²와 Cufflinks 모형에 대해 추정된 α_i를 보여준다. 이러한 도면은 참 α_i 및 5tie-Mix² 모형에 의해 추정된 α_i 사이에 우수한 합치를 보여주고, 반면 Cufflinks 모형에 의해 추정된 α_i은 정확한 해답으로부터 강하게 일탈한다. 참 α₁의 작은 값의 경우에만 Cufflinks 모형은 5tie-Mix² 모형에 근접하는 추정값을 산출한다. 이러한 경우에, 전사체 1은 매우 낮은 농도를 갖고, 그리고 이런 이유로, 전사체 2와 3으로부터 단편만 EM 알고리즘에서 일정한 역할을 수행한다. 이들 단편은 EM 알고리즘에 의해 전사체 2와 3 사이에 거의 정확하게 분포되는데, 그 이유는 전사체 1이 부정확한 끝으로 인해 전사체 2에 상실하는 단편의 숫자가 미미하기 때문이다. 하지만, α₁이 증가함에 따라서, Cufflinks 모형은 α₂를 심각하게 과대평가하고, 그리고 α₁과 α₃의 과소평가가 더욱 확연해진다. 이러한 효과는 도면 37에서 훨씬 강한데, 여기서 α₂에 대한 Cufflinks 추정값은 짧은 초기 감퇴 후, 0.5에서 거의 일정한 것으로 보인다. 대조적으로, 5tie-Mix² 모형에 의해 추정된 α_i는 이번에도, 가중치의 전체 범위에 걸쳐 매우 정확하다.

도면 38과 39에서 5tie-Mix² 모형의 추정값은 다시 한 번 매우 정확하고, 반면 Cufflinks 모형은 α₁과 α₃을 심각하게 과소평가한다. 이에 더하여, 도면 38과 39에서 Cufflinks 모형은 0.4와 0.28의 알파2에 대해 각각 거의 일정한 추정을 산출한다. 표 7로부터 목격될 수 있는 바와 같이, 완전한 60개 파라미터 세트에 걸쳐 Cufflinks 모형은 0.21977의 평균 KL 발산을 획득하고, 반면 5tie-Mix² 모형은 0.014482의 평균 KL 발산을 획득한다. 따라서 평균 KL 발산의 면에서 5tie-Mix² 모형의 정확도는 Cufflinks 모형의 정확도를 15의 인자로 초과한다. 실제로, 5tie-Mix² 모형은 정확한 주해에서 Cufflinks 모형보다 훨씬 정확한데, 상기 Cufflinks 모형은 표 3에 따르면, 0.12368의 평균 KL 발산을 산출한다. 따라서 평균 KL 발산의 면에서 부정확한 전사체 주해에서 5tie-Mix² 모형의 정확도는 정확한 주해에서 Cufflinks 모형의 정확도를 8의 인자로 초과한다.

4.2.3 1par-Mix² 모형과의 비교

1tie-Mix² 모형이 섹션 4.1.2에서 정확한 전사체 주해로 실험에서 Cufflinks 모형보다 훨씬 우수하였기 때문에, 5tie-Mix² 모형 및 1tie-Mix² 모형 둘 모두를 부정확한 전사체 주해에서 비교하는 것은 가치가 있다. 도면 40은 1tie-Mix²와 5tie-Mix² 모형에 대해 추정된 α_i 및 0.2의 참 α₃의 값을 보여준다. 도면 36에서처럼, α₂는 고정된 부정확한 전사체 주해를 이용하는 1tie-Mix² 모형에 의해 심각하게 과소평가된다. 비록 도면 40에서 1tie-Mix² 모형의 α_i에 대한 추정값이 도면 36에서 Cufflinks 모형의 추정값보다 약간 더 정확하긴 하지만, 이들은 여전히, 5tie-Mix² 모형보다 훨씬 나쁘다. 이러한 경향은 또한, 도면 41, 42와 43에서 도시되는데, 이들은 1tie-Mix² 모형에 대한 추정값이 Cufflinks 모형에 대한 추정값에 비하여 약간 향상되긴 하지만, 5tie-Mix² 모형의 추정값보다 훨씬 나쁘다는 것을 보여준다. 전반적으로, 이것은 표 8에서 1tie-Mix² 모형에 대한 평균 KL 발산을 야기하는데, 이것은 KL 발산의 면에서 5tie-Mix² 모형의 정확도가 1tie-Mix² 모형의 정확도보다 약 12 배 높다는 것을 보여준다. 요약하면, 본 섹션에서 실험은 정확한 전사체 주해의 부재에서 α_i에 대한 신뢰할 만한 추정을 획득하기 위해, 정확한 전사체 위치를 학습할 수 있는 모형을 이용하는 것이 중요하다는 것을 보여준다.

실시예 5: 결론

마지막 섹션에서 실험은 1tie-Mix²와 5tie-Mix² 모형을 정확하고 부정확한 전사체 주해 둘 모두에서 Cufflinks 모형과 비교하였다. 60개 파라미터 세트 α₁, α₂, α₃의 세트에서 실험의 결과는 표 9에서 요약된다. 이들 결과는 정확한 전사체 주해에서 1tie-Mix² 모형이 Cufflinks 모형보다 훨씬 뛰어나고, 반면 부정확한 전사체 주해에서 5tie-Mix² 모형이 1tie-Mix²와 Cufflinks 모형 둘 모두를 능가한다는 것을 보여준다. 이것은 적절한 Mix² 모형의 이용이 존재비 추정값의 정확도를 상당히 향상시킨다는 것을 암시한다.

실시예 6: 선행 기술과의 비교

Cufflinks는 전사체 어셈블리와 전사체 존재비 추정을 위한 방법을 실행하는 프로그램이다. 이러한 실행의 상세한 설명은 Trapnell et al., 2010에서 발견될 수 있다. Cufflinks에서, 본 발명에 따른 α_i에 상응하는 전사체 존재비는 방정식 (13)에서 인수분해를 통해 추정되는데, 여기서 p_trans(s(r)|t=i,l(r))은 균일하고 p_trans(l(r)|t=i)는 전사체 t=i와 관계가 없다. Mix² 모형과 대조적으로, Cufflinks는 데이터로부터 단편 시작 부위 s(r)의 분포를 학습하지 못하고, 그리고 게다가, 정확한 전사체 주해의 이용가능성에 의존한다.

Cufflinks 존재비 추정에 확장을 실행하는 Roberts et al., 2011에서, p_trans(r|t=i)는 방정식 (13)에서처럼 인수분해되고, 그리고 p_trans(s(r)|t=i,l(r))는 전사체 좌표에서 각 위치에 대해 규정되는 서열 특정한 위치 가중치의 정규화된 결과물이다. 이것은 심지어 단일 전사체에 대해서도 훈련될 필요가 있는 다수의 파라미터를 유발한다. Roberts et al., 2011의 실험에서 모형을 연산적으로 다루기 쉽게 만들기 위해, 위치 가중치는 단계 함수에 한정되어야 하고, 이들 중에서 단지 5개만 모든 전사체 사이에 추정되고 공유된다. Roberts et al., 2011에서 모형과 대조적으로, 본 발명의 모형은 모든 확률 p_trans(s(r)|t=i)을 개별적으로 추정하지 않고, 혼합체의 파라미터 (16)를 추정한다. 결과적으로, 훨씬 적은 파라미터가 본 발명의 모형에 대해 추정될 필요가 있고, 그리고 이의 파라미터 추정은 이런 이유로, Roberts et al., 2011에서 모형의 것보다 더욱 견실하고 연산적으로 다루기 쉽다. 이에 더하여, (Roberts et al., 2011)에서 모형은 전사체 존재비에 대한 정확한 시작 추정값을 필요로 하고, 그리고 이의 가중치는 이런 이유로, 단일 동종형 유전자에서 추정된다. 따라서 Roberts et al., 2011에서 모형은 본원에서 도시된 실험에서 이용된 데이터에 적용가능하지 않다. 본 발명의 모형은 다른 한편, 전사체 존재비에 대한 정확한 시작 값을 필요로 하지 않고, 그리고 이런 이유로, 다중-동종형 유전자에서 훈련될 수 있다. 이에 더하여, Roberts et al., 2011에서 모형은 본 발명의 모형과 달리, 정확한 전사체 주해를 필요로 한다.

Wu et al., 2011에서 단편의 분포에서 편향은 단일 가중치를 각 엑손과 전사체에 배정함으로써, 유전자의 엑손에 대해 모형화된다. 본 발명의 모형과 대조적으로, 이것은 p_trans(s(r)|t=i)가 각 엑손에서 일정하도록 제한한다. 이에 더하여, Wu et al., 2011에서 가중치는 확률적 프레임워크 외부에서 휴리스틱 방법에 의해 부분적으로 추정된다. Wu et al., 2011에서 작업은 이것은 엑손 상에서 리드 수치의 확률에 대한 포아송 분포를 이용하고, 반면 본 발명의 모형은 전사체를 고려할 때 단편의 확률에 대해 혼합체 모형을 이용한다는 점에서 본 발명의 모형과 더욱 상이하다. 이에 더하여, Wu et al., 2011에서 작업은 정확한 전사체 주해에 의존한다.

Li et al., 2010에서 작업은 서열 특정한 편향을 교정하고, 그리고 이런 이유로, Roberts et al., 2011에서 서열 특정한 가중치와 유사하다. 이것은 또한, 수치 데이터를 모형화하기 위한 포아송 분포의 용법에서 Wu et al., 2011에서 작업과 유사하다. Wu et al., 2011과 대조적으로, Li et al., 2010에서 수치 데이터는 엑손 기초되지 않지만 커버리지에 의해 제공된다. Li et al., 2010에서 작업은 이것이 전사체 존재비를 추정하지 않는다는 점에서, 본 발명의 모형, 그리고 Roberts et al., 2011과 Wu et al., 2011에서 모형과 상이하다.

Glaus et al., 2012에서 베이스 모형이 전사체에 의해 산출된 리드를 관찰하는 확률에 대해 이용된다. 이러한 모형은 Roberts et al., 2011에서와 동일한 서열 특정한 위치 가중치, 그리고 전사체 존재비를 추정하기 위한 Gibbs 표본추출 절차를 이용한다. 따라서, 본 발명의 모형과 비교하여, Glaus et al., 2012에서 파라미터 추정은 최고 우도 프레임워크보다는 베이스 프레임워크 내에서 수행된다. Glaus et al., 2012에서 모형은 또한, 이것이 p(r|t=i)에 대해 임의의 혼합체를 이용하지 않고 정확한 전사체 주해를 필요로 한다는 점에서, 본 발명의 모형과 상이하다.

Li et al., 2010 & 2011에서 모형은 Glaus et al., 2012에서 모형과 유사한 구조를 갖는 베이스 모형이다. 하지만, Glaus et al., 2012와 대조적으로, Li et al., 2010 & 2011에서 전사체 존재비는 특히, EM 알고리즘을 갖는 최고 우도 프레임워크 내에서 훈련된다. 나머지 모형 파라미터, 예를 들면, p(s(r)|t=i)는 사전에 도출되거나 또는 휴리스틱 방법으로 추정된다. 이것은 모든 파라미터를 최고 우도 프레임워크 내에서 추정하는 본 발명의 모형과 대조적이다. 이에 더하여, Li et al., 2010 & 2011에서 모형은 p_trans(r|t=i)을 모형화하기 위해 혼합체를 이용하지 않고, 그리고 정확한 전사체 주해를 필요로 한다.

본 발명에 따라서, 기대값 최대화 알고리즘이 본 발명의 모형을 추정하는데 이용되었다. 기대값 최대화 알고리즘의 일반적인 프레임워크가 Dempster et al., 1977에서 개발되긴 했지만, 본 발명의 모형에 대한 구체적인 EM 갱신 공식은 본원에서 설명된 바와 같이 도출되었다.

실시예 7: Mix ² 모형으로 전사체 군 묶기

방정식 (20) 다음의 단락에서, Mix² 모형의 파라미터의 묶기는 상이한 전사체

에 대한

사이에 일정한 유사성을 암시하는 것으로 언급되었다. 이런 이유로, Mix² 모형에 대한 파라미터는 이러한 유사성을 표시하는, 전사체

사이에서만 묶여져야 한다. 전사체가 단지 군 내에서 전사체만 그들의 파라미터를 공유하도록 상이한 군으로 분리되면, Mix² 모형의 EM 갱신 공식 (27), (36)과 (37)은 변형되어야 한다. 아래에서, 각 전사체

는 연관된 군

를 갖는데, 이것은 함수

를 통해 검색될 수 있다. 이후, 군

내에서 파라미터

와

의 EM 갱신 공식은 아래와 같이 제공된다.

Lagrange 방법을 이용하여 제약을 강제하고:

(51),

그리고

에 대하여 도함수를 취하는 것은 다음을 야기하고:

(52),

이것은 일부 재배열 후, 다음을 결과한다:

(53),

여기서

(54)

전사체 사이에 묶이는 나머지 파라미터, 다시 말하면,

와

의 경우에, 갱신 공식 (36)과 (37)은 전사체

의 완전한 세트에서 합계를 군

내에서 전사체에서 합계로 대체함으로써 변형되어야 한다, 다시 말하면,

(55)

(56)

유전체 좌위 내에 각 전사체를 이의 자체 군으로 배치하는 것은 모든 파라미터 묶기를 제거하고, 그리고 앞서 언급된 바와 같이, 단편 분포와 전사체 존재비에 대한 부정확한 추정을 야기한다. 전사체를 상이한 군에 배치할 때, 이런 이유로, 상이한 단편 분포

의 분리 사이에 올바른 균형을 잡고, 그리고 Mix² 모형의 안정성을 담보하기 위해 충분한 숫자의 파라미터의 묶기를 유지하는 것이 중요하다. 합리적인 요건은 예로서, 각 군이 최소한 2개의 전사체를 내포하는 것이다. 대안으로, 3개 또는 그 이상 군의 존재에서 기껏해야 하나의 군이 단일 전사체를 내포하도록 하는 것이 합리적이다.

실시예 8: Mix ² 모형으로 단편 길이 분포의 훈련

섹션 4에서 논의된 실험에서 전사체 특정한 단편 확률

은 (48)에서처럼 인수분해되었고, 그리고 단편 길이의 확률 분포

는 단지, 전사체

의 단편 시작

과 길이

에만 의존하는 것으로 가정되었다. 이에 더하여,

는 정해진 것으로 가정되었다. 이들 가정은

가 데이터 세트

로부터 추정되면 불필요한데, 이것은 Mix² 모형의 프레임워크 내에서 행위될 수 있다. 이런 목적으로, 단편 시작

의 분포에 대해,

는 혼합체로서 기재된다, 다시 말하면,

(57)

여기서

는 앞서

에 의해 표시된 숨겨진 변수이다. 여기에서

는 "단편 시작의 빌딩 블록"에 대한 연상기호이고, 반면 숨겨진 변수

은 "단편 길이의 빌딩 블록"에 대한 연상기호이다.

가

와 관계가 없고, 그리고 따라서 (57)이

(58)로 환산된다고 가정하는 것은 합리적이다.

(15), (48)과 (58)을 결합하는 것은 단편의 확률에 대해 하기 표현을 산출한다:

(59)

여기서 다음이 유효하다:

(60)

따라서 (59)는

의 혼합체 가중치가

와

의 결과물이고, 그리고

의 혼합체 성분이 (60)에서 조건부 확률 분포의 결과물인 Mix² 모형이다.

(15)와 (25)에 더하여, 다음이 (59)로부터 도출될 수 있다:

(61)

그리고 결과로서, (22)와 (27)과 유사하게,

는 아래와 같이 EM 알고리즘으로 추정될 수 있다:

(62)

여기서

(63)

그리고

(64)

단편 길이 빌딩 블록

의 경우에, 예로서, 평균이

과

또는

과

사이에 등거리로 분포되고, 그리고 구별된 또는 연속적 1-차원 확률 공간 내에서 정규화되는 가우스가 이용될 수 있다. 연속적 1-차원 확률 공간에서 분포가 선택되면, 이들의 내부 파라미터, 예를 들면, 평균, 표준 편차, 이동과 환산 파라미터는 단편 시작

의 확률 분포에 대해 이들 파라미터와 유사하게 추정될 수 있다.

실시예 9: 고정된 전사체 종결점으로

추정

5tie-Mix² 모형은 이동과 환산 파라미터

와

를 추정한다.

가 0에 세팅되고 EM 알고리즘 동안 갱신되지 않으면, 전사체

의 시작은 변화가 없다. 전사체의 끝점이 고정되면, 5tie-Mix² 모형은 약간 변형되어야 한다.

이 주해에 따른 전사체

의 길이이면, 전사체

의 고정 끝점에 대해 다음이 유효하다:

(65)

이것은

에 대해 하기 EM 갱신 공식을 산출한다:

(66)

이것으로부터

가 (65)에 따라 도출될 수 있다, 다시 말하면,

(67)

실시예 10: Mix² 모형을 다른 편향 모형과 결합

섹션 4에서 실험에서 논의된 Mix² 파라미터 묶기는 위치 단편화 편향, 다시 말하면, 전사체 내에서 단편 시작에 관련된 편향에 대한 모형을 실행한다. 다른 종류의 편향, 예를 들면, 서열 특정한 편향은 다른 모형, 예를 들면, [6]에서 변수 길이 숨겨진 마르코프 모형 (VLMM)으로 설명될지도 모른다. 전형적으로, 비-위치 편향에 대한 모형은 뉴클레오티드 서열의 관찰된 빈도를 편향 없는 데이터의 귀무 가설 하에 그들의 빈도에 비교한다. 이러한 비교는 단편

의 단일 사본이 편향 있는 데이터에서 관찰된다는 점을 고려하면, 임의의 편향의 부재에서 단편

의 다중성

위에서 확률 분포

를 도출하는데 이용될 수 있다. 분포

는 이후, Mix² 모형의 EM 갱신 공식에서 각 단편

을

의 예상된 다중성으로 가중함으로써, 비-위치 편향을 데이터로부터 연산적으로 제거하는데 이용될 수 있다. 존재비

의 EM 갱신 공식의 경우에, 예로서, 이것은 다음을 야기한다:

(68)

분포

는 Mix² 모형 하에 데이터 세트

의 우도를 최대화함으로써 추정될 수 없는 것으로 유의되어야 하는데, 그 이유는 이것이

의 무한한 기대값을 야기할 것이기 때문이다. 본 섹션에서 논의된 바와 같이, Mix² 모형과 다른 편향 모형의 조합의 경우에, 이런 이유로 Mix² 모형의 최고 우도 프레임워크 외부에서

를 추정하는 것이 중요하다.

실시예 11: FPKM과 RPKM 값을 교정하기 위해 Mix ² 모형을 이용

RNA-Seq에서 전사체의 농도는 통상적으로, FPKM (백만 리드당 킬로 염기쌍마다 단편) 또는 RPKM (백만 리드당 킬로 염기쌍마다 리드) 척도로 계측되는데, 여기서 전사체

에 대해 후자는 다음에 의해 제공된다:

(69)

FPKM 척도의 경우에 (69)에서 전사체

의 길이

는 조정된 전사체 길이

에 의해 대체된다 [7]. 따라서, 정확한 FPKM과 RPKM 값을 획득하기 위해, (69)에서 정확한 전사체 길이

를 이용하는 것이 중요하다. 5tie-Mix² 모형은 부정확한 전사체 주해를 교정할 수 있고, 그리고 이런 이유로, 더욱 높은 정밀도의 전사체 길이를 산출하는데, 이것은 RPKM과 FPKM 값의 향상된 정확도를 유발한다. 상이한 방법이 5tie-Mix² 모형에 기초하여 전사체 길이에 대한 추정값을 획득하는데 이용될 수 있다. 전사체 시작은 예로서, 위치

로서 추정될 수 있는데, 여기서

(70)이고,

이때

이다. 전형적으로,

는 작은 양성 값, 예를 들면,

일 것이다. 유사하게, 전사체 끝은 (70)이 1에 가까울,

의 값

, 예를 들면,

에 의해 대체된

에서 유지되는 위치로서 추정될 수 있다. 유사하게,

와

를 갖는 (70)이 임의의 혼합체 성분

에 적용될 수 있다. 이러한 경우에, 전사체 시작을 결정하는데 이용된 혼합체 성분은 전형적으로, 작은

의 영역에서 집중될 것이고, 반면 전사체 끝을 결정하는데 이용된 혼합체 성분은 큰

의 영역에서 집중될 것이다. 5tie-Mix² 모형으로부터 전사체

의 길이에 대한 추정값을 획득하는 다른 방법은 전사체 주해의 길이를 Mix² 모형의 수렴 후

로 곱하는 것이다.

참고문헌

[1] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society, Series B, 39(1) : 1 38, 1977.

[2] Peter Glaus, Antti Honkela, and Magnus Rattray. Identifying differentially expressed transcripts from RNA-scq data with biological variation. Bioinformatics, 28(13) :1721- 1728, 2012.

[3] Bo Li and Colin Dewey. RSEM: accurate transcript quantification from rna-seq data with or without a reference genome. BMC Bioinformatics, 12(1) :323, 2011.

[4] Bo Li, Victor Ruotti, Ron M Stewart, James A Thomson, and Colin N Dewey. Rna-seq gene expression estimation with read mapping uncertainty. Bioinformatics, 26(4) :493- 500, Feb 2010.

[5] Jun Li, Hui Jiang, and Wing Wong. Modeling non-uniformity in short-read rates in RNA-seq data. Genome Biology, 11 (5) :R50-, 2010.

[6] Adam Roberts, Cole Trapnell, Julie Donaghey, John L Rinn, and Lior Pachter. Improving rna-seq expression estimates by correcting for fragment bias. Genome Biol, 12(3) :R22, Mar 2011.

[7] Cole Trapnell, Brian A Williams, Geo Pertea, Ali Mortazavi, Gordon Kwan, Marijke J van Baren, Steven L Salzberg, Barbara J Wold, and Lior Pachter. Transcript assembly and quantification by RNA-Seq reveals unannotatcd transcripts and isoform switching during cell differentiation. Nat Biotechnol, 28(5) :511 515, May 2010.

[8] Zhengpeng Wu, Xi Wang, and Xuegong Zhang. Using non-uniform read distribution models to improve isoform expression inference in RNA-Scq. Bioinformatics, 27(4) :502- 508, Feb 2011.

[9] Roberts et al: Genome Biology, 12 (3) (2011): R22

[10] Wen-Ping et al: Genome Biology, 8 (6) (2007): R98

Claims (26)

1. 다음 단계를 포함하는, 전사체 존재비를 추정하는 방법:
   a) 관심되는 유전자 좌위의 전사체의 잠재적 혼합체로부터 전사체 단편 염기서열결정 데이터를 획득하는 단계,
   b) 상기 단편 염기서열결정 데이터를 관심되는 상기 좌위의 유전자 좌표에 배정하여, 단편 유전자 좌표 커버리지의 데이터 세트를 획득하되, 각 유전자 좌표에 대한 상기 커버리지가 통합되어 커버리지 포락선을 형성하는 단계,
   c) 상기 혼합체의 전사체의 숫자를 세팅하는 단계,
   d) 각 전사체 i에 대한 모형화된 유전자 커버리지의 확률 분포 함수를 미리 세팅하되, i가 전사체에 대한 수치 식별자를 나타내고, 상기 확률 분포 함수가 최소한 2개의 확률 부함수 j의 합계와 곱셈된, 상기 전사체 i의 가중 인자 α_i에 의해 규정되며, j가 확률 부함수에 대한 수치 식별자를 나타내고, 각 확률 부함수 j가 가중 인자 β_i,j에 의해 독립적으로 가중되며, 여기서 확률 부함수 j는 비주기 함수이고, 그리고 여기서 전사체에 대한 확률 부함수 j는 양성 값을 갖는 전사체의 전체 길이를 커버하도록 유전자 좌표에서 배치되거나 또는 이동되는 단계,
   e) 각 전사체의 확률 분포 함수를 더하여 합계 함수를 획득하는 단계,
   f) 합계 함수를 커버리지 포락선에 적합시키고, 그것에 의하여 α_i와 β_i,j에 대한 값을 최적화하여 적합을 증가시키는 단계,
   g) 미리 세팅된 수렴 규준이 실현될 때까지 단계 e)와 f)를 반복하고, 그것에 의하여 수렴 규준이 실현되었을 때의 최적화된 가중 인자 α_i에 의해 제공된 혼합체의 각 전사체에 대한 추정된 전사체 존재비를 획득하는 단계.
2. 청구항 1에 있어서, 전사체 단편 염기서열결정 데이터는 최소한 5개의 전사체 단편 서열을 포함하는 것을 특징으로 하는 방법.
3. 청구항 1 또는 2에 있어서, 관심되는 유전자 좌위는 하나 또는 그 이상의 유전자 또는 유전자 원소의 하나 또는 그 이상의 동종형을 포함하는 것을 특징으로 하는 방법.
4. 청구항 3에 있어서, 관심되는 유전자 좌위는 하나의 유전자 또는 유전자 원소의 최소한 2개의 스플라이스 변이체를 포함하는 것을 특징으로 하는 방법.
5. 청구항 1 또는 2에 있어서, 전사체의 숫자를 세팅하는 단계는 관심되는 유전자 좌위로부터 미리 주해된 서열 데이터를 획득하고, 그리고 전사체의 숫자를 최소한, 상기 관심되는 유전자 좌위로부터 예상되는, 스플라이스 변이체를 포함하는 상이한 동종형의 숫자로 세팅하는 것을 포함하는 것을 특징으로 하는 방법.
6. 청구항 1 또는 2에 있어서, 확률 부함수 j는 각 유전자 좌표에 대해 양성 값으로 구성되는 것을 특징으로 하는 방법.
7. 청구항 6에 있어서, 확률 부함수 j는 밀도 함수인 것을 특징으로 하는 방법.
8. 청구항 1 또는 2에 있어서, 확률 부함수 j는 가우스 함수, 사각형 함수, 또는 삼각형 함수인 것을 특징으로 하는 방법.
9. 청구항 8에 있어서, 확률 부함수 j는 가우스 함수인 것을 특징으로 하는 방법.
10. 청구항 1 또는 2에 있어서, 유전자 좌표는 관심되지 않는 유전자 영역이 제외되도록 임의선택적으로 변환된, 유전체 내에 뉴클레오티드 위치에 상응하는 것을 특징으로 하는 방법.
11. 청구항 10에 있어서, 관심되지 않는 상기 유전자 영역은 상기 전사체 단편 염기서열결정 데이터에 의한 커버리지를 내포하지 않는 것을 특징으로 하는 방법.
12. 청구항 1 또는 2에 있어서, 스플라이스 접합부를 갖는 유전자 좌표 위치를 상기 커버리지 포락선으로부터 제거하는 것을 포함하는 단계 b2)를 더욱 포함하는 것을 특징으로 하는 방법.
13. 청구항 1 또는 2에 있어서, 상기 단편 유전자 좌표 커버리지는 유전자 좌표에 배정된 각 단편 서열에 대한 최소한 하나의 뉴클레오티드의 수치를 내포하는 것을 특징으로 하는 방법.
14. 청구항 13에 있어서, 최소한 하나의 뉴클레오티드는 단편 시작 부위 또는 전체 단편 서열을 포함하는 것을 특징으로 하는 방법.
15. 청구항 1 또는 2에 있어서, 전사체에 대한 확률 부함수는 각각 상이한 유전자 좌표에서 최대값을 포함하는 것을 특징으로 하는 방법.
16. 청구항 1 또는 2에 있어서, 상기 전사체 단편 염기서열결정 데이터를 제공하기 위해, 최소한 하나의 전사체의 서열 리드를 결정하는 것을 포함하고, 여기서 상기 리드는 상기 전사체의 단편의 서열을 포함하는 것을 특징으로 하는 방법.
17. 청구항 16에 있어서, 상기 전사체는 mRNA인 것을 특징으로 하는 방법.
18. 청구항 1 또는 2에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 5 내지 800개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
19. 청구항 18에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 6 내지 600개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
20. 청구항 18에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 7 내지 400개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
21. 청구항 18에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 8 내지 200개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
22. 청구항 18에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 9 내지 150개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
23. 청구항 18에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 10 내지 100개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
24. 청구항 18에 있어서, 상기 전사체 단편 염기서열결정 데이터의 전사체 단편 서열은 12 내지 70개 뉴클레오티드의 길이를 갖는 것을 특징으로 하는 방법.
25. 청구항 1 또는 2에 있어서, 전사체 i에 대한 각 확률 부함수에 대한 반치전폭 값은 동일한 것을 특징으로 하는 방법.
26. 컴퓨터에서 청구항 1 또는 2의 방법을 수행하기 위한 컴퓨터 프로그램을 포함하는 컴퓨터 판독가능한 기억 장치.

Images