RNASEQR-a streamlined and accurate RNA-seq sequence analysis

ISB에서 나온 논문

Evaluation of Algorithm Performance in ChIP-Seq Peak Detection

제목 그대로 ChIP-Seq 프로그램 performance를 비교한 것인데.. 음 이런 논문이야 말로 짐 회사에서 내기 좋은 주제가 않을까 싶은데.. 아숩다.

ChIP-Seq 데이터 분석에서 필요한 peak finding을 위한 프로그램이 31개나 있단다. introduction에 보면 ChIP-Seq 분석 프로그램의 대략적인 알고리즘 개요가 나온다. NGS 특성상 5' 쪽의 tag만 읽기 때문에 생기는 strand-dependent bimodality를 보정하기 위한 방법(paired-end 는 몇개의 프로그램에서만 지원된단다), read가 많이 mapping된 genomic region을 찾는 방법, peak region을 찾기 위해 threshold를 정하는 방법(background signal model 이용:1.manual threshold,2.Poisson or negative binomial model 이용,3.control data 이용), peak의 significance를 정하는 방법에 대한 여러 알고리즘의 간략한 소개가 있다.

일단 여기서는 11개의 peak calling algorithm을 3개의 transcription factor ChIP-Seq 데이터를 가지고 비교한다. 이것이 이 논문의 목적.

Results
overview
3가지 dataset은 NRSF(human neuron-restrictive silencer factor), GABP(growth-associated binding protein), FoxA1(hepatocyte nuclear factor 3a)의 ChIP-Seq데이터. 
테스트한 프로그램들 11개의 리스트는 아래와 같고 각 프로그램의 option들은 default를 이용(control data를 이용할 수 있는 프로그램만 선택).

Sensitivity
3개의 dataset에 대해 11개의 program에서 찾는 peak의 수는 차이가 있다. 

How to Interpret a Genome-wide Association Study

GWAS는 candidate gene analysis(supervised analysis)와 family linkage study, 그리고 HapMAP Project의 성과를 바탕으로 이루어진 것이다.
GWAS는 common disease¹, common variant 의 가정하에 있는데 이는 common disease의 유전적인 영향 혹은 원인이 제한된 수의 allelic variant (SNP 혹은 indel 등)에 의한 가정인데 이 제한된 수의 allelic variant 라는 것이 전체 인구에서 1% 혹은 5%이상의 사람들이 가지고 있는 allelic variant를 의미한다.


Overview of GWA Studies
GWAS는 NIH에서 관측된 특성(질병등)의 유전적인 연관성을 찾기 위해 사람 전체 genome에 걸쳐 common genetic variation을 연구 하는 것이라고 정의 되어 있다. genome wide 라는 것의 정확한 기준은 없지만 이 논문에서는 최소한 1,000,000 SNP를 assay 한 연구에 대해서만 언급한다.
GWAS는 크게 4가지 부분으로 이루어진다(PLINK manual에서는 크게 6가지 단계로 나눈다).
1.특정 trait(질병등)를 갖는 집단군(사람들)과 그렇지 않은 집단군을 선택한다.
2.위 단계에서 뽑은 모든 사람을 genotyping 하고 genotyping quality에 대해서 review를 한다
3.2번 단계에서 quality threshold를 넘는 SNP들 중 어느 SNP가 trait와 연관이 있는지 통계적 테스트를 진행한다.
4.3번단계까지 해서 뽑은 genetic variant의 증명단계로 완전 새로운 집단군을 뽑아서 동일한 association test를 진행하거나 아니면 실험적으로 기능적인 효과가 있는지 테스트 한다.


