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 옵션이 같이 있지 않으면 당연히 무시된다.

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

우분투 10.x에 Crossover 모니터(30Q3 PRO) 설치

설치에 성공했다는 것은 아니고, 이제 해 보려 한다는 뜻이다. 거의 일주일 넘게 blast2go 작업을 하느라 컴퓨터를 재부팅할 기회가 없었다.

현재는 다음과 같이 듀얼 모니터를 쓰고 있다.

LG L226W 1680 x 1050
LG L1954 1280x1024

비디오카드 드라이버는 NVIDIA-Linux-x86_64-295.53. 비디오카드 제품명이 뭐였더라? Lantic 어쩌고 하는 것이었는데.

비상 사태를 대비하여 /etc/X11/xorg.conf를 오늘 날짜로 백업해 두었다.

DVI cable을 연결하고 재부팅. 어라? 아무런 손을 대지 않았는데 알아서 잘 잡히고 2560x1600의 해상도가 나온다. 모니터 명칭은 DND LGi_DIGITAL이라는 것으로 표시된다.

이렇게 쉽게 화면 설정이 되어 버리니 정말 허무하다... 옆에 놓은 19인치 모니터와 배열을 맞추기 위해 nvidia-settings에서 layout만 조금 손을 보았다.

tophat mapping statistics와 관련한 수수께끼

tophat-cufflinks pipeline을 이용하여 RNA-seq data에 대한 분석을 시도하였다. 총 여섯 개의 Illumina data file을 매핑하였는데, 이 중에서 3개 결과가 좀 이상하다. 예를 들어 3400만개의 read가 투입되었는데 mapped read의 수는 이보다 많은 5400만개라는 결과가 나온 것이다.

상식적으로 이해가 되지 않는다. fastq file 안에 들어 있는 read의 수를 아무리 헤아려 봐도 34,165,827개이다. tophat 결과 디렉토리에 들어 있는 정보 파일도 거짓말을 하지 않는다.

$ cat prep_reads.info
min_read_len=40
max_read_len=40
reads_in =34165827
reads_out=34129283

그런데,

