Widespread RNA and DNA Sequence Differences in the Human Transcriptome 2

예전에 이 논문을 소개한적이 있는데 이제야 제대로 봐 볼까 한다. 음.. 사실 제대로 다 볼건 아니고 RDD(RNA-DNA Difference)를 어떻게 찾는지 방법에 대한 것에 초점을 맞추겠다. 이것이 SNP를 찾는 방법과 동일한지 확인해 보고자.


일단 아래 flowchart가 전체적인 분석 과정을 나타낸다. 잘 그렸다. 이것만 봐도 알겠다. 최곤데.

Genotype and SNP calling from next-generation sequencing data

NGS 데이터를 만들면 하는 기본 작업중의 하나가 snp calling & genotyping. 어떻게 보면 이것이야 말로 NGS의 가장 큰 장점이 아닌가 싶다(예전에 랩 모임에서 작은 지웅이가 했던 말이 인상 깊었다. expression estimation을 굳이 RNA-Seq으로 해야 하나하는 문제였는데.. 물론 chip을 이용한 실험보다는 NGS를 이용한 transcriptome 분석에 많이 장점이 있기는 하나 단순 expression 양을 측정하는데 NGS를 사용하는 것은 다른 적용에 비해 얻을 수 있는 장점이 많지 않다는 것). 이 논문은 아래 그림에 SNP calling & genotyping 까지의 과정의 각 단계를 대략적으로 설명한다. 큰 틀을 잡는데 상당히 도움이 되는 논문이다.

review라 대략적으로 알고있는것들이 있는데 그런것은 넘어가고 새로운 것만 기록한다.

Base calling and alignment

base calling and quality scores

illumina 에서 보통 miscall error rate가 1%정도 인데 이것의 원인이 무엇인가 하면 기기 유리판에 증폭되어 같은 dna template가 cluster를 이루는데 이 cluster 안의 서로 다른 DNA copy들의 synthesis process의 desynchronization에 의한것. 그러니까 다르게 말하면 cluster를 이루는 동일한 DNA 가닥들이 한 base씩 합성하면서 fluorescence signal을 나타내게 되는데 이 합성되는 단계의 sync가 맞지않아서 signal이 동일해지지 않는다는 것. 이러한 현상은 synthesis의 cycle이 뒤로 갈수록 심해진다. 

기기 회사에서 제공하는 base calling program 이외에 

• 454 platform : Pyrobayes
• Illumina platform : Ibis, BayesCall
• SOLID platform : Rsolid

의 프로그램이 더 5~30 % 정도 improvement가 있단다.

Alignment and assembly

아시다시피 압축 알고리즘을 이용한 프로그램들(Bowtie, SOAP2,BWA) 가 상당히 빠르나 MAQ. Novoalign, Stampy 등 hash-based algorithm에 비해 sensitivity가 떨어진다, Novoalign과 Stampy 가 현재 가장 정확도가 높다한다니까 한번 봐볼만 하다.

Recalibration of per-base quality scores

SOAPsnp과 GATK의 recalibration에 대해 소개한다. 이것만 봐서는 이해가 힘들다.

Genotype and SNP calling

• SNP calling : polymorphisms의 base position을 찾기 위함
• genotype calling : 각 개인의 genotype을 결정하는 process. SNP calling이 된 위치에서만 가능

예전에는 그냥 단순 allele counting으로 cutoff를 정해서 SNP나 genotype을 calling 했는데 최근에는 probabilistic framework에 uncertainty를 포함하고 더 나아가 allele frequency나 LD pattern 같은 것도 고려사항에 추가 하기도 한다.

Early methods for calling genotypes

NGS를 이용한 초창기 연구에서는 filtering step(>=Q20)을 통해 high confidence base만을 이용해서 고정된 cutoff를 기준으로 SNP를 calling 했다. 예를 들어 Q20 filter를 하고 나서 non-reference allele이 20~80% 사이면 heterozygous genotype으로 그 외에 떨어지면 homozygous로 판단. 이는 deep sequencing일때는 상당히 유용하다. 