Study Designs Used in GWA
여태껏 가장 많이 사용되었던 GWAS의 design은 case-control design으로 환자군과 정상군의 allel frequency의 비교였다. 이는 가장 간단한 design인데 아래 표와 같이 많은 가정 하에 design 된 연구고 만약 이 가정이 충족되지 않게 된다면 연구의 결과는 상당한 bias가 있게 된다. 그러나 역학적인 디자인의 원리에 충실하게 연구가 디자인 되었다면 case-control은 rare disease에 대한 효과적인 연구가 될 수 있다. 하지만 이게 쉽지 않다는 것.
trio study은 환자와 함께 환자의 부모를 포함한 design이다. 환자인 자식과 그리고 부모의 genotyping을 해서 transmission frequency를 측정한다. 이게 무슨 뜻이냐면 만약 특정 SNP가 질병과 관련이 없다면 부모에서 자식으로의 이전률이 50%일거지만 질병과 관련 있는 SNP는 그 이상일거라는 가정하에 test를 하는 것이다.
또 다른 design은 cohort  study이다.

Selection of Study Participants

Genotyping and Quality Control in GWA Studies

Analysis and Presentation of GWA Results

Replication and Functional Studies

Limitations of GWA Studies

Clinical Application of GWA Findings



-----------------------------------reference---------------------------------
1.common disease : mendelian disease이외의 것, mendelian disease 혹은 mendelian disorder 라는 것은 멘델의 유전법칙을 따르는 질환으로 DNA 상의 하나의 mutation으로 인한 질환.

De novo Transcriptome Sequencing in Anopheles funestus using illumina RNA-Seq technology

Bioconductor: open software development for computational biology and bioinformatics

- primary motivations -

• transparency : entire process가 확실하게 노출되어야 한다.
• pursuit of reproducibility : algorithmic work에도 standard가 필요.
• efficiency of development : 기존의 code의 extension과 novice의 발전을 위해 필요.

- seven topics important to establishment of a scientific open source software project -

1. Language selection : 왜 R을 선택 했냐?
-prototyping capabilities : 빠르게 prototype을 만들수 있다. 물론 나중에 더 빠르게 run 할 수 있도록 re-implement도 가능하다.
-packaging protocol : package 형태로 제작, 테스팅, 배포가 가능하다.
-object-oriented programming support : To secure reliable package interoperability
-WWW connectivity : http와 같은 web resource롤 통해 데이터와 package에 접근 가능하며 XML처리하는 package도 있어서 다양한 데이터를 다룰 수(?perceive) 있다.
-statistical simulation and modeling support : R에서 이미 있는 numerical algorithm의 사용이 용이하다.
-visualization support : graphical tool로서의 기능이 좋다.
-support for concurrent computation : parallel computation을 위한 tool들이 있다.
-community : active user and developer communities

2. Infrastructure base 
Bioconductor project 에서 첫 2년은 software infrastructure의 투자에 집중한다. 이 infrastructure는 reusable data structure & software 형태로 만든다. 
이 software infrastructure concept의 두 예로 Biobase package의 "expreSeq" class와 Bioconductor metadata package 중의 하나인 hgu95av2를 들 수 있다.
expreSeq 은 three-tier architecture 를 용이하게 한다. 뭔말인고 하니 low-level processing software designer는 expreSeq instance만 생성하는데 focusing 하면 되고, 분석가는 low-level processing에 신경쓰지 않고 expresSeq 자료구조만 focusing 해서 분석만 하면 된다. 
hgu95av2는.. 잘 모르겠네..

3. Design strategies and commitents
-designing by contract
-object-oriented programming
-modularization
-multiscale and executable documentation
-automated software distribution

4. Distributed development and recruitment of developers
distrubted development의 강조. CVS를 통한 같은 component에 대해 여러 developer가 개발에 참여하게 됨으로 다양한 viewpoint와 experience가 project에 속하게 된다. 한 사람의 개발자에 의한 code의 변화가 다른 코드를 망가지게 하지 않는 것을 원칙으로 한다. 이는 사실 packaging화로 가능하다. 그리고 이 R package의 규격화된 testing system을 제공함으로 인해 안정적인 개발이 가능케 한다.

