TopHat 옵션 관련한 잡설

글 하나만 올리고 퇴근하자...

TopHat 매뉴얼을 숙독하면서 나름대로 중요하다고 생각되는 사항을 잊어버리기 전에 정리해 본다.

... The software is optimized for reads 75 bp or longer.

괜히 trimmimg해서 길이를 들쑥날쑥으로 만들지 말아라!

TopHat이 어떻게 junction을 찾는가?

TopHat은 유전자 구조에 대한 사전 정보가 없이도 엑손-인트론 경계를 검출할 수 있다. 현재의 NGS 장비는 100 bp 이상의 read를 만들어 내지만 많은 exon은 이보다 작다. 그렇다면 최초의 매핑 과정에서 이를 놓칠 가능성이 많다. TopHat은 입력 read를 잘게 쪼개어서 이를 각각 매핑하고, 최종적으로 이를 다시 종합하면서 end-to-end read alignment를 만든다. 물론 초기 Illumina data와 같이 read length가 ~40 bp에 불과하다면 이렇게 read를 쪼개는 것은 어렵다. 쪼개지는 단위는 기본이 25 bp이고, --segment-length num 옵션으로 변경할 수 있다.

If the read length is greater than or equal to 45bp, we strongly recommend that you decrease --segment-length to about half the read length because TopHat will work better with multiple segments.

TopHat은 다음의 두 가지 증거를 이용하여 splice junction에 대한 정보를 수집하여 일종의 데이터베이스를 만든다.

1. 하나의 read에서 나온 두 segment가 동일 유전체 서열의 서로 다른 부분에 정렬할 때, 그러나 내부 segment는 붙지 않을때. 이는 매우 강력한 정보이다. 그러나 read length < 45 bp라면 적용하기 어렵다. 이 증거에 의해서 GT-AG, GC-AG, AT-AC intron이 찾아질 것이다.
2. 두번째 증거는 소위  coverage island를 이요하는 것이다. 이는 초기 매핑 과정에서 read가 두텁게 쌓여 있는 distinct region을 의미한다. 서로 인접한 island는 splicing에서 서로 만나게 되는 exon으로 볼 수 있다. 이 경우 오직 전형적인 인트론인 GT-AG만 찾을 수 있다. read length가 짧거나(얼마나?) 10만 이하의 read만 있다면, 이 기능을 켜도록 하라.

 어느 증거를 사용할 것인가는 바로 --coverage-search 옵션의 활용에 달렸다. 매뉴얼을 읽어 보아도 확실하지 않은 것은, read의 길이에 따라서 이 옵션이 자동으로 작동하느냐의 여부이다.  그 동안 기록한 노트를 보면, read length가 75 bp 이상이면 자동으로 coverage-search 기능이 꺼지는 것으로 생각된다.

실제로 40 bp RNA-seq를 이용하여 mapping을 실시할 때, 아무런 옵션을 주지 않은 상태에서 다음과 같은 메시지가 나왔었다.

Coverage-search algorithm is turned on, making this step very slow. Please try running TopHat again with the option (--no-coverage-search) if this step takes too much time or memory.

[coverage-search 알고리즘이 작동하므로, 이 과정이 매우 누리다. 만약 너무 시간이 많이 걸리거나 메모리가 많이 소요되어 못참겠다면, --no-coverage-search 옵션을 주고 다시 실행하렴]

정리해 보자.

tophat --no-coverage-search: read가 75 bp도 되지 못하여 coverage-search를 작동하려는데, 마음에 들지 않으면 꺼라!

tophat --coverage-search: read가 충분히 길어서 coverage-search가 작동하지 않는데, sensitivity를 최대한 높이고 싶다면 기능을 활성화시켜라!

하지만 아무리 그래봐야 유전자 annotation 정보가 존재하는 상황이라면 이는 의미가 없으리.

유전자 annotation 정보가 있다면

유전자 annotation은 -g (or --gtf) gtf/gff_file의 형태로 주어진다. 여기에서 공급한 junction에 걸치는 read만을 찾고 싶다면, --no-novel-juncs 옵션을 주어라. -g 옵션이 같이 있지 않으면 당연히 무시된다.

이것 외에도 논할 것이 더 많은데... 일단 이 정도로 하자.