Probabilistic methods

moderate 혹은 low depth sequencing에서 초창기 방식을 이용하면 heterozygous genotype이 under-calling 될수도 있으며 Q20 과 같은 filtering을 하게되면 가뜩이나 없는 read가 더 줄게 된다. 또 다른 문제는 단순 cutoff를 정해서 하는 것은 genotype inference에 대한 uncertainty 를 측정할 수 없다는 것이다. 그래서 각 genotype에 quality score를 이용하여 posterior probability를 제공하는 확률적인 방법론이 나오게 된다.

• posterior probability = prior probability * likelihood / constant
• p(G|X) = p(G) * p(X|G) / constant , G : genotype
• prior probability : the probability of genotype으로 reference data에서 이 값을 가져온다.

바로 아래에서 설명하는 방식으로 genotype likelihood를 계산하고 genotype prior(P(G))을 이용하여 Bayes' fomula로 각 genotype G에 대한 posterior probability (P(G|X))를 계산해서 최대의 확률을 갖는 genotype,G1을 그 개체의 genotype으로 할당한다. 이때 P(G1|X) 혹은 P(G1|X)와 P(G2|X) 사이의 비율(G2는 posterior probability가 두번째로 높은 genotype)을 genotype confidence로 사용한다.

Calculating genotype likelihoods

각 read의 quality score를 이용해서 genotype likelihood(P(X|G)) 를 구할 수 있다. 한 개체에서 나온 특정 genomic site의 read i의 genotype G에 대한 확률 p(Xi|G)는 Xi의 quality score의 rescaling으로 계산이 가능하다. 그리고 genotype likelihood p(X|G)는 이 p(Xi|G)들의 product로 계산되어진다. 곱을 이용한다는 건 Xi 사이에 dependence가 없다는 가정이 들어간 것인데  alignment error나 PCR artefact가 있다면  이 가정은 사실이 아니게 된다. 이를 위해 weighting scheme가 제안되었는데.. (이는 잘 안나왔네..)

Martin이 낸 논문에서는 quality score 대신에 read에서 직접 error rate를 측정하는 방식을 제안했는데 이는 quality score의 정확도에 상관없이 error rate를 구할 수 있는 장점이 있지만 이게 얼마나 제대로 작동하는지는 아직 검증이 안되었단다.

Assigning priors using single samples

prior probability는 각 genotype 마다 동일한 확률로 할당하는 방법이 있고 아니면 외부 정보를 이용하여 할당하는 방법이 있다. SOAPsnp의 경우 dbSNP를 이용하는데 특정 genomic site에 dbSNP에서 G/T polymorphism이 알려져 있을때 GG와 TT에 대해서는 0.454를 GT에 대해서는 0.0909 확률을 prior probability로 할당한다. 나머지 genotype에 대해서는 10^-4보다 작은 값을 할당한다.

Assigning priors using multiple samples

Incorporating LD information

A comparison of genotype-call accuracies

SNP calling

여태까지 genotype calling에 대해서만 설명했는데이는 SNP calling과는 약간의 차이가 있다. 예전 NGS 논문들에서는 한사람에 대해서만 분석을 했기 때문에 genotyping이 곧 SNP typing과 동일했다. 여러 사람의 샘플을 갖는 large data set에서는 여느 사람에게 나온 non-refence allele을 SNP라고 할 수 있는데 이는 사람 수가 늘수록 false-positive rate가 늘어날 수 있다는 문제점이 있다. 그담 뭐라 뭐라 하는데.. 당췌 이것만으론 이해 할 수 없다.

Filtering

Incorporating genotype uncertainty in anlayses of NGS data

Recommendations

Conclusions

Evaluation of Gene Structure Prediction Programs

1996년도 논문으로 꽤 오래된 논문인데.. cufflink에서 cuffcompare의 결과로 stats을 뽑아주는데 그 내용이 이 논문과 연관이 있다. 이 논문 완전 보기 싫게 되어 있다. 
이 논문을 보는 목적은 어디까지나 stats을 이해하기 위한것. 


