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