만화로 배우는 미분적분

통계학 입문책 저자가 맘에 들어서 그리고 다음 스텝인 수리통계를 대비하기 위해 그냥 미적분이 본지 너무 오래 되서 다시 한번 되새김질을 하기 위해 구매한 책.
생각보다 쉽지 않게 되어 있다. 다른건 모르겠고 정석에서 배웠던거랑은 약간 다르게 공식에 접근한다. 공식만 따라가다 보면 내가 지금 뭐하고 있지란 생각이 들지만 다시 처음과 끝을 보게 되면 신기하게도 말이 되게끔 해놨다.


다른건 모르겠고 확률 밀도 함수와 확률 분포 함수의 관계를 나타낸 섹션이 인상적이였다. 확률 분포 함수를 구할때 확률 밀도 함수를 적분하게 되는게 아무생각 없이 그냥 면적을 더할려고 한다고 생각했는데 두 그래프를 모아놓고 보니 느낌이 오는...


이 책.. 내가 보기엔 정석보다 어렵다. 헐.. 미적분 모르는 사람이 보면 힘들고 아는 사람이 보면 좀 다른 각도로 접근하기 때문에 신선하다. 절대 쉽게 미적분을 배우려는 사람에게는 비추하는 책. 

1000 genome supplementary note

nature에 나온 1000 genome 논문의 supplementary information을 정리한다.


2. Samples
YRI (Yoruba in Ibadan, Nigeria), CEU (ancestry from Northern and Western Europe), CHB (Han Chinese in Beijing, Chian), JPT (Japanese in Tokyo, Japan), LWK (from the Luhya in Webuye, Kenya), TSI (Toscani in Italia), CHD(the Chinese in Metropolitan Denver, CO, USA)


4. Read mapping and generation of BAM files
-quality recalculated -> remap -> merge lanes from the same library (Picard MergeSamFiles) -> remove duplicate (samtools : rmdup for paired end, rmdupse for single end) -> merge libraries to the plaform level -> remove duplicate (Picard MarkDuplicates)


4.1 Reference genome
-NCBI36, revised Cambridge reference sequence instead of mtDNA. sex-specific reference (Y chr only for male, psudoautosomal region masked in Y chr)


4.2 Mapping of Illumina Data
-Maq v0.7 -u -a 1000


4.5 Recalibration of Base Quality Values
-recalibrate qulity after initial alignment. this algorithm(covariate-aware base quality recalibration algorithm) is implemented in GATK software.
http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration 
-effect of recalibaration : the total number of variants called decreased by 2.8%. changing Ti/Tv ratio from 1.07 to 1.96 (true variants around 2. random 0.5)
Ti/Tv ratio : 
http://www.cbs.dtu.dk/staff/dave/roanoke/genetics980415f.htm
http://paup.csit.fsu.edu/paupfaq/paupans.html