Intro
이 논문에서는 gene structure prediction 프로그램들을 비교 분석한다. coding nucleotide, exon structure, protein product level에서의 프로그램들의 accuracy를 측정하기 위해 gene structure가 알려진 DNA sequence를 테스트 하고 여러가지 performance matrix를 적용한다.


Method 

2.2 Measures of Prediction Accuracy
intro에서 이야기 했듯 accuracy를 3개의 level에서 측정한다.


2.2.1. Coding nucleotide
test sequence에 대한 predicted coding value와 true coding value를 비교한다. coding value라는 것은 coding 혹은 noncoding 두 값중 하나를 값게 되는 binary variable이고 predicted와 reality의 coding value를 표현하기 위해 위 그림의 table과 같이 2X2 contingency table을 이용한다. table의 각 cell안의 값은 bp이다. 일단 이 table 안에 reality와 prediction 사이의 관계에 대한 많은 정보가 담겨 있는거고 이를 single scalar metric로 표현하는 데는 많은 방법이 있다. sensitivity(Sn)와 specificity(Sp)가 대표적인 표현식.

쉽게 말하면 Sn은 실제 coding nucleotide중 얼마나 제대로 prediction 됐나의 비율이고 Sp는 noncoding nucleotide중 얼마나 제대로 noncoding이라고 prediction 되었나의 비율.

상대적으로 noncoding 영역이 많기 때문에 TN이 FP에 비해 상당히 크기 때문에 위의 식으로 구한 Sp는 크게 의미 없다. 그래서 대체식으로 아래 식을 이용한다.

위식의 의미는 coding region이라고 prediction 한것 중에 제대로 prediction 한 비율이 얼마만큼 되나. 뭐 이런 의미라 할수 있겠다. 이를 조건부 확률식으로 표현하면 아래와 같다.

• Sn = P(F(x)=c|x=c), Sp = P(x=c|F(x)=c)  #x는 실제 상태, F(x)는 예측된 상태를 의미

그런데 Sn과 Sp 하나로만 accuracy를 측정하기에는 완벽하지 않다(Sn은 높고 Sp는 낮은 경우 혹은 반대. 예를 들어 어떤 프로그램이 대부분의 nucleotide를 coding region이라고 prediction 했으면 이는 Sn은 상당히 높게 되고 Sp는 상당히 낮게 된다). 그래서 이 둘을 summarizing 하는 single scalar value가 도입되는데 이것이 Correlation Coefficient(CC). 

그런데 CC는 contingency table의 cell 중 하나라도 0 값을 갖으면 구할 수 없다는 단점이 있다. 이를 보안하기 위한 한 방법이 SMC(simple matching coefficient) 값을 이용하는 것. 

다르게 표현하면 SMC = P(x=c)P(F(x)=c|x=c) + P(x=n)P(F(x)=n|x=n)이 된다. 일반적으로 P(x=n)이 P(x=c)에 비해 값이 상당히 크기 때문에 SMC는 noncoding에 대한 sensitivity에 영향을 많이 받게 되는 단점이 있다(1000bp중 100bp가 coding region인 sequence가 있다고 하자. 첫번째 prediction은 100bp를 걸치는 250bp의 coding region을 prediction 했고, 두번째 prediction은 noncoding 지역에다가 50bp인 coding region을 prediction 했을때 두 prediction의 SMC 값은 0.85로 동일하다).

결론적으로 CC처럼 bias가 없으면서 SMC 처럼 어떤 상황에서든 구할 수 있는 single scalar value가 필요한데 그것이 바로 ACP(average coefficient probability)이다. 

• ACP = 1/4[ TP/(TP+FN) + TP/(TP+FP) + TN/(TN+FP) + TN/(TN+FN) ]

즉, 용어 그대로 ACP는 각 conditional probability의 평균값을 의미한다. 이 ACP를 확률 값이므로 0과 1 사이의 값을 갖으므로 CC처럼 사용하기 위해서는 AC(approximate coefficient)를 이용하는데 그 식은 아래와 같다.