5. Reuse of exogenous resources  
다른 project의 software를 adapting 하는데 있어서의 3가지 쟁점
-가능하면 re-implementation 하지말고 있는거 갖다 쓰자.
-CBB(computational biology & bioinformatics)는 다양한 분야를 아우르기 때문에 다른 많은 프로젝트와의 공동의 노력이 필요하다. 그렇기 때문에 다른 언어나 시스템에서 쓰여진 데이터나 알고리즘을 사용하기 위한 구조화된 패러다임이 필요하다.
-standardization and reuse of existing tools

6. Publication and licensing of code

7. Special concerns
CBB에서 생기는 4가지 challenge
-reproducible research
-dynamics of biological annotation
-training
-responding to user needs




- Using Bioconductor (example) -


ALL(Acute lymphocyte leukemia)

Enrichment of differentially methlyated regions with MethylMiner fractionation and deep sequencing with the SOLID System

invitrogen 에서 MethylMiner 라는 kit의 application note로 publish 한 것(여기). 이걸 본 이유는 elution의 농도에 따른 bias를 어떻게 해야 하는가에 대한 의문에서인데.. 
생각보다 MBD-Seq의 장점을 많이 알 수 있게 되서 posting 한다(물론 invitrogen에서 자기네 kit 홍보성으로 만든 것이기에 공정성은 좀 떨어지지만 그래도..).




MBD-Seq의 장점 (MeDIP-Seq 대비)


antibody로 DNA methylation을 detecting 할려면 DNA 를 denature 시켜야 한단다. antibody가 mCpG에 fully access해야만 하기에 denaturation 시키는 과정이 들어가게 되고 이는 실험적인 소요가 요구된다. 음.. 전혀 생각 못했던거다. 워낙에 논문들 대세가 MeDIP-Seq 이라서 그냥 생각없이 antibody를 이용한 precipitation이 좋을 줄 알았는데.. 
글고 여기서 보여지금 그림 하나.

위 그림이 MeDIP 처리 한거랑 MethylMiner 처리한거랑 dye를 붙여서 chip으로 찍었을때 어느 것이 enrichment되나를 본건데 MethylMiner를 이용한것이 훨씬 많은 array hit에서 enrichment 가 되어 있음을 알수 있다. 단순하게 생각한다면 MethylMiner가 훨씬 sensitive 한것처럼 보인다.

elution salt 농도는 어쩌라는 거냐?

음.. 여기 있는거 가지고 굳이 결론을 내려보자면(엄연히 내 생각이고.. 실험에 대한 이해 부족으로 잘못된 결론일 수 있다는 가능성은 농후하다)..

"Greater than 90% of the captured DNA is sequentially eluted with 500 mM and 1M NaCl"

이란 언급에서 500mM이랑 1M만으로 sequentially elution 했음에도 거의 대부분의 captured DNA가 elution 됨을 알수 있다. 그렇기 때문에 gradient로다가 elution을 하게 되면 거의 모든 captured DNA 가 elution 될 것이므로 gradient에 따른 특정 CpG density 영역의 read들만 elution 되는 일은 없을 것으로 생각되어 진다.


그리고 반드시 elution 농도를 섞어에 하는게 아닌가 싶다. 이 application note에서는 500mM과 1M 만을 비교 했는데.. 언뜻 보면 1M 이 좋아보이는가 싶지만 마지막 그림에서 보이듯이 500mM에서만 유독 enrichment 되는 부위가 있다. 곧 모든 elution 농도로다가 뽑아낼 수 있는 DNA는 전부 뽑는게 가장 좋은 방법이 아닐까 한다만... 모르겠다.

근데 또 이 논문을 보긴 봐야 할거 같은데.. 언뜻 method 만봐서는 잘 이해가 안가긴 하지만..

Comparative methylome analysis of benign and malignant peripheral nerve sheath tumors

이번 보건원 미팅 저널 발표 자료.


음.. 뭐랄까.. DNA methylation landscape의 진화라고 할까. 분석의 대부분이 genomic structure의 methylation 변화를 관찰한거였는데.. satellite repeat을 또 subtype으로 나눠서 관찰했다. 


