[WMS tech] 2. QC와 트리밍

출처

참조한 글은 문장의 뒤에 링크를 포함한 출처를 표기해두었습니다.

---

0. QC와 트리밍 개요

데이터 분석에서 데이터 자체의 품질, 데이터가 생산되기까지 거치는 과정, 그리고 메타데이터 등 해당 데이터의 여러 속성들을 면밀하게 그리고 충분히 많은 시간을 들여 살펴보는 것은 추후 분석을 위해서 가장 중요한 과정 중 하나이다. 매우 긴 시간을 들여 확인하고 또 확인해도 부족하다고 생각한다.

QC는 Quality Control의 줄임말로, 이 QC과정을 거쳐 sequencing data 품질을 평가하고, 또한 트리밍 (trimming)이라는 과정을 통해 해당 데이터의 품질을 개선한다.

---

1. sequence data에서 quality란?

언뜻 보기에는 sequence data에 quality란게 어디 있다는 건지 감을 잡기 힘들 수도 있다. 보통 sequence data라 함은 3’-AGCGGCCCC-5’와 같이, 연구자가 PCR (polymerase chain reaction)로 증폭시킨 DNA 라이브러리 내 각각의 리드를, 시퀀싱 기계를 통해 읽어 내고 난 뒤의 데이터를 말한다.

일루미나 사의 NGS sequencing을 예로 들어 더 자세히 설명하자면, 이렇게 PCR로 증폭시킨 각각의 리드들이 담긴 DNA 라이브러리가 있다고 가정해보자. 각각의 리드에 대해서 서열을 읽기 위해, 형광을 붙인 dNTP (deoxynucleotide)를 “흘려” 보낸다. 흘려보낸 뉴클레오타이드가 상보성에 의해 각 리드의 서열에 결합할 때마다 형광이 방출되고, 기계는 그 형광을 읽어 결과적으로 리드별 sequencing data를 기록한다 (NGS sequencing에 대한 더 자세한 영상링크).

이때 다양한 sequencing error가 발생할 수 있고, 기계는 그러한 오류율을 계산하여 읽은 서열의 품질로서 기록한다. 품질이 떨어지는 일부 사례에 대한 관련 링크를 올려두었으니 도움이 되었으면 좋겠다.

이러한 오류율을 “Quality score”라 한다. 즉 염기가 잘못 호출되어 읽힐 확률이다. 다음 표를 보자.

Phred Quality Score	Probability of incorrect base call	Base call accuracy
10	1 in 10	90%
20	1 in 100	99%
30	1 in 1000	99.9%
40	1 in 10000	99.99%
50	1 in 100000	99.999%

여기서 어떤 하나의 염기에 대한 Phred Quality Score가 30이라면, 이는 오류율이 0.1% (정확도 99.9%)라는 것을 의미하고, 흔히 이 Phred Quality Score에 33을 더하여 ASCII code로 변환한 Sanger format을 이용하여 다음에 설명할 FASTQ 형식의 파일에 기록해둔다.

---

2. FASTQ

FASTQ 형식의 파일은 너무도 잘 알려져 있기 때문에, 간단하게 아래 사진으로 정리해두었다.

---

3. FASTQC 및 관련 코드

FASTQ 내 수많은 리드들의 품질을 여러 방면에서 검사하여 report 형식으로 output을 내주는, 가장 유명한 툴이 바로 FASTQC이다.

QC할 파일이 너무 많을 경우를 생각해 본인은 다음과 같은 코드를 작성해두었다.

파이썬의 상세 코드 설명은 추후에 천천히 다루려고 한다. 지금은 우선 주석만 표기해두었다.

현재 폴더 내 모든 파일들을 인식하게 한 뒤, 그 파일 하나하나를 fastqc 명령어 적용하게끔 하는 단순한 코드이다. 이 코드도 멀티프로세싱 관련하여 단점이 있는데, 그건 글의 요지와 맞지 않으므로 일단 나중에.. 하하

import os

def main():
    input_path = "/input_path/"
    output_path = "/output_path/"
    file_list = os.listdir(input_path)
    processing_count = 0  # 프로세스 진행 카운트 기록

    for f in file_list[start_index]:
        processing_count += 1
        print(processing_count)

        # fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam]  [-c contaminant file] seqfile1 .. seqfileN
        os.system(f"nohup fastqc -o {output_path} -f fastq {input_path}/{f}")  