• AC = (ACP - 0.5) X 2       #AC는 CC 처럼 -1에서 1 사이 값을 갖는다.

2.2.2. Exonic structure
test sequence 내에 predicted exon과 실제 exon을 비교한다. 다만 exon이 correctly match 됬다는 기준이 좀 애매한데 여기서는 완전히 exactly match 된 것만을 match 됬다고 여긴다(원한다면 threshold를 정해서 그 이상으로 overlap된 exon은 match 됐다고 해도 된다). accuracy를 측정하기 위한 방법으로 specificity와 sensitivity를 사용하는데 coding nucleotide level에서의 사용했던 것과 동일하다(단 coding nucleotide level에서는 bp가 단위였는데 exonic structure level에서는 exon 갯수가 단위). 위 그림 참조.


coding nucleotide level은 프로그램의 content component를 exonic structure은 signal component에 의한 것이라 볼 수 있겠다(말 그대로 coding nucleotide는 bp 단위로 얼마나 잘 content가 커버 됐나를 나타내는 것이고 exonic structure는 start stop codon과 같은 boundary가 정확히 맞는 것을 보기 때문에 sequence의 atomic signal을 측정한다는 것).
2.2.3. Protein product

Identification of novel transcripts in annotated genomes using RNA-Seq

CuffLinks 논문중 하나. CuffLinks의 한 옵션인  -g 옵션의 설명을 위한 논문. RABT(reference annotation based transcript assembly)를 설명하는 논문이다.




INTRODUCTION
RNA-Seq의 탄생으로 좀 더 정교한 transcript abundance estimation이 가능해졌다. 그런데 그것 뿐만 아니라 RNA-Seq을 EST로 봄으로써 genome annotation도 가능함을 이야기 한다. 이런 RNA-Seq으로 annotation하는데 있어서 발현양이 적은 transcript의 경우에는 transcript의 부분만이 커버가 된다는 문제점이 있다. 아래 그림에서와 같이 95%이하로 cover된 transcript가 64.4%를 차지함을 볼 수있다. 이렇듯 naive assembly 방식으로는 제대로된 transcript를 구축하는데 문제가 있을 수 있다. 

기존의 transcript assembly를 아래 3가지로 나눌수 있는데 그 어떠한 방식도 기존의 annotation을 이용하지 않는다.

1. De novo assembly
2. Genome reference based transcript assembly : genome에다가 mapping 하고 assembly 한다. Scripture 프로그램에서 했던 것. cover된 영역 사이의 gap을 채우는 것. 하지만 기존에 존재하던 annotation 파일을 이용하지는 않았다는 한계점이 있단다.
3. RNA-Seq assisted protein coding gene annotation : ab initio gene finding program의 증거로다가 RNA-Seq데이터를 이용하는 것.

그래서 이 논문에서 기존의 annotation을 assembly에 이용하는 방식인 RABT assembly 방식을 소개한다.

METHODS

위 method를 이해하기 위해서는 먼저 cufflinks의 transcript assembly를 봐야한다.
Cufflinks Transcript Assembly
cufflinks의 transcript assembly는 다음과 같은 목적을 갖는다.

1. 모든 fragment(paired-read)는 최소 한개의 assembled transcript와 일치한다
2. 모든 transcript는 read들의 tiling에 의한 것이다
3. transcript의 갯수는 위의 목적 1번을 만족하게 하는 최소한의 transcript 수이다.
4. the resulting RNA-Seq models are identifiable

결국 fragment를 설명하기 위한 최소한의 갯수의 transcript를 찾아내는 것이 cufflinks의 목적.

reference annotation을 assembly algorithm에 이용하기 위해 3가지 접근을 채택했다(위의 그림이 overview).

