FASTQ file QC 및 전처리기

제12차 KOBIC 차세대 생명정보학 교육의 내용을 되새기기 위해 FASTQ file QC 및 전처리기에 대한 지식을 정리해 보기로 한다. 2일차 강의인 de novo assembly를 열심히 진행해 준 정원형 박사께 감사를...

나는 거의 모든 작업을 CLC Genomics Workbench에서 실시하고 있으므로, 어쩌면 공개형 소프트웨어를 일부러 찾아서 시험해 볼 필요성이 별로 없는지도 모르겠다. 그러나 상호 비교를 통해 부족한 점을 찾아내고 각 프로그램의 특성을 이해하는 것도 바람직하다고 생각한다.

QC용 프로그램은 말 그대로 시퀀싱 결과물이 얼마나 양호한지에 대한 그림이나 보고서를 생성해 주는 역할을 한다. 예전에는 염기 서열 위치에 따른 quality의 box plot을 보여 주는 정도였는데, 요즘은 매우 다양한 수치를 제공한다. 이에 따른 용어도 매우 세련되어졌다. 'per base sequence quality', 'per sequence quality scores', 'per base GC content' 등...

전처리(preprocessing) 프로그램은 de novo assembler나 reference mapping program에 바로 투입할 수 있는 양질의 read를 선별하거나(filtering) 혹은 low-quality/adapter sequence를 제거해 주는(trimming or clipping) 기능을 한다.

내가 가장 먼저 접한 프로그램은 FASTX-Toolkit이다. 설치는 약간 까다로운 편이고(libgtextutils를 미리 설치해야 한다거나), quality encoding을 자동적으로 검출하여 해결해주지 않는다. QC와 관련한 가장 핵심적인 프로그램은 fastx_quality_stat인데, FASTQ file을 입력물로 받아서 텍스트 형태의 리포트를 만든다. 이를 재료로 이용하여 quality box plot을 그리면 된다.

문제는 phred 33으로 인코딩된 최근의 Illumina fastq file을 입력하면 에러가 난다는 것.


$ fastx_quality_stats -i 2962_1.fastq -o report.txt
fastx_quality_stats: Invalid quality score value (char '/' ord 47 quality value -17) on line 4

최신 버전인 0.0.13을 기준으로 보면 소스코드(src/libfastx/fastx_args.c)에서 fastq_ascii_quality_offset의 값을 64에서 33으로 수정하면 된다. 하지만 이렇게 하면 phred 64로 인코딩된 옛날 데이터 파일을 다루지 못하게 된다. 구글링을 통해 알아낸 바에 의하면 fastx_quality_stats 실행 시 -Q 33 혹은 -Q 64 옵션으로 quality encoding을 바꿀 수 있다고 한다. 불행한 현실은 이러한 노하우가 매뉴얼에는 나와있지 않다는 것. 최근의 프로그램들은 인코딩을 자동으로 검출해 주는데 말이다.

KOBIC 강좌를 통해 알게 된 첫번째 프로그램은 FastQC이다. 다양한 플롯을 제공하고, 처리 속도도 빠르며, quality encoding도 자동으로 검출해 준다. 단, 전처리 기능은 없는 것으로 보인다.

KOBIC 강좌를 통해 알게 된 두번째 프로그램은 prinseq이다. html file로 매우 훌륭한 리포트를 만들어 준다. 그러나 사전에 설치해야 할 Perl module이 많고, 더욱 치명적인 현실은 실행 속도가 느려도 너무 느리다는 것.

마지막으로 소개할 것은 내가 가끔 사용하고 있고, 또 작년 5-6차 교육에서 내가 미생물유전체해독 강의를 하면서 실습용으로도 사용한 바 있었던 SolexaQA이다(기억을 더듬어 보니 Fastx-Toolkit도 소개했었다). 이 프로그램은 입력 파일을 포맷을 자동으로 인지하며, DynamicTrim.pl과 LengthSort.pl 스크립트를 통해서 트리밍된 서열을 출력하게 된다.  리포트 생성을 하려면 R과 matrix2png가 설치되어 있어야 한다. 

