bowtie2에서 mismatch 허용하기

Bowtie2에서 mismatch 허용하기.

bowtie2에서는 bowtie와 차이점은 단순히 read length에 따른 최적화나, gap 허용 외에도 

bowtie에서는 최대 3개까지의 mismatch만을 허용하는데 반해 bowtie2에서는 mismatch 또는 indel의 각각의 페널티 점수를 입력하여 read length의 일정 비율만큼 mismatch 허용할 수 있다는 차이점이 있다.

bowtie2에 대한 자세한 옵션은 http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml 를 참고하면 된다. 

이 포스팅에서는 mismatch에 대해서만 다루고자 한다.

For an alignment to be considered "valid" (i.e. "good enough") by Bowtie 2, it must have an alignment score no less than the minimum score threshold. The threshold is configurable and is expressed as a function of the read length. In end-to-end alignment mode, the default minimum score threshold is -0.6 + -0.6 * L, where L is the read length. In local alignment mode, the default minimum score threshold is 20 + 8.0 * ln(L), where L is the read length. This can be configured with the --score-min option. For details on how to set options like --score-min that correspond to functions, see the section on setting function options.

Scoring options

--ma <int>	Sets the match bonus. In --local mode <int> is added to the alignment score for each position where a read character aligns to a reference character and the characters match. Not used in --end-to-end mode. Default: 2.
--mp MX,MN	Sets the maximum (MX) and minimum (MN) mismatch penalties, both integers. A number less than or equal to MX and greater than or equal to MNis subtracted from the alignment score for each position where a read character aligns to a reference character, the characters do not match, and neither is an N. If --ignore-quals is specified, the number subtracted quals MX. Otherwise, the number subtracted is MN + floor( (MX-MN)(MIN(Q, 40.0)/40.0) ) where Q is the Phred quality value. Default: MX = 6, MN = 2.
--np <int>	Sets penalty for positions where the read, reference, or both, contain an ambiguous character such as N. Default: 1.
--rdg <int1>,<int2>	Sets the read gap open (<int1>) and extend (<int2>) penalties. A read gap of length N gets a penalty of <int1> + N * <int2>. Default: 5, 3.
--rfg <int1>,<int2>	Sets the reference gap open (<int1>) and extend (<int2>) penalties. A reference gap of length N gets a penalty of <int1> + N * <int2>. Default: 5, 3.
--score-min <func>	Sets a function governing the minimum alignment score needed for an alignment to be considered "valid" (i.e. good enough to report). This is a function of read length. For instance, specifying L,0,-0.6 sets the minimum-score function f to f(x) = 0 + -0.6 * x, where x is the read length. See also: setting function options. The default in --end-to-end mode is L,-0.6,-0.6 and the default in --local mode is G,20,8.

bowtie2 manual 홈페이지에서 가져온 score 계산 방법이다.

우선 위의 end-to-end score threshold는 -0.6 + -0.6 * readlength 라고 적혀있다.

read length를 100으로 놓으면 비율로 생각하기 쉬우니 그렇게 계산해보면 -60.6 보다 점수가 낮으면 align하지 않겠다는 뜻으로 해석 가능하다.

각 base가 match, mismatch, gap일 때의 score를 살펴보면 마찬가지로 end-to-end일때 mismatch score는 read의 base quality에 따라 점수가 다르게 측정 되는 것으로 보인다. (maximum mismatch penalty와 minimum mismatch penalty가 존재)

preprocess과정에서 low quality base를 자르기도 하고 요즘 sequencing을 하면 대체적으로 high quality read가 많으니 일단 quality가 좋다고 가정할 때의 score인 6으로 계산한다.

(read의 quality가 나쁘면 penalty는 작다. 예를들어 read의 base가 A이며 quliaty score가 낮다면 실제로 이 base는 T이고 mismatch가 아니라 match일 수도 있기 때문이다.)

default 설정의 경우

- mismatch가 10번 생기면 -60점이니 read length의 10%라고 보면 된다.

- gap같은 경우는 open과 extend가 각각 다르게 적용되니 gap이 하나만 있다면 18개 까지 생길 수 있다.

사용할 때는 dafult 값인 mp와 rdg등은 가급적이면 안건드리고 --score-min L,-0.6,-0.3 등으로 바꿔서 (5%의 mismatch를 허용) 해보는 것을 추천 하지만 data 특성에 따라서 mismatch적절한 값을 주고 사용해야 한다.