1. reference transcript로 부터 가상으로 faux-read를 tiled 되게 만들어서 low coverage transcript의 missing 된 부분을 찾을 수 있도록 한다. reference transcript의 15bp마다 405bp 길이의 faux-read를 생성.
2. read와 1번 step의 faux-read 모두를 이용하여 Cufflinks로 parsimonious assembly를 한다. 
3. 2번 step에서 생긴 transfrag들을 reference transcript와 비교해서 다음 다섯가지 모든 조건과 일치하면 제거한다.
• 5' endpoint가 reference transcript에 포함된 경우
• 3' endpoint가 reference transcript의 600bp 이상으로 뻗지 못한 경우(그리고 이 영역에 intron이 없는 경우)
• reference transcript에 없는 intron을 가지지 않은 경우
• a
• a

RESULTS


DISCUSSION

deFuse: An Algorithm for Gene Fusion Discovery in Tumor RNA-Seq Data

RNA-Seq 분석 과정의 일환으로 fusion gene을 찾고자 선택된 프로그램 논문. 프로그램 다 돌려보고 나서 논문 읽어 볼려고 했는데 프로그램이 중간에 error가 나서 stop 한 관계로.. error 원인을 찾다가 run_adaboost.R의 input file에서 NA가 존재하는 것이 원인인 것을 알고 이 record를 지워도 되는가 판단하고자 급하게 논문을 읽어본다.


abstract를 보자면..
defuse는 unique read align만 쓰는게 아니라 모든 align을 쓰고 exon의 끝부분의 fusion만 보는것이 아니라 모든 부분에서의 fusion을 본다고 한다. 그래서 더욱 sensitive 하다고. specificity를 높이기 위해서 RT-PCR로 confirm 한 novel feature로 adaboost classifier를 train 한다는데 이건 영뭔소린지 본문을 봐야 알겠다. 뭐 그리그리 해서 ovarian cancer에서 gene fusion을 찾았단다.


gene fusion은 double stranded DNA breakage이 DNA repair error에 의한 것이라는데 이건 이 논문을 봐보자.


뭐 여튼.. defuse는 다른 프로그램과는 달리 ambiguously aligning read 도 쓰고, split read와 discordant read를 이용하는데 그 순서를 기존 프로그램과 달리 discordant read 분석을 먼저 하고 그 담에 dynamic programming-based split read analysis가 들어간단다. 이 순서가 좀더 sensitive 하다고.


The deFuse algorithm


일단 용어 정리부터 하자면

• fragment : a size selected cDNA sequence during RNA-Seq library construction
• read : fragment에서 sequencing 된 부분
• insert sequence : fragment에서 paired-read를 제외한 sequencing 안된 부분
• fusion boundary : nucleotide 단위의 genomic 위치로 gene fusion이 일어나는 양쪽의 break point, 그러니까 두 유전자가 gene fusion이 일어날때 합쳐지는 그 bp위치
• spanning reads : paired-ends 사이의 insert sequence에 fusion boundary가 있는 reads
• split read : read 안에 fusion boundary가 있는 read
• discordant alignment : spanning reads와 split read의 alignment을 의미. spanning reads의 경우 paired-ends가 서로 다른 유전자에 align이 될거고, split read의 경우 한쪽 부분이 유전자의 끝부분에 align될거고 나머지 끝부분은 align이 안됨 

defuse의 개요가 아래 그림과 같다. 총 4단계로 구성되어 있는데.. 

1. read를 reference에 align 한다. 이 때 reference로는 spliced & unspliced gene을 모두 사용. 어떤경우에는 unspliced region이 intron인 gene fusion에 의해 발현되는 경우도 있기에. 이때 두가지 기준을 정해서 동일한 fusion event를 나타내는 discordant alignment를 clustering 한다.
2. 가장 그럴듯한 fusion event를 선택한다. 아래 그림 a
3. 각 fusion event의 fusion boundary를 찾기 위한 dynamic programming based solution에 이용될 split read를 찾는다. 아래 그림 b
4. spanning & split reads 증거들을 증명을 위한 test. split read로 찾은 fusion boundary를 바탕으로 paired-end의 fragment의 putative length를 계산하여 fragment length distribution을 기준으로 차이가 있는지를 확인한다. 그리고 나서 quantitative feature를 구한뒤 adaboost classifier로 진짜와 가짜를 구분한다. 그림 c