4.6 Comparison of Read Data to known HapMap Genotypes
-genotype log likelihood  (samtools pileup -g) was used for matching expected genotype, and if the best genotype did not seperate well from the others(1.2 separation), then removed.  
likelihood:
http://www.aistudy.co.kr/math/likelihood.htm
genotype likelihood : maq paper(http://graphy21.blogspot.com/2011/01/maq.html)


5.SNP calling 
-maq에서 나온 genotype likelihood(GLij(g) = P(Bij,Qij | Gij =g))를 이용해서 snp를 call한다. 
P(Gij = g|Bij,Qij) = P(Bij,Qij | Gij = g) P(Gij=g) / Kij , Kij = Σg P(Bij, Qij| Gij = g) P(Gij = g)  
말로 풀어서 다시 말하면 maq에서 나온 공식으로 genotype likelihood를 구하고 bayesian공식으로 poterior probability, 즉 데이터가 나왔을때 어떤 genotyp이냐를 추즉한다.
이렇게 snp가 call되면  post-processing step으로 false positive를 제거하고 VCF (variant call format) 형식으로 저장한다.
-post-processing filtering {
--expected depth보다 너무 낮거나 높은거(평균 depth의 반 or 두배), 아마도 CNV에 의한 paralog로 잘못 mapping된거라 생각되서
--snp call 부위의 local realignment, indel에 의한 misalignment를 방지 하기 위해(보통 gap open penalty가 mismatch보다 크다)
--poor mapping quality 제거 , reference 자체가 완벽하지 않기때문에 unrepresented region에서 나온 read가 잘못 맵핑될수 있다(경험상 잘못 mapping되는 region을 6.1에 있다). }


5.1 Low-Coverage SNP calling


5.3 Exon project SNP calls 

통계학 입문

요즘 보고 있는 책이다. 우연히 영건씨네 집에서 보고 다시 통계학을 공부할 생각으로다가 2010년 마지막날 회사가 조금 일찍 끝나 시간이 좀 비는 김에 서점에가서 살짝 보고 구매한 책.
너무 다 아는 내용이 아닌가 싶어서 망설였는데, 책의 한 섹션을 보고 망설임 없이 책을 샀다. 그것이 무언가 하면 여러가지 평균에 대한 이야기. 고1수학에서 보면 산술, 기하, 조화 평균이 나오는데 공식만 알고 이런 여러가지 평균을 어따 쓰는지 제대로 생각해 본적이 없다. 사실 6개월 전쯤 이런 생각을 해봤는데 우연찮게도 이 책에서 그런 평균들의 활용을 언급한부분이 있었다. 비록 입문자를 위해서 쉽게 쓰여졌나 개념만큼은 그 어떤책보다 몸으로 깨달게 만들려는 노력이 보이는 책이다. 
1부는 도수분포표부터 구간 추정까지 개념만 잘 정립하기 위해 준비해 놓은 섹션이고 2부에서 조금더 자세하게 설명한다. 괜찮은 책이다. 


다 읽고 나서 생각해보니 무엇보다 맘에 들었던건 카이분포가 정규분포에서 유래 되었다는 것과 t분포가 어떻게 나오게 되었는지이다. 그리고 모집단이 정규분포일 때 표본들의 대표값이 어떤 분포를 따르는지 설명하게 그 분포로 하여금 모집단의 어떠한 대표값의 단서가 없이도 모집단을 추측할수 있는 구간을 추측하게 된다는 식의 전개가 매끄럽게 연결되었다. 


맘에 드는 책이다.

maq

short read alignment의 대표 프로그램, 1000 genome project에서 사용되었다. 사실 이 논문을 볼생각은 없었는데, 1000 genome 논문 정리하다가 
http://genome.gov/Pages/Research/DER/1000GenomesProjectTutorials/DataDescription-GaborMarth.pdf
genotype likelihood라는 용어가 나오는데 이것의 수식이 maq 논문에서 나왔기 때문에 어찌보면 가장 먼저 봐야 할논문이 아닌가 해서 정리한다. 완벽 이해는 아직 한계가 있지만 최대로 끌어 올려본다.



maq 설명


--alignment stage :
1. ungapped match with lowest mismatch score(the sum of  qualities at mismatching bases) in first 28bp(less than 2 mismatch)
2. failed in matching but whose mate pair is mapped are searching gapped alignment
3. maq assigns individual alignment a phred-scaled quality score
4. if read aligned multiple position, then maq peak one randomly and that mapping quality is 0.


--SNP calling stage :
1. produce a consensus genotype sequence inferred from a bayesian statistical model (그리고 이런 genotype consensus에다가 틀릴수 있는 확률을 phred quality로 할당한다)
2. compare consensus genotype sequence with reference to detect SNP
3. filter called SNP in step 2 by predefined rules(for compensation of simplificaton and assumptions used in the statistical model, treating neighbor positions independently)


Method

<--single end read mapping-->
1. build multipe hash tables to index the reads 
2. scan the reference against the hash tables to find the hits(potential read alignment postion)
In More Detail
1. indexing 
-모든 read를 메모리에 올린뒤, read를 indexing (1st template, 2nd template(complementary sequence of 1st template)의 sequence를 24bit integer로 hashing한뒤 이를 interger를 키로 해서 sorting & grouping), 이 정보를 hash table이 갖고 있음, 이런 hash table을 6개 만든다(sequence를 4조각내서 combination for allowance of 2 mismatch).
2. scan reference
-forward, reverse로 base-by-base로 28bp subsequence를 취해서 위의 1step에서와 마찬가지로 indexing한뒤 위 1번 step의 hash tables 에서 hit를 찾는다. 
-hit가 발견되면 28bp(seed) 뒤로 전체 시퀀스를 gap 없이 align했을때 mismatch의 quality를 합한다(q).
-coordinate of hit과 read ID를 또 다른 24-bit integer (h) 로 hash한다. 이는 나중에 mutiple alignment되었을때 같은 최저 q를 갖을때 random하게 position을 고를수 있게 한다.




<--Single end mapping qualities-->
mapping quality(Qs) = -10 log10 Pr{read is wrongly mapped}
Qs가 30 이란건 1/1000확률로 잘못 alignment될 확률을 뜻함
x를 reference, z를 read, u 를 position이라고 할때 
p(z|x,u)는 reference x의 posiotion u에서 read z가 나올 확률인데 이는 단순하게 생각하면 mismatch의 quality의 곱과 같다. 여튼 구하고자 하는 mapping quality,Qs(u|x,z) = -10 log10[1-ps(u|x,z)] 로 p(z|x,u) 를 베이지안으로 이용하면 구할수 있다 (x때문에 헷갈릴수 있는데 x는 reference 시퀀스이니 항상 주어져 있다고 생각하면 그냥 x를 제거해서 생각하면 된다). 그런데 이 베이지안을 구하기 위해서는 p(z|x) 를 구해야 하는데 이는 reference x의 모든 position에서의 p(z,x)의 합을 구해야 하는데 이는 실질적으로 불가능하다. 그래서 나온게 
Qs = min{q2 - q1 - 4.343 logn2, 4+(3-k')(qave -14)-4.343logp1(3-k',28)}
이는 supplementary에 유도되어 있음.
문제는 snp가 mapping시 base error와 같이 취급된다는 것인데, we should set the minimum base error probability as the rate of differences between the reference and the reads(이해한걸론 reference에 비해 실제 sample이 snp를 갖을 비율이 0.001이면, 곧 1000베이스당 1나의 베이스가 snp일 수 있다면 0,001 보다 큰 error확률을 갖는 mismatch base 만을 mapping quality에 사용해야 한다는 것). 그러나 이 역시 approximation이므로, 첫번째 mapping하고 reference를 변형한뒤 다시 mapping해서  mapping quality를 구하길 권장 




<--Paired-end read alignment-->


<--Detecting short indels-->


<--Consensus genotype calling-->
MAQ assumes the sample is diploid.
1.before consensus calling,  MAQ combines mapping quality & base quality (base quality used in SNP calling cannot exceed the mapping quality of the reads) and reasigns the quality as the smaller value between mapping quality & base quality.
2. genotype likelihood
계산시 최대 많이 나온 2 종류의 nucleotide에 대해서만 고려(나머지는 error로 간주)
이 역시 supplementary 참조
genotype likelihood = Pr(D|genotype) , the probability of data given each possible genotype.
3. prior of genotypes
heterozygote의 확률을  r이라고 가정하면 각각의 homozygote는 (1-r)/2 로 한다. known SNP에 대해선 r은 0.2로 novel SNP에 대해선 0.001을 설정
4. posterior probability
실질적으로 구해야 하는것 P(g|D), 데이터(reads) 가 주어졌을때 genotype의 확률로 g^ = argmax P(g|D) 를 genotype으로 취하고 그 observing quality를 Qg = 10log10 [1 - P(g^|D)] 

위랑은 크게 관련 있다고 봐야 할지는 모르겠지만. 그냥 참조
hash : http://internet512.chonbuk.ac.kr/datastructure/hash/hash2.htm

parallel with python

multi core를 활용해야 한다는 생각에 여러 프로세스를 만들어서 돌릴 생각을 해본다 병렬로다가.. 사실 python에서 thread를 이용하는 것이 multicore 를 이용하는 것인줄 알았다. 그런데 여러 블로그를 보다 보니.. GIL 인가 때문에 한 프로세스 내의 여러 쓰레드는 결국  하나의 cpu를 잡고 돌아가는 것이란걸 알았다. 


IBM의 python 관련 글 (첫번째 목록에 보면 multiprocess에 관한것이 있다)
https://www.ibm.com/developerworks/mydeveloperworks/bookmarks/html?tag=python&base=http://www.ibm.com/developerworks/aix/&dwapp=AIX%20and%20UNIX%20zone&lang=en


IBM에 thread에 관한 것
http://www.ibm.com/developerworks/aix/library/au-threadingpython/index.html


장혜식 님의 블로그(찾기로 thread를 검색)
http://openlook.org:625/blog/tag/python/




http://docs.python.org/library/multiprocessing.html


http://www.parallelpython.com/


http://www.benjaminlog.com/entry/Python-for-Unix-and-Linux-System-Administration-Noah-Gift-Jeremy-Jones

language benchmark

python을 주로 쓰는데 비슷한 작업을 친구가 perl로 짯는데.. 속도가... 너무 차이난다. file IO가 많아서 그런지.. 좌절.. perl도 공부해야 겠다는 생각이 든다. 아래는 언어 performance 비교한거.
http://www.bioinformatics.org/benchmark/results.html

그리고 보통 프로그램 시간 체크 하기 위해서 time 명령어를 쓰는데 default로 real, user, sys 모드의 시간이 나온다. 거기에 대한 설명은 아래 링크
http://stackoverflow.com/questions/556405/what-do-real-user-and-sys-mean-in-the-output-of-time1


그렇다면 두 프로그램의 퍼포먼스를 비교할려면.. user + sys 로만 비교해야 하는건가... 음..


-이번에 이런저런 check하면서 알게된 tip
--어떠한 key값이 사전형 변수에 key로 있는지 없는지 확인하고 대입할경우 예를 들어 

key = 'good'
d_some = {#something in here#}
if key in d_some.keys():
      d_some[key] = "something"

이렇게 하는 것보다

try:
      d_some[key] = "something"
except:
      pass

식으로 하는것이 더 빠르다. 특히 d_some의 key가 많을때.


--그리고 itertools 괜찮은거 같다. 활용해 볼 만하다.
http://www.builderau.com.au/program/python/soa/Faster-smaller-clearer-Python-iterator-tools/0,2000064084,339280043,00.htm
위 내용을 간단하게 정리 하자면 zip()을 써서 두 list 합쳐서 거기서 하나씩 iterative하게 뽑나낼려면 우선은 두 리스트를 합치는 과정이 들어가기 때문에 느린데 izip을 쓰면 합치는 과정 없이 그냥 iterative object를 생성한다는 것 같은데..  그렇다면 고스란히 메모리에 올려놓는건가? 두 리스트를? 음.. 


--python에서 프로세스 시간을 재는 방법
http://www.tutorialspoint.com/python/time_clock.htm

web server

이것저것 고려해보다가.. 결국 웹 서버를 만드는데 옵션 사항들을 결정했다. 웹서버를 만드는데 들어가는 대부분의 사항들이 내가 약한 부분이라 많은 시간을 들여 공부해야 할것이다. 우선적으로 hanil 서버에다가 web server를 구축 할 것이 아니라 google app engine을 사용하기로 한다. python을 이용하며 google app engine에 django가 포함되어 있으므로 장고 식으로 갈거고 이 참에 jquery까지 공부해 보도록 한다.또한 버젼 관리를 위해 전에 해보지 않았던 svn을 이용하기로 한다. (아.. 근데 문제가 있긴 있구나.. 만약 내가 python에서 특정 모듈을 사용할려면 그게 서버에 깔려 있어야 하는데.. 이런 문제는 어떻게 해결하지.. 음.. biopython모듈을 사용해서 사이트를 만든 것이 있는걸로 보아 기본적으로 외부 모듈을 지원하는거 같은데 조금더 알아봐야 한다.예 :http://biosqlweb.appspot.com/ )


-결론적으로 봐야할 list
django, jquery, svn


초보자로서 큰 틀을 다시 그려보면 hanil 서버에 svn 저장소를 생성하고 내 계정에 권한을 줘서 계속 소스를 갱신해서 웹 프로그래밍을 한다. 웹 프로그래밍의 테스트 및 나중에 google app 에 올리는 것은 google_appengine을 이용한다.


-setting & study-
-google app engine SDK 설치
google_appengine_1.4.0.zip 다운로드 하고 path와 PYTHONPATH에 goolge_appengine 위치를 잡아 놓는다.
-svn
http://www.oneone.kr/?document_srl=2030
-html
http://deadfire.hihome.com/html/html01.html
-css 
http://blog.naver.com/rusdudtn2?Redirect=Log&logNo=140047835464
-javascript 
http://deadfire.hihome.com/jscript/jscript01.html (참고 : http://ej5811.blog.me/80094536125 )
-django
http://www.wikidocs.net/read/955 (참고 : http://www.wikidocs.net/read/828)