첨에 읽을 땐.. 이 논문 별룬데 하고 생각했는데. 그래도 이 정도까지 이야기꺼리를 만들었다는게 참 훌륭하달까.


아.. 맘에 들었던 점은 cnv를 찍고 methylation을 관찰했다는 것. genomic copy number에 의한 methylation enrichment가 아니라는 것을 확인했다는게 맘에 들었고.. 남의 expression data 가져와서 clustering 한것도 괜찮았다. 아.. 글고 저자가 답장을 굉장히 빠르게 해줘서 음.. 맘에 든다.

A Bayesian deconvolution strategy for immunoprecipitation-based DNA methylome analysis

MBD-Seq 을 분석하기 위해서 결정한 MEDIPS R package에서 많이 언급하는 논문이다. 일명 Batman이라는 프로그램을 소개한 논문. 그리고 보건원에서 발표할 다음 논문에서도 이용한는 프로그램이다. immunoprecipitation-based DNA methylation profiling에서 absolute methylation level을 표현하기 위한 algorithm이 소개된다.




서론서는 immunoprecipitated-based DNA methylation profiling이 가격대비 성능이 젤 좋다고..

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 하다.. 뭐 이런.. 

The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants

FASTQ 형식이 여러가지다. bowtie를 돌릴 때도 --phred33-quals, --phred64-quals, --solexa-quals, --solexa1.3-quals 네가지를 인자로 받는다(default 가 --phred33-quals이고, solexa1.3과 phred64는 동일하다). 이 논문을 전에 본적이 있는데 대충 봐서.. 기억이 안나는 관계로 확실하게 하기 위해 다시 보기로 한다.


<PHRED scores and the qual format>
PHRED 라는 프로그램이 있는데 이는 DNA sequencing trace file을 input으로 받아서 base calling을 하고 거기에 대한 quality를 할당한다. 

여기서 Pe는 estimated probability of error. 이를 PHRED는 QUAL format 파일로 저장하는데 이는 fasta와 비슷한 형식으로 아래와 같다.

이 PHRED score는 사실상 base quality를 표현하는데 표준이 되었다. 이는 SAM, ACE, FASTQ에서 사용된다.


<Sanger FASTQ format>
FASTQ format은 sanger 센터에서 만들어 졌단다. 뭐 파일 형식이 어떻고를 떠나서 여기서는 어떻게 quality score가 string으로 encode 되는지에 초점을 맞춘다. 초창기 FASTQ format을 사용했던 Sanger capillary sequencing은 바로 위의 PHRED quality score를 따른다. 그리고 나서 이 score를 single character로 변환해서 파일에 담는데 이때 charater로 ASCII 33-126 번까지의 문자를 이용한다. 곧 표현 가능한 quality score가 0-93. OBF project (Open Bioinformatics Foundation)은 이 format의 이름을  fastq-sanger 라 한다.


<Solexa FASTQ format>
2004년에 Solexa Inc. 에서 소개한 FASTQ format. solexa quality score는 아래와 같다.

phred quality와 solexa quality 간의 변환 공식은 다음과 같다.

solexa score는 -5부터 값을 갖을수 있기에 64를 offset으로 정해서 ASCII 59-126 까지를 쓰는 것으로 한다. OBF projects에서는 이 format을 'fastq-solexa' 라고 한다.

<Illumina 1.3+ FASTQ format>
solexa가 Illumina로 팔리고 나서 GA(Genome Analyzer Piepine) 버젼 1.3 이후로는 solexa score 대신에 phred score를 쓰기 시작한다. 그러나 다른 점이 64를 offset으로 한다는 것. phred score로 범위가 0-62 값을 갖을 수 있으나 현재는 0-40 까지의 값만 사용한다.OBF project에서는 이를 'fastq-illumina'라고 한다.

Fast and SNP-tolerant detection of complex variants and splicing in short reads