Conditions for considering discordant alignments to have originated from reads spanning the same fusion boundary.
defuse 알고리즘의 가장 첫번째는 fusion event가 일어나는 region을 걸치고 있는 spanning read를 찾는것. 그렇기 때문에 spanning read 의 선택 조건을 정하는 단계이다.
concordant paired-end로 fragment length distribution인 P(L)을 알 수 있는데, 여기서는 [lmin,lmax] 길이 안의 fragment만 고려한다. 

• lmin : a/2 percentile of P(L)
• lmax : (1-a/2) percentile of P(L)
• a : proportion of paired end reads that are not guaranteed by the algorithm to be assigned to the correct fusion event

서로 다른 discordant alignment가 동일한 fusion boundary를 포함하는 spanning read인지를 판단하는 2가지 조건. 

• c1 : overlapping boundary region condition : 동일한 gene fusion event에서 나온 두 paired-end의 fusion boundary region은 반드시 overlap 되어야 한다. 아래 그림 c
• c2 : similar fragment length condition : 두 fragment의 길이차이(=dx+dy)가 lmax-lmin보다 작아야 한다.

아래 그림을 좀더 설명하자만 a와 d는 fusion boundary를 알고 있을때의 그림이고 b,c,e는 정확한 fusion boundary를 모를때의 위 조건에 대한 설명이다. 그러니까 실질적으로 a와 d와 같은 model에서 paired-end가 나왔을텐데 아직까지 boundary를 모르니까 b,c,e와 같은 그림이 그려질거고 거기서 부터 조건 c1과 c2가 나왔다고 생각면 된다.

Assigning a unique discordant alignment to each spanning read.
spanning read에 의해 가능한 모든 discordant alignment 중(=valid cluster)에서 split read를 찾기 위한 후보가 되는 valid cluster를 찾는 단계이다.
ambiguous alignment : spliced & unspliced gene sequence를 reference로 mapping 할때 gene들의 homology에 의해서 ambiguous alignment가 나타날 수도 있고 똑같은 유전자의 alternative splicing에 의한 동일한 exon들이 multiple splice variant(=isoform)에 나타남에 의할 수도 있다.
여튼 이러한 ambiguous alignment들에서 부터 제대로 된 alignment를 찾아야 한다.

• valid cluster : discordant alignments set으로 이 set 안의 모든 두 paired-ends 는 조건 c1과 c2를 만족한다. discordant 하게 alignment된 paired-end 의 집합인데 이 안의 모든 원소들은 서로 overlap(조건 c1)되고 fragment 길이 차이가 최소한보다는 작다(조건 c2).

ambiguous alignment에 의해 동일한 paired-end가 여러개의 valid cluster 속할 수 있다. 이렇기 때문에 하나의 ambiguous aligned paired-ends를 하나의 valid cluster에 할당함으로 해서 fusion event를 최소화 한다(maximum parsimony solution 이라고 하는데 supplementary를 참조할 필요가 있다).  모든 read와 모든 valid cluster를 unselected라고 해놓고 알고리즘의 각 step에서 read를 가장 많이 갖은 valid cluster를 selected 라고 해놓고 그 안의 read들도 assigned라고 해놓는다. 이 같은 step을 반복해서 valid cluster를 뽑는데 이를 maximal valid cluster라 한다.

Split read boundary sequence prediction.
찾아진 valid cluster 가 영유하는 영역(=approximate fusion boundary)에 align이 될거라 예상되는 split read를 찾고 그 split read와 approximate fusion boundary를 dynamic programming으로 align 해서 fusion boundary를 찾는 단계이다. 
위 단계까지 해서 fusion event가 일어났을 것이라 예상되는 region을 찾으면 bp 단위의 fusion boundary를 찾기 위해 targeted split read analysis를 한다. 