$ samtools flagstat accepted_hits.bam
54565324 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
54565324 + 0 mapped (100.00%:-nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

이건 도무지 말이 되지 않는다. spliced alignment를 했으므로, 혹시 하나의 read가 여러 곳에 나뉘어 정렬한 경우 이것이 SAM file 안에서 별도의 line으로 취급되나 싶어서 다시 한 번 점검을 해 보았다. read name을 추출하여 만일 중복이 있다면 이를 제거하고 그 수를 센다.

$ samtools view -F 0x004 accepted_hits.bam | awk '{print $1}' | uniq | wc -l
54565305

으헉! 미치겠군.

2019년 8월 6일 업데이트

당시에는 multiple mapping에 대해서 별로 진지하게 생각해 보지 않았었다. 부끄러운 일이다!

또 LVM의 디스크 오류 발생

LVM을 이용하여 묶어 놓은 다섯 개의 디스크(씨게이트 바라쿠다 2TB) 중 또 하나가 문제를 일으켰다. 동일 HDD에서 벌써 세번째의 오류던가... 다행히 인텔리전트한 OS께서는 입출력 오류가 나는 즉시 문제의 로지컬 볼륨을 read-only 상태로 만들어 두어서 데이터의 손실은 일어나지 않았다. 이 상태에서 pvdisplay를 해 보면 오류가 난 디스크는 보이지 않는다. 아직까지는 재부팅하여 파일시스템을 체크하면 원상태로 복구가 되지만, 반복적으로 오류가 발생하고 있으므로 해당 디스크를 아예 제거하거나 다른 것으로 교체할 필요가 있다. 다행히 비교적 최근에 전체 백업을 해 두었으므로 데이터를 유실할 염려는 없다.


LVM에 물리 볼륨을 추가하는 것은 쉽지만, 제거하는 것은 매우 까다롭다. 볼륨 그룹에 여유 익스텐트가 하나도 없이 논리 볼륨을 꽉 채워놓은 상태인데다가, 디스크 베이조차 비어 있는 것이 없다. 가장 정석적인 방법은 다음과 같을 것으로 생각된다.

1. 새 디스크를 장착하고 물리 볼륨으로 설정한 뒤 볼륨 그룹에 추가한다. 물론 나는 이렇게 하기 어렵다. 빈 베이가 없으므로. 또는 파일 시스템을 줄이고(resize2fs) 논리 볼륨을 줄여서(lvreduce) 빼 낼 디스크 하나 이상의 공간을 확보한다. 이 과정에서 파일이 무사히 남아 있을 것으로는 확신하기 어렵다. 
2. 제거할 디스크에 있는 익스텐트를 여유 공간으로 옮긴다. pvmove 명령을 이용한다.
3. 디스크를 제거한다.

어차피 전체 데이터를 다 백업해 두었으므로, 아예 논리 볼륨을 싹 지워버리고 새로 LVM을 설치하는 것이 더 빠를지도 모른다. 논리 볼륨을 줄이는 과정에서 파일시스템 체크를 하기 때문에 시간이 너무 많이 걸린다. 다음의 방법을 따르면 될 것이다.

1. 마운트 해제.
2. # lvremove -f /dev/dataVG/dataLG
3. # lvdisplay로 논리 볼륨이 없음을 확인.
4. # pvremove /dev/sab1 (물리 볼륨의 제거)
5. # vgremove 볼륨그룹
6. # rm /dev/vg##/group
7. # rmdir /dev/vg##
8. # strings /etc/lvmtab
9. 물리 볼륨 라벨이 제거된 디스크는 파티션 도구를 사용하여 새로 파티션을 한다.

마지막으로는 가장 과격한 방법이다. 디스크에 영구적인 오류가 발생하여 아예 인식조차 되지 않는 경우에는 백업이 되어 있지 않는 경우 해당 디스크의 자료가 유실될 것을 각오하고 시스템을 복구해야 할 것이다. 그렇다면, 전원을 내리고 문제의 디스크를 아예 빼 버린다음 재부팅을 하면 어떻게 될까? 고장이 나서 특정 디스크가 인식이 되지 않는 것과 매우 흡사한 상황이다.  이는 Red Hat Enterprise Linux 6.4 Beta Logical Volumne Manager Administration 자료의 6.6 Removing Lost Physical Volumes from a Volume Group 항목을 참고하면 된다.

[인용] If you lose a physical volume, you can activate the remaining physical volumes in the volume group with the --partial argument of vgchange command. You can remove all the logical volumes that used that physical volume from the volume group with the --removemissing argument of the vgreduce command. It is recommended that you run the vgreduce command with the --test argument to verify what you will be destroying.

BioPerl을 찾아서...

근무 부서를 옮기고, 설치와 유지 보수 때문에 신경을 쓸 필요가 없는 고급(?) 컴퓨터에 계정을 개설하였다. Sun Grid Engine 환경에서 예전과는 비교할 수 없이 빠른 속도로 작업을 수행하게 되니 정말 놀랍기만 하다.

관리자와 일반 사용자가 엄격히 분리되어 있는 환경이라서, 직접 리눅스 컴퓨터를 관리하면서 관리자 권한으로 설치했던 Perl module이나 수많은 응용 프로그램들을 어떻게 다루어야 하는지 새롭게 공부를 해야만 하는 실정이다. 상당한 수의 생물정보학 응용 프로그램이 깔려 있기는 하지만 이것으로 내 수요가 다 충당되지는 않는다. 관리자가 아직 설 연휴에서 돌아오지 않아서 오늘 당장은 물어볼 수 있는 상황은 아니다.

먼저 이 시스템의 @INC 배열을 들여다 보았다.

$ perl -e 'print join "\n", @INC; print "\n"'
/usr/local/lib64/perl5
/usr/local/share/perl5
/usr/lib64/perl5/vendor_perl
/usr/share/perl5/vendor_perl
/usr/lib64/perl5
/usr/share/perl5
.

매우 상식적인 결과이다. perldoc Bio::Seq 실행에서 에러가 나는 것을 보면 BioPerl은 설치되어 있지 않은 것으로 보인다. 그런데 CPAN 자체도 설치가 되어 있지 않다. CPAN shell 자체를 띄울 수 없는 환경인 것이다. 일반 사용자 모드에서 CPAN-1.9800을 직접 설치하려고 하니 참으로 마땅치가 않다! CPAN.pm을 비롯하여 여러개의 모듈이 필요하고, 또 cpan 스크립트도 설치해야 한다.

내부 교육용 자료를 읽어 보면 perl-5.14.1.tar.gz 소스를 받아다가 설치하라고 되어 있다. 그렇다면 홈 디렉토리에 perl interpreter가 이중으로 설치되고 만다. .bash_profile에 홈 계정에 perl interpreter 디렉토리를 PATH로 지정해 놓는 것으로 해결이 될까?

오랜만에 간단한 목공을...

예전에 코팅합판을 구입하여 건반 받침을 만든 일이 있었다. 그 뒤에 철제 스탠드를 사게 되면서 합판을 해체하여 발코니에 쌓아 두었었다. 며칠 전에 건반을 거실로 내놓으면서 스피커를 놓을 공간이 필요하게 되었다. 코팅합판을 잘라서 목공본드로 접합면을 보강한 다음 나사못으로 고정하였다. 공간이 좁아서 스피커의 앞부분은 건반에 걸치게 만들었다.

건반은 Korg X2이고 피아노 음원은 Alesis NanoPiano를 사용 중이다. 액티브 스피커는 Thonet & Vander의 제품이다.

DIY 중에서 톱질만큼 어려운 일이 없다... 20 cm 혹은 그 아상을 손으로 톱질하여 단면이 정확한 직선 및 직각이 되게 하는 것이 과연 가능할까?

코팅합판은 주문 시 단면에 코팅을 해 주므로 재활용을 위해 톱질을 하면 단면이 노출되고 만다. 집에서 하는 공작 수준에서는 사포로 대강 마무리하는 것 말고는 별 방법이 없다. 물이 닿지 않게 조심하는 것 말고는.

합판을 버릴 생각을 하지 말고 재활용을 위해 아이디어를 내 보자!

옴니아팝에 스카이프 설치하기

국외 출장을 갔을 때를 대비하여 옴니아팝에 스카이프를 설치하기 위해 갖은 애를 쓴 끝에 어렵게 성공하였다. 아무리 좋은 앱이라 해도 내 단말기에 설치를 하지 못한다면 무슨 소용이 있겠는가. 더군다나 윈도우 모바일 6.5용의 설치 파일을 찾는 것은 정말 어렵다. 마치 2013년 현재 윈도우 3.1 시대의 응용 프로그램을 찾는 것이나 다를 바가 없다.

스카이프 홈페이지에는 당연히 이런 구닥다리 운영체제를 지원하는 설치 파일이 있을 리가 없다. 검색을 거듭한 결과, 누군가가 알집으로 분할 압축을 해서 올린 것을 발견하였다. 이런 파일을 내려받아 사용하는 것은 상당한 주의가 필요하다. 왜냐하면 이를 올린 사람이 의도했던 아니든 악성 코드 같은 것이 숨어 있을지 모르니까.

옴니아팝과 PC 양쪽에 프로그램을 설치한 뒤 통화를 시도해 보았다. 잘 들린다. 다음으로는 아내의 엑스페리아에도 스카이프 앱을 깐 뒤 통화. 약간 지연이 느껴지지만 불편함은 없는 정도이다. 스카이프가 좋은 점은 스마트폰, 아이패드, 컴퓨터, 일반 전화 등 대상을 가리지 않는다는 것.

요즘은 해외 사용이 가능하도록 언락된 단말기의 경우 해당 국가 통신사의 선불 유심을 구입해서 끼우면 비싼 로밍 요금을 물지 않고도 국내와 통화할 수 있는 것으로 안다. 아마 내 옴니아팝도 해외 유심을 꽂을 수 있는 것으로 안다. 공항에서 삼만원짜리 전화카드를 사던 일이 엊그제 같은데, 환경이 참 많이도 변했다.

카카오톡이나 라인을 쓰면 와이파이 접속이 된 상태에서 무료 음성 통화가 가능하다. 설정 방법도 스카이프보다는 좀 더 간단하다. 그런데 왜 이렇게 수고스럽게 스카이프를 구닥다리 단말기에 까는가? 카카오톡을 쓰지 않는 개인적인 철학에 대해서는 어제 포스팅하였다.