그러면 SolexaQA를 활용하여 paired end read의 trimming을 실시하는 예를 들어 보겠다. (SolexaQA의 진정한 가치는 tile에 따른 quality를 보여준다는 것인데, 이는 Illumina 장비를 직접 다루는 사람이 아니라면 크게 관심을 보일 내용이 아니므로 여기에서는 건너뛰기로 한다). 안타깝게도 quality trimming 전략에 대해서는 좋은 자료가 매우 적다. quality score라는 것은 마치 언덕을 올라갔다 내려오듯 nucleotide position에 따라 평탄하게 변하는 것이 아니므로, phred score(or quality value) 얼마를 기준으로 자른다고 쉽게 이야기할 수 있는 것이 아니다. 가장 합리적인 quality trimming 전략은 modified Mott algorithm인데, phred documentation에 실린 글이 거의 원전인 것으로 생각된다. 왜냐하면 이보다 더 오래된 문서를 발견하지 못했으므로.

    The modified Mott trimming algorithm, which is used to calculate the
    trimming information for the '-trim_alt' option and the phd files,
    uses base error probabilities calculated from the phred quality
    values. For each base it subtracts the base error probability from an
    error probability cutoff value (0.05 by default, and changed using
    the '-trim_cutoff' option) to form the base score. Then it finds the
    highest scoring segment of the sequence where the segment score is
    the sum of the segment base scores (the score can have non-negative
    values only). The algorithm requires a minimum segment length, which
    is set to 20 bases.

간단히 설명하자면 이러하다. 기준 error(보통 5%, 즉, 0.05)에서 각 염기의 basecalling error를 뺀 뒤 이를 주변부와 더해 나간다. quality가 좋은 영역에서는 error의 합이 주변으로 확장할수록 커지겠지만, quality가 떨어지면 합산한 수치가 점차 줄어든다. 이 값이 최대가 되는 subsequence의 위치를 trimming position으로 삼는다는 것이다. DynamicTrim.pl에서 BWA (-b option) trimming 방법이라고 설명한 것은 바로 이것과 동등한 것이라고 생각된다. 매우 빠른 read aligner인 BWA에서 구현한 trimming algorithm이기에 이렇게 부르는 것 같다.

perl DynamicTrim.pl -b K-BC_1.fastq K-BC_2.fastq 

Info: Using default quality cutoff of P = 0.05 (change with -p or -h flag)

Automatic format detection: Sanger FASTQ format
Info: K-BC_1.fastq.trimmed: mean segment length = 96.0, median segment length = 101
Automatic format detection: Sanger FASTQ format
Info: K-BC_2.fastq.trimmed: mean segment length = 94.0, median segment length = 101

에혀~ 무진장 느리다. 결과가 나온뒤 게시를 하려 했으나, 성급한 나는 미리 게시부터 한다! -b option 없이 기본 트리밍 조건으로 하면 훨씬 빨리 끝난다. 아무튼 두 개의 paired end sequence file로부터 .trim이라는 확장자가 달린 새로운 파일이 생긴다. 이 파일을 열어 보면 서열의 길이가 제각각이고, 또 어떤 것은 헤더만 남고 염기 서열이 전혀 없는 것도 있다. 이제 LengthSort.pl을 사용하여 일정 길이 이상의 것만 필터링하되, pair가 깨진 것과 그렇지 않은 것을 따로 추려내 보자.

perl LengthSort.pl K-BC_1.fastq.trimmed K-BC_2.fastq.trimmed -l 51
(샘플로 사용한 파일이 바뀌었음을 유의할 것. 이 글을 쓰는 현재 -b option으로 시작한 DynamicTrim이 아직도 끝나지 않았기 때문이다)

결과는? XYZ.trimmed.paired1, XYZ.trimmed.paired2, XYZ.trimmed.single, XYZ.trimmed.discard의 4개 파일이 새로 만들어진다. 설명은 필요가 없을 것이다.

<종합 평가>

• 간편한 QC report의 생성은? FastQC
• FASTQ/FASTA 서열의 다양한 조작을 하려면? Fastx-Toolkit
• Trimming
• 그러나 뭐니뭐니해도 CLC Genomics Workbench만큼 간편하고 빠른 솔루션이 없다. 돈만 있으면... 이것이 진리다! FastQC를 제외하면 나머지 프로그램들은 CPU를 달랑 하나만 사용하여 작업을 한다. 그러다 보니 시간이 매우 오래 걸린다. CLC는 올해 업데이트된 버전에서는 graphical QC reporting 기능이 매우 화려해졌다.