• approximate fusion boundary : 같은 valid cluster 안에 속하는 discordant alignment들의 fusion boundary region 의 intersection region. 아래 그림 a
• candidate split read : 한쪽 부분이 approximate fusion boundary 에 anchored된 paired-end read
• mate alignment region : candidate split read는 paired-end 중 한쪽만 anchored 된 read의 반대쪽 discordant read 가 고려 대상인데 한쪽이 align되었을때 align이 안된 read가 align 될거라고 예상되는 영역, 아래 그림 b

candidate split reads는 mate alignment region과 approximate fusion boundary가 겹치는 곳에 위치하게 되는 read들이 된다.
the split read analysis는 candidate split read를 approximate fusion boundary에 align 함으로써 진행된다.  candidate split read의 fusion boundary에 의해 split 되었을 꺼라 예상되는 read(그러니까 anchored read 말고 반대쪽 read)를 transcript X의 approximate fusion boundary(=Sx)와 align 하고 transcript Y의 approximate fusion boundary(=Sy)와는 reverse로 align 한다(transcript X, Y는 gene fusion이 일어났을 거라 예상되는 trascript). 이때 align을 dynamic programming을 이용하기 때문에 X,Y 각각에 대해 matrix가 생성될 것인데 이를 Dx,Dy 라 한다. 이 allign에서 split을 (ix,iy,j)로 표현 할수 있는데 ix와 iy는 Sx와Sy에서 fusion boundary의 bp position을 뜻하고  j는 read 상에서의 bp 위치를 뜻한다. 이 (ix,iy,j)는 아래 공식으로 찾는다.

이때 Dx와 Dy는 threshold m_anchor를 넘어야 하는데 m_anchor = m*n_anchor 이다(m은 match score, n_anchor은 Sx 혹은 Sy에 align되어야 하는 최소한의 bp 갯수, 곧 Sx나 Sy 한쪽에만 너무 align 되는 것을 방지하기 위함이다).
여기서 나타날 수 있는 문제점이 여러(ix,iy,j)가 같은 최대값을 갖는 경우이다.  이는 여러 split read들의 fusion boundary (ix,iy)를 clustering 해서(=동일한 (ix,iy)를 갖는 split read를 clustering) 이 cluster 중 read의 anchoring score의 합의 최대인 것을 뽑는것으로 해결한다.




Corroborating spanning read and split read evidence.
최종적으로 찾아진 fusion boundary를 confirm 하는 단계로 split read의 align으로 찾아진 fusion boundary를 바탕으로 spanning read의 fragment length distribution과 concordant read의 fragment length distribution과의 비교를 통해 유의한지 판단한다.
spanning read evidenced와 split read evidence 사이의 일치성을 확인한다. 일단 split read에 의해 예측된 fusion boundary를 이용하여 spanning read의 fragment length를 추론한다. 이렇게 나온 spanning read의 fragment length의 분포(={li})와 concordant paired-end에서 나온 fragment length의 길이 분포(P(L))를 z-test 한다({li}가 P(L)에서 나왔다는 가정하에). 


----------------------------reference--------------------------------
solid tumor : cysts와 liquid에 상관없는 비정상적인 조직덩어리. 정확히는 모르겠지만 말그대로 액체에 의해 커진게 아니라 조직 자체가 비정상적으로 커져버린 덩어리를 뜻하는것 같다. 악성일수도 있고 양성일 수도 있단다. 3가지 type(sarcoma, carcinoma, lymphomas)의 solid tumor가 있다. 
sarcomas는 뼈와 근육같은 connective tissue에서부터 생성된 tumor type.
carcinomas는 glandular and epithelial cell에서부터 형성된것(샘선이나 상피세포). 이런 cell들은 공기의 유출입이나 위장과 관련된 cell. lymphomas는 림프절에 관련된 tumor 일것이고.

Human DNA methylomes at base resolution show widespread epigenomic difference

전에 포스팅한 논문에서 quality assessment에서 이번 논문의 방식을 따랐는데.. incomplete bisulfite conversion에 의한 putative non-CpG methylation을 구별하고자 binomial test랑 false discovery rate를 이용했는데.. 이거 또 봐야하지 않겠나해서 본다.