GSNAP 논문. 아.. BWA-SW도 봐야 하는데.. 뭐 읽어도 읽어도 끝도 없고 모르는게 너무 많다는 생각밖에 안들고.. 에라이..


요즘 short read aligner는 몇가지 고려해야 할 사항이 있는데 speed와 sequence variant 그리고 splicing event. speed는 suffix tree와 Burrows-Wheeler Transform을 활용한 방법들이 많이 나왔다. sequence variant의 경우 SNP이 1000bp 당 하나가 있고 또한 human polymorphisms의 7~8%가 indel이며 이 coding indel 중 25%는 3nt보다 길단다. 이 같은 sequence variant는 read가 길어짐에 따라 더 심각해진다. splicing event를 찾는 방법으로 exon-exon을 이어서 인공적인 sequence를 만들어서 mapping 하는것이 한 방법이 될 수 있다. 아니면 tophat 처럼 exon 주변의 splice site junction을 찾는것. 그러나 이것들은 exon 정보를 미리 알고 있거나 아니면 expression 이 많이 일어나는 exon에만 적용이 가능하다는 한계점이 있다.
뭐 이와 같은 문제점을 고려해서 만든것이 GSNAP(Genomic Short-read Nucleotide Alignment Program). 아래 그림이 GSNAP이 찾을 수 있는 complex variant 의 예. 또한 GSNAP은 single reference sequence 뿐 아니라 dbSNP 같은걸 포함하는 reference, 여기서 표현하는 것을 빌리자면 'space' reference를 이용할 수 있다고 한다.

Overview
alignment는 search problem과 같다고 보고 searching은 generating, filtering, verifying을 포함한다. efficiency는 generating과 filtering에 의존적이다. MAQ과 같은 기존의 프로그램은 read를 먼저 pre-processing 하고 나서 이 read index를 genome에 대해 generating과 filtering 해서 candidate genomic region을 찾는다. genome이 큰 경우에는 genome을 먼저 preprocessing 하는 것이 보다 효율적이다.

Extensive genomic and transcriptional diversity identified through massively parallel DNA and RNA sequencing of eighteen Korean individuals

GMI에서 나온 논문이다. 완전 부럽다. GMI이고 싶을 뿐이다. 뭐 상황이 그렇지 않으니 할 수 없는거고.. 
여튼 음.. 힘은 빠진다만.. RNA-Seq 분석에 추가할 사항들과 분석 방법들을 뽑아내고자 선택한 논문. 내가 여기서 봐야 할 것들은 SNP detection, indel detection, annotated SNP, alternative splicing, gene fusion에 관한 것. 정확하게 그것들을 어떻게 수행했는지 파악하는데 목적을 둔다.

이 논문에서 10명의 사람의 whole genome sequencing을 하였고 추가적으로 8명의 사람의 exome sequencing을 하였다. 그리고 이 18명의 사람 중 17명의 사람의 transcriptome sequencing을 하였다.
대략적인 개요는 아래 그림과 같다.

SNP and short indel identification
10명에 대한 whole genome sequencing 정보는 아래 표와 같다.

solexa read의 경우 GSNAP으로 SOLID 데이터의 경우 Bioscope로 hg18에 mapping.
SNP detection을 위한 방법은 예전 논문을 따른단다. 사실 이 논문을 제대로 본적이 없어서 여기서 정리한다. 아래 reference를 참조
Rare and population-specific variants
Large deletions with breakpoints
Transcriptome sequencing analysis
Comparison of DNA and RNA sequence
New sequences from de novo assembly



