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