기본 조건으로 미생물 유전체의 일루미나 시퀀싱 데이터를 227 kb짜리 contig 상에 정렬해 보았다. 아래 그림은 tablet으로 BAM alignment를 시각화한 것이다.
Reference 위에 군데군데 coverage가 비정상적으로 높은(아래로 내려갔으니 깊은?) 부분이 보인다. 자세히 확대를 해 보면 mismatch가 많은 매우 짧은 read를 reference에 억지로 덕지덕지 붙인 듯한 모습이 나타난다. 빨갛게 표시된 염기가 바로 reference와 불일치하는 곳에 해당한다. 아래 그림에서 반투명 팝업으로 나타난 read info를 보자(맨 위의 빨간 박스로 둘러친 read에 해당). Read pair가 서로 짝을 이루었지만 insert의 길이는 겨우 25 bp에 지나지 않는다. 이는 적절한 mapping이 절대로 아니다. k-mer word의 길이를 더 키우면 하면 이러한 mismatched alignment의 생성이 개시조차 되지 않을 것 같다는 생각이 얼핏 든다.
이 염기서열의 정체가 무엇인지는 도대체 알 수가 없다. 어댑터 서열도 아니고, NCBI blast 검색에서도 나타나지 않는다. Mismatch가 이렇게도 많음에도 불구하고 mapping을 해 주는 smalt에게 고마움을 느껴야 하나? 복잡한 변이를 찾으려고 한다면 고마울 수도 있겠다.
일단 시퀀싱 라이브러리 제작 과정에서 나타난 문제라고 여기기로 하고, smalt의 매핑 파라미터를 건드려서 reference와 일치도가 높은 read를 붙여 보도록 한다. Smalt의 매뉴얼 파일(링크)를 찾아서 읽어보니 다음의 두 파라미터가 건드려 보기에 가장 적당해 보인다.
-c mincover Only consider mappings where the k-mer word seeds cover the query read to a minimum extent. If mincover is an integer or floating point value > 1.0, at least this many bases of the read must be covered by k-mer word seeds. If mincover is a floating point value <= 1.0, it specifies the fraction of the query read length that must be covered by k-mer word seeds. This option can be used only when the ’-x’ flag is also set.
-y minid Filters output alignment by a threshold in the number of exactly matching nucleotides. minid is a positive integer or a floating point number <= 1.0 specifying a fraction of the read length.
여기에서 인용한 글귀는 매뉴얼 파일의 것이다. 'smalt map -H' 명령으로 터미널 화면에 표시되는 것과는 조금 다르게 표현되어 있다. 최소 score의 threshold를 설정하는 파라미터도 있지만 복잡하게 생각하기 싫어서 이는 제외하기로 한다.
'-y minid'는 이해하기 쉽다. Read length 전체에 대하여 reference와 일치하는 염기 수 threshold를 fraction으로 지정하는 것이다. 기본값은 놀랍게도 '0'이다.... 그러니 이렇게 대충 붙는 read까지 전부 결과에 남는 것 같다.
'-c mincover'는 이해하기에 약간 까다롭다. 1보다 큰 수가 주어질 경우, read에서 k-mer seed가 뒤덮어야 하는 염기의 수에 대한 threshold를 의미한다. 0과 1사이(1을 포함)의 숫자를 지정할 때에는 fraction에 해당한다. 이 옵션은 '-x'(triggers a more exhaustive search for aligment at the cost of speed)과 같이 쓰여야만 한다.
'-c 3' 조건으로 매핑을 실행해 보았다. 이제는 reference 전체 영역에 걸쳐서 고른 coverage 분포를 보인다.
그러나... 매뉴얼 파일의 제 13 절 'Tuning performance'에 의하면 word length(-k wordlen)와 step size(-s stepsiz)가 sensitivity & accuracy를 결정하는 중요한 요소라 한다. 이 목표를 달성하려면 속도와 메모리 효율은 떨어진다. 하지만 이 두 가지의 파라미터는 매핑 단계가 아니라 인덱스 생성 과정에서 작용한다. 그렇다면 인덱싱을 적절하게 해 두었을 경우 매핑 파라미터는 기본 수치로 두어도 된다는 뜻일까? 내가 지금까지 사용한 인덱스 조건은 '-k 11 -s 1'이었다. 시험적으로 '-k 20 -s 5'로 인덱스를 만든 뒤 매핑을 재실시하였다. Read length가 충분히 길기 때문에(2 x 301 nt cycle) word length를 기본 조건('-k 13')보다 더 크게 해도 된다. 매뉴얼의 설명에 의하면 stepsiz가 더욱 critical한 인자라고 한다.
 |
| 인덱스 작성 방법을 '-k 13 -s 1'에서 '-k 20 -s 5'로 바꾸었더니 mapping parameter에 손을 대지 않아도 양호한 결과가 나왔다. |
(결론) 매핑 파라미터는 손을 대지 않아도 된다. 인덱스 단계를 최적화하는 것이 일치도가 높은 read alignment를 생성하는 지름길이다. 물론 최종 목표가 변이 발굴이라면 불일치를 수용하는 read alignment를 만들어야만 한다!
동일한 alignment score를 갖는 multiple alignment의 처리 방법이라든가('-r seed') partial alignment의 표시('-p') 등에 관한 약간 복잡한 설정도 있으니 매뉴얼을 잘 읽어야 한다. Bowtie2만큼 설정 가능한 파라미터가 장대하지는 않으니 그나마 다행이다.
Insert size의 측정
Read pair overlap을 사용하여 insert size를 산출하는 것은 미처 생각하지 못했다. Reference mapping을 사용하지 않기 때문에 가장 빠르게 계산이 가능하지만, insert가 큰 library molecule의 경우 join이 일어나지 않게 되므로 평균 insert 길이의 결과 값은 실제 크기보다 약간 작은 쪽으로 편향이 일어날 것이다. 실용적으로는 큰 문제가 될 것이라고 생각하지는 않는다.
#1) Via mapping, which requires a reference: $ bbmap.sh in1=r1.fastq in2=r2.fastq ref=ref.fasta ihist=ihist.txt reads=2m pairlen=2000 #2) Via overlap, which requires overlapping reads (they probably overlap given you ran at 2x250): $ bbmerge.sh in1=r1.fastq in2=r2.fastq ihist=ihist.txt reads=2m #3) Via assembly, which requires sufficient read depth and memory to assemble the genome: $ bbmerge-auto.sh in1=r1.fastq in2=r2.fastq ihist=ihist.txt extend2=200
댓글 없음:
댓글 쓰기