---------------------------checklist------------------------------
1.Sequence Alignment :  GSNAP을 이용해서 align. 5% mismatch까지 허용해서 highest scoring alignment를 선택한다. 200bp 까지의 read는 GSNAP 사용, 그 이상은 GMAP을 사용하길 권장. 자세한 GSNAP의 내용은 여기 참조.
2.SNP detection : korean genome 논문에 보면 Alpheus software system으로 SNPs랑 indel을 detection 했다고 나오는데.. 아무래도 이거 상용인거 같다는 생각이 든다. 여튼 SNPs call이 된 것들을
autosome의 SNP의 경우 4개 이상의 unique read가 있고, 20% or higher aligned reads 그러니까 mapping된 read 중 SNP으로 나온 read의 비율이 20%로 이상이 될때 이를 SNP로 보고 그 비율이 90%가 넘어가면 이를 homozygous SNP라고 여긴다.
3.indel detection
3.annotation of SNP
4.alternative splicing
5.gene fusion

Optimization of de novo transcriptome assembly from NGS data

제목이 상당히 인상적인라서 함봐야 할거 같단 생각에 선택한 논문.
요즘 NGS data assemblers는 대부분 De bruijn graph 방식으로 assembly 하는데 선택 옵션 중의 하나가 k-mer size. 과연 어떠한 k-mer size를 정해야 하는가에 대해 고민하게 된다.


velvet의 경우 1.palindrome을 피하기 위해 홀수로 2.당연한 것이지만 read length보다는 작게해야 한다고 한다. 이 k-mer의 size는 specificity와 sensitivity의 trade-off 라고, 즉 k-mer의 길이가 길어지면 specificity가 높아지는 반면 sensitivity가 떨어진다. 이들이 경험적인 해결책으로 내놓은 것은 다음 공식을 따른다.
C_k = C*(L-k+1)/L
여기서 C는 standard coverage, C_k는 k-mer coverage, L은 read length, k는 k-mer size로 경험상 C_k가 최소 10은 되야 하고 20이 넘어가면 이는 "wasting" coverage가 된다고 한다(위 공식에서 풀어서 생각하면 일반적인 coverage를 구하기 위해 read count * read length / genome size를 하는데 read length 대신에 read에서 나올 수 있는 k-mer의 수를 곱해준 것으로 각 base에서 몇 개의 k-mer가 cover되는가, 즉 말그대로 k-mer coverage를 구하게 되는 것이다). 


abstract를 보자면..
de novo assembly of transcriptome의 최적화를 위해 여기서는 두가지 방법을 제시한다. 한가지가 Multiple-k method, 다른 한가지가 STM(Scaffolding using Translation Mapping) method. Multiple-k method는 말그대로 k-mer size 다양하게 해서 assembly 한다는 것이고 STM method는 가장 유사하다고 여겨지는 이용가능한 reference proteome을 이용해서 같은 protein에 mapping되는 contig들을 scaffolding 하는것. 이 두 방법을 사용해서 catfish, 매기에 transcriptome 을 분석하였단다.


음.. 굉장히 별거 아닌거 같이 보인다.. k-mer 의 길이를 다양하게 하는 것은 velvet등의 assembler에서의 권장 사항이며 reference based assembly 역시 velvet에서 columbus라는 모듈로 나온것인데.. 뭐 끝까지 함 봐보자.


intro가 상당히 긴데 기억할 만한 것은 higher expressed 되는 transcript의 경우 k-mer 길이를 길게 함으로 해서 더욱 긴 contig를 만들수 있는데 반해 poorly expressed 되는 transcript경우에는 짧은 k-mer size가 났단다. 결국 이 k-mer의 길이의 선택은 어느것에 초점을 맞추느냐인 주관적인 문제라고. 음.. abstract 만 보고 reference based assembly에 대해 뻔하고 했는데.. 여기서도 다양한 선례를 보인다. 다만 선례들은 close relative genome을 이용했으므로 제한적이라는 것. genome을 사용했다는 건 reference로 사용한 genome의 evolutionary distance가  멀수록 nucleotide difference가 클 수 밖에 없으므로. 이 문제의 해결 방법으로 amino acid sequence를 이용한다는 것.


