Tuesday, April 5, 2011

approximate string matching

전에 이야기 했던 당착한 문제점 중 하나. read에서 특정 문자열이 포함된 substring을 찾아서 이를 clipping 하는 일. 사실 바로 전 포스팅, sff 파일 파싱도 그 작업의 일환으로다가 찾아본것. 
그럼 최종적으로 해야 하는 일이 뭐냐. sff 파일에서 MID 별로 분류한 뒤 adaptor 시퀀스를 찾아서 clipping point를 변형해서 MID 별 각각의 sff들을 생성하는것. 
그럼 Roche에서는 이런 프로그램이 없냐? 있단다. 그런데 듣기로 adaptor clipping인가를 5' 쪽만 한다고 했나 뭐래냐. 여튼 필요하단다.

규칙은 간단하다. 
1.MID는 perfect matching으로 찾고 
2.adaptor clipping을 위한 string matching은 n 개의 mismatch (insertion, deletion, mismatch) 를 허용해서 찾는다.

2번 항목을 어떻게 할 것이냐.. 알고리즘에 심하게 취약한 나로서 검색, 이해, 적용을 해보기로 한다.

검색,이해,적용
역시 위키다. 기부금 내길 잘했다. 이런걸 approximate string matching이라고 하는군. edit distance는 몇번 고쳐야 exact match가 되느냐, 고치는 횟수를 의미(여기서 고치기 위한 primitive operation중  swap 이라는 방식은 나에겐 필요 없을 것).
완전 노가다로 찾는 brute-force 방식은 아닌거 같고, 그 다음으로 언급되는 것이 dynamic programming을 이용한 것. 이건 preprocessing이 필요없는 on-line 방식이란다. 그리고 다음에 언급하는 것이 suffix tree. 아.. 역시나.. 진짜 내가 이거 공부하고 만다.

그래서 선택한 것은 bitap algorithm.
요건 levenshtein distance term으로다가 pattern 이 text에 있는지 확인한다. 그러니까 levenshtein distance는 두 string을 줬을때 거리가 얼마인지 그 값(아마도 insertion, deletion, mismatch 갯수를 세는게 아닐까 싶은데) return 하고 ... 아 근데 결정적으로 fuzzy searching의 c 코드를 이해 할수가 없다. 망했다. 


http://en.wikipedia.org/wiki/Levenshtein_distance
이건 그냥 insertion, deletion, substitution까지 고려해서(transposition 까지 포함하고 싶다면 Damerau-Levenshtein distance 참고) 문자열 A와 B의 거리를 판별하는 알고리즘. dynamic programming으로 하는건데 이해는 쉬우나 내 아쉬운 능력으로는 원하는 걸(문자열 A에서 pattern B를 error k 까지 허용했을때 match 되는 substring 모두 찾기) 만들기 위해 변형하기가 힘들다(생각하는 능력이 많이 모자람). 


아 근데 왠걸.. 
http://www-igm.univ-mlv.fr/~lecroq/seqcomp/node3.html
최고다. 요거다. 딱. 한줄기의 빛이 되는 사이트. 아 근데 이거 그냥 Levenshtein distance랑 같다. 다만 text의 -1 index들을 전부 0으로 놔서 시작의 insertion 값을 안줬다는 차이만 있을 뿐. 역시 똑같은 걸 봐도 이해도에 따라 응용도가 달라진다. 다만 이 사이트에서 이해가 안되는건 back trace. 7개의 결과가 나오는데 첫 3개랑 마지막 꺼는 이해가 되는데 가운데 3개는 어떻게 해서 나온건지.. 아.. 이것 참.. 여튼 중요한건 clipping point를 찾는거기 때문에 Levenshtein distance로 나온 matrix의 마지막 행에서 cutoff 보다 작은 최소 값을 찾아서 그 곳을 clipping point로 해주면 될 것같다. 

sff 파일 파싱해보기(?)

바이너리 파일 핸들링 두번째 시간으로(첫번째는 여기) 저번에 어떻게 binary file 핸들링 하는지 대충 감잡았으니까 이번에는 실전으로다가 지금 필요한 일을 하면서 적용해 보려한다.


이번 미션의 목표 : FLX의 sff 파일을 변경하기. 정확하게 이야기 하자면.. sff 라 함은 the Standard Flowgram File의 약자인데 여기에는 GS FLX 기기의 read와 그것의 trace data가 들어 있다. 보통 read를 trim할때 sff 파일을 fastq 파일로다가 전환해서 quality trimming을 하는데, 회사 과장님 말씀으로는 fastq 파일로 newbler로 assembly 한거랑 sff 파일로 assembly 한거랑은 결과가 다르다는데(fastq가 아닌가.. 그냥 fasta 파일로 했을 때 다르다는건가. 음 근데 다르다면 sff 파일에 뭔가 더 들어 있어서 그 정보를 이용한다는 건데.. 음).. 해서 sff 파일의 format을 이해하고 원하는 방향(trimming)으로다가 변형을 하고자 한다.


sff file format 
http://sequence.otago.ac.nz/download/GS_FLX_Software_Manual.pdf page 528부터
여기에 따르면 sff 파일에는 각 homopolymer의 길이에 대한 측정값이 있단다(이건 좀 더 뒷 내용을 봐야 할 것 같다).
일단 sff 파일은 3 section으로 구성
1. common header section : 파일에 한번 나옴, byte size는 31 + number_of_flows_per_read + key_length보다 크거나 같은 가장 가까운 8의 배수. number_of_flows_per_read 가 800이고 flow_char를 보면 TACG가 200번씩 있는걸 보아하니.. 이 flow 라는게 read 그러니까 emulsion PCR이 된 bead를 시퀀싱 할때 TACG 순서대로 각각 200 번씩 nucleotide를 보내줬다는 의미. 
2. read header section : byte size는 16 + name_length 보다 크거나 같은 최소의 8배수, 이 section에서 나타내는 position 숫자는 1-based indexing.
3. read data section : byte size는 2*number_of_flows + 3*number_of_bases 보다 크거나 같은 최소의 8배수, number_of_flows는 common header section에 number_of_bases는 read header section에 있다.


% uint8_t, uint16_t, uint32_t, uint64_t 는 각 각 1,2,4,8 byte numeric value를 뜻함, big endian byteorder, character field는 sinlge-byte ASCII character임, 모든 section은 "eight_byte_padding" 이라는 field로 끝난다(이 필드는 0 to 7 bytes of padding, 이게 뭔소리냐, 각 section의 byte가 8 배수가 되게끔 8배수가 되기에 모자라는 byte 만큼 0 으로 padding이 되어 있다.. 뭐 이런 소리). 나머지는 pdf의 표 참조.


보니까 각 section별로 byte size 만큼 읽어 들여서 필요한 정보만 빼오는데.. 내가 원하는게 read clipping(manual 보니까 trimming이란 용어보다 clipping이란 용어를 사용하는거 같다) 이므로 건드려야 할 부분이  read header section의 clip_qual_left, clip_qual_right, clip_adapter_left, clip_adapter_right 부분이다.