abstract를 보자면 처음으로 single-base resolution DNA methylation을 본 논문이란다. 동시에 mRNA, small RNA, histone modification, DNA-protein interaction도 봤단다. 돈 많이 썻다. 두개의 genome(human embryonic stem cell, IMR90 fetal lung fibroblast)의 methylation을 비교했는데 stem cell의 경우 methylation의 1/4 가량이 non-CpG 에서 나타났단다. 엄청나다. 그래서 이것이 아마도 stem cell에서는 완전 다른 methylation mechanism이 있는게 아닐까 한단다. 재밌는건 non-CpG methylation이 gene body 부분에 enrichment되어 있고 protein binding site과 enhancer에는 depletion 되어 있다는 것. 더욱 재미있는건 이와 같은 non-CpG methylation이 differentiation 되면서 잃어가는 반면 induced pluripotent stem cells에서 다시 나타난다는 것.

내용이 재밌어 보이기는 하지만 시간 관계상 일단은 내가 보고자 하는 것만 본다. supplementary information을 보자. 


Data Analysis에 보면 처음에 read의 3단계의 preprocessing 과정과 또 3단계의 post-processing 과정이 나온다. 이건 패스 하고.. 


Identification of methylated cytosines를 보면 내가 찾고자 했던 binomial test에 대해 나온다. 일단은 binomial test¹ 에 대한 기본적인 내용은 아래 reference를 참조하고.. 이 논문에서는 unmethylated lambda genome을 spiked(이거 직역이 뭔지 모르겠네) 했는데, 이 lambda genome에 align 되는 read를 가지고 error rate(그러니까 unmethylated genome을 BS-Seq 했으니까 read들의 모든 cytosine이 T로 변환되어 있어야 정상인데 그렇지 않은 cytosine 을 전체 cytosine의 갯수로 나누어서 error rate을 구함)를 구하고 이를 binomial distribution B(n,p)의 확률 p로 사용한다. 그리고 n은 각 cytosine position에서의 read depth. 음 여기서 확실하게 이해가 안되는 부분이 있는데.. 여기 표현을 빌리자면 0.01 FDR corrected P-value를 이용해서 methylation 의 threshold를 구하고 methylation 됐다 안됐다를 따지는데..


다시 말하면 B(n,p) n은 각 cytosine position에서의 read depth, p는 위의 error rate. 한 cytosine position에서 cytosine 이라고 sequencing 된 read(k; cytosine이라고 나왔다는 건 conversion이 안됐다는 거고 곧 methylated cytosine을 의미) 의 p-value가 M보다 작아야 하는데 M은 M*(number of unmethylated cytosines) < 0.01*(number of methylated cytosines) 를 만족하는 수. 내가 이해하기로 곧 p-value가 0.01*(# of methylated cytosines)/(# of unmethylated cytosines) 보다 작아야 이를 methylated 됐다고 판단 한다는 것. 이해가 안가는게 바로 이점. 여기서 0.01을 곱했다고 해서 0.01 FDR corrected P-value라고 표현한 것 같은데 이게 이해가 안된다.. 아.. 


아니면 내가 전혀 딴소리 하고 있는건가..




--------------------------------------reference------------------------------------------
1.binomial test : 일단은 두개의 category가 일어날 확률이 비슷하다는 것이 귀무가설일때 사용하는 test 란다. 예를 보자. 주사위를 이용한 보드 게임을 한다고 할때, 235번 주사위를 던져서 51번 6이 나왔을때 과연 이 주사위가 제대로 된거냐(fair dice)? binomial test를 이용해서 이 문제를 해결하기 위해 이항분포 B(235,1/6)을 이용한다. 이 이항분포에서 6이 52 번 이상이 나올 확률을 더한다. 그니까 이항 분포니까 평균이 np 분산이 npq, 곧 235*1/6=39 가 평균, 235*1/6*5/6 이 분산인 정규 분포를 따르게 되고 one-tailed test로 하자면 51을 이상에서의 확률 밀도 함수의 면적이 0.026544 가 나오니까 이는 상당히 significant 하다.. 뭐 이런..