main()

---

4. FASTQC output

명령어를 돌린 후에 나오는 output html 코드를 예시로 가져와봤다. 이해하기 어려울 수 있는 그래프들만 한번 살펴보자.

---

Basic Statistics

Measure	Value
Filename	SRR5579973_20V5_P4_1.fastq
File type	Conventional base calls
Encoding	Sanger / Illumina 1.9
Total Sequences	10560182
Sequences flagged as poor quality	0
Sequence length	151
%GC	45

해당 파일 (SRR5579973_20V5_P4_1.fastq)의 총 시퀀스 수 (total sequence)가 대략 1000만개를 넘는 것을 확인할 수 있고, GC 함량 (%GC, 즉 G와 C의 비율, GC 함량이 높을 수록 DNA 밀도가 높으며 변성이 잘 이루어지지 않는다고 함(참조: 식품의약품안전처))이 45%인 것을 확인할 수 있다. GC 함량의 경우, 추후 식별된 다른 taxon의 GC 비율과 비교하거나, 시퀀스별 bias들을 확인함으로써 quality의 일부 요소로 살펴볼 수 있다.

---

Per base sequence quality

한 파일 내에서 존재하는 모든 리드들에서, 리드 내 포지션별 서열의 품질 점수를 boxplot으로 나타낸 것이다. 구글에 fastqc manual pdf를 검색하면 바로 관련 설명들을 볼 수 있다. 메뉴얼에서 설명하기를, 리드의 후반부로 갈 수록 quality가 떨어지는 것은, 대부분의 플랫폼에서 장비가 실행됨에따라 quality는 저하되므로 주황색 (warning)영역으로 quality가 떨어지는 것을 보는 것이 일반적이라고 한다.

---

Per base sequence content

모든 리드내 포지션에서 염기의 비율이 어느 정도인지를 보여준다. 메뉴얼에 의하면 각 염기의 함량은 게놈에 있는 염기들의 전체 양을 반영해야 하지만 어떤 경우에도 서로 크게 불균형해서는 안되고, 만약 강한 bias가 있을 경우 이는 일반적으로 오염된 시퀀스를 나타낸다고 한다. 만약 어떤 위치에서든 A와 T, 또는 G와 C의 차이가 10% 이상이라면 warning 혹은 failure이 반환된다고 한다.

---

Per sequence GC content

파일 내 리드별 GC 함량 분포가 어떠한지 표시한다. 이후 이를 GC함량의 정규 분포 (normal distribution)와 비교한다.

해당 정규분포는, 지금의 raw data가 아닌 무작위의 라이브러리에서의 GC 함량의 분포를 말하고, 이 정규분포와 우리의 raw data 분포를 서로 비교함으로써 혹시 어떠한 이름 모를 원인으로 인해 본인의 raw data가 해당 정규분포와는 다르게 편향되어 있거나, 쌍봉 분포 (bimodal distribution)을 이루는 등 상당히 좀 다른 분포를 보인다면 해당 데이터 자체에 문제가 있지는 않은지 의심해봐야 한다.

---

Per base N content

시퀀스 장비가 충분하게 확신을 가지지 못할 경우 시퀀스 장비는 일반적으로 sequencing data에 시퀀스를 기록하는 대신 ‘N’ 이라는 문자로 대체하여 기록한다. 따라서 이 비율이 일정 퍼센트 이상 증가하면 시퀀스 장비가 데이터를 충분히 읽지 못했다는 것이다.

---

Sequence Duplication Levels

여기에는 수 많은 리드들 중 sequence가 중복되는, 즉 겹치는 리드들이 얼마만큼의 비율을 차지하는지, 그리고 중복되지 않는 리드들은 얼마나 되는지가 표기된다. 메뉴얼의 핵심 부분만 설명하자면, 이러한 중복 수준이 낮으면 낮을 수록 이러한 sequence data가 라이브러리내 대부분의 리드를 읽은 기록일 가능성이 높다 (coverage가 넓다). 하지만 만약 높은 수준의 중복성이 있을 경우 (그래프를 보며 설명하자면 만 개 이상의 duplication level을 가질 경우), 이것은 일종의 enrichment bias (PCR)을 나타낼 가능성이 높다.