아하. 이 논문의 차별성은 보통은 다양한 k-mer 길이를 이용하여 테스트를 해보기는 하나 결국은 하나의 k-mer size를 정해서 분석하는데 반해 이 논문은 다양한 k-mer에서 나온 결과를 전부 이용한 다는 것이며, 보통 reference based assembly에서 genomics sequence가 사용되는데 반해 여기서는 amino acide sequence를 사용한다는 것.


아래 두 표에서 보듯이 k-mer size를 늘려가면 확실히 sensitivity는 줄어드나(table1에서 reference coverage 감소) contig의 average coverage는 늘어남을 알 수 있다(table2의 rpkm 값 증가). 여기서 집고 넘어가는 것 중에 하나가 k-mer size가 늘어날 수록 rpkm의 mean은 커지고 동시에 SD도 커진다는 것. 즉 k-mer size가 커지면 적은 수의 다양한 expression level을 보이는 transcript의 assembly가 된다는것. 곧 k-mer size마다 assembly 결과 특징이 다르다.
아래표는 모기의 trascriptome 데이터 가지고 evaluation을 해본 것이다.

그래서 이 논문에서는 다양한 k-mer size를 이용하는 두가지 방법을 제시. 하나가 subtractive Multiple-k 와 additive Multiple-k.
subtractive Multiple-k : 큰 k-mer size로 assembly 하고 나서 nonassembled read만 모아서 작은 k-mer size로 assembly
additive Multiple-k : subtractive Multiple-k와는 달리 큰 k-mer size로 assemble 했을때 contig 생성에 참여한 read를 제거하지 않고 다시 중복해서 작은 k-mer size로의 assembly에 사용. 이렇게 되면 contig에 redundancy가 생기게 되는데 CD-HIT-EST를 이용해서 clustering을 한뒤 가장 긴 contig만 남긴다.
위 표에서 볼수 있듯이 subtractive Multiple-k는 그다지 나아진 점이 없는걸로 보인다. 반면 additive Multiple-k의 결과는 상당히 개선됨을 보인다. 이 논문에서 찝어서 이야기 하는 것중 하나가 No. of transcript로 보면 Multiple-k가 19 k-mer보다 나은게 없어보이지만 100이상인 contig의 수가 2배 이상이라는 것(커버하는 transcript의 수는 비슷할지라도 그 transcript의 cover 면적은 Multiple-k가 훨씬 좋다는 것. 그렇다면 cover 했다는 definition이 뭐냐? blastn을 돌렸을때 query의 95%가 align에 참여하고 reference의 99% identity가 있을 경우를 의미 한다. 곧 query 대비 길이가 95%가 넘어가면 그 transcript를 커버했다고 보는 것).


STM method : initial assembly 결과 contig를 translation해서 reference proteome에서 orthologous region을 찾는것. 만약 서로 다른 contig의 translated amino acid가 하나의 protein에 mapping된다면 이를 scaffolding 하는 것이다. 이 과정에서 initially unassembled read를 사용하느냐 마느냐에 따라 STM+ 와 STM-로 나뉜다. 자세한 과정은 아래 그림을 참조한다.

아래 표가 simulated read를 STM method을 이용해서 assembly 한 결과. 에러는 적은데 반해  maximum length랑 N50이 커진 것을 알 수 있다. 그리고 이 STM method가 transcript length나 adundance에 bias가 있는지도 확인했는데 그런건 없단다.

마지막으로 STM method의 한계점에 대해서 이야기 하는데.. 어떻게 보면 STM method가 EST로 orthology 찾는거랑 비슷한데 다른 논문에서 말하길 비교하자고 하는 두 종 중 하나의 종의 transcriptome이 완성됐다면 그 정확도가 상당히 높다는 것. 또한 비교하고자 하는 종의 evolutionary distance가 멀어짐에 따라 prediction 할 수 있는 ortholog는 줄어들지만 accuracy는 거의 똑같다고. 곧 이를 STM method에 적용해서 보자면 complete transcriptome을 사용한다면 reference의 evolutionary distance에 상관없이 low error rate의 결과를 얻을 것이라는 것.