여기서는 중복된 시퀀스가 라이브러리 내 50% 이상일 경우 평가를 failure로 칭하니, 사실상 cutoff가 그리 엄격하지는 않다고 생각한다. 따라서 success가 떴어도 그래프를 보고 연구에 맞게 판단하는 것이 중요하다고 생각한다.

이외에도 adapter content 항목이 있는데, 이것은 그 다음 항목인 MULTIQC 및 Trimming에서 설명하려고 한다.

---

5. MULTIQC

지금까지 설명한 것은, 단 하나의 FASTQ 파일 QC report에 관한 것이었다. 하지만 연구를 수행할 때에는 샘플이 많기 때문에 수많은 FASTQ 파일 (여러 샘플)들을 QC해야한다. 그렇기 때문에 이 수많은 각각의 QC report를 한데 모은 종합 report가 필요하다. 그럴때 주로 쓰는 툴이 MULTIQC이다. 위 제목 링크를 클릭하여 더 자세히 확인해보자.

본인은 설명을 위해 일부 샘플 파일 (대락 1700개)들을 다운로드 받고, 툴을 돌려 결과를 확인했다. MULTIQC는 기본적으로 FASTQC에서 보는 항목과 같기 때문에, 소수 몇몇 항목만 살펴보려 한다.

---

한 개의 FASTQC report만 봤을 때는 전의 그림들처럼 큰 문제 없어보였는데, 다양한 파일들을 한데 모아서보니 quality가 좀 떨어지는 파일 (노란색 및 빨간색)도 있다.

---

또한 비록 모든 파일이 success레벨이더라도, 몇몇 파일에서 어댑터가 일정 비율로 존재하는 것을 확인했다. 어댑터는 NGS 시퀀싱시 구성요소로써, 증폭을 위해 DNA 조각에 부착되어 프라이머와 결합하게 한다 (출처 링크). 즉 실제 분석에는 그다지 쓸 데 없는 서열이다. 제거가 필요하다.

---

6. Quality & Adapter trimming

따라서 이렇게 리드 (read) 내 quality가 안 좋은 서열, adapter 서열들은 제거(trimming)해야한다. 파일 내 리드 하나를 완전히 제거하는 것이 아니다. 리드 내에 있는 서열을 제거하는 것이다.

trimming tool은 trimmomatic, cutadapt, fastp외에도 굉장히 다양하다. 여기서는 fastp tool을 이용했다. 본인에게 있어서 제일 사용하기 편리하고, 다른 사람들도 많이 이용하며, 무엇보다도 본인의 연구에서 가장 적합하다 싶은 tool을 선택해서 사용해야 한다. 관련된 tool 논문을 직접 읽는 것이 가장 좋고, tool끼리 서로 비교하여 benchmarking하는 논문을 읽는 것도 좋다고 생각한다.

본인은 다음과 같은 cutoff를 사용하여 trimming 진행했는데, 이 역시 cutoff 하나하나를 자신의 연구에 맞게 잘 선택하고, 그 이유를 기록하여야 한다. 플랫폼에 따라 존재하는 어댑터 서열도 다 다르기 때문에, 이 또한 고려하자.

Quality Score >= 20  
min read length >= 35bp  
fastp 어댑터 제거 cutoff는 Default 사용  

---

7. Trimming 후 Quality check

Trimming 후에는 다시 Quality check하여 quality가 어느정도 개선되었는지, 어댑터는 잘 제거되었는지 확인해야 한다.

본인의 데이터에서는 일부 failure이 보였던 파일들이 trimming 후 그 quality가 어느정도 개선된 것을 확인할 수 있다.

adapter는 trimming 후 multiQC 시 그래프는 따로 안나오고 모든 adaptor 함량이 0.1% 미만이라는 문자만 표기된 것을 확인했다. Trimming 전 함량에 비해서는 상당히 제거된 것을 확인할 수 있다.

---

이전글: WMS tech 1- 개요
다음글: WMS tech 3- human genome filtering