Short reads 의 QC 및 전처리 #

NGS 데이터의 Quality control을 하기 위한 다양한 툴들이 이미 소개되었다. NGSQCToolkit, FastQC등 이 있는데 여기서 퀄리티 스코어에 대해 잠시 알아보고자 한다.

NGS데이터를 분석함에 있어 선행되어야 하는 것이 있는데 바로 시퀀싱 장비에서 생산된 데이터의 QC이다. QC를 통해 Good reads들을 추출하고 이를 분석해 사용해야 한다.

Quality score는 통상적으로 phred scoring scheme를 사용한다. 보통 read quality value = −10 * log10 (error probability) 이렇게 표현을 한다.

예를 들어

Quality = 10 => error rate = 10% => base call has 90% confidence
Quality = 20 => error rate = 1% => base call has 99% confidence
Quality = 30 => error rate = 0.1% => base call has 99.9% confidence

라고 표현한다.

바이오파이썬에서 사용했던 NGS reads들의 데이터를 분석할 때 사용한 예제이다. 여기에서 보면 if min(rec.letter_annotations["phred_quality"]) >= 20) 라는 구문이 있는데 여기서 20이라 함은 Q20을 말하는 것이다. 즉 에러율 1%로 99%의 비율로 신뢰할 수 있다는 말이 된다. 즉 100에 한개 정도는 에러 일 수 있다는 말이다. 따라서 NGS를 이용한 SNP 분석 같은 경우에 보다 고퀄리티의 서열을 가지고 분석을 해야만 좀 더 나은 결과를 기대할 수도 있다는 말이 된다.

from Bio import SeqIO
good_reads = (rec for rec in \
          SeqIO.parse("SRR020192.fastq", "fastq") \
          if min(rec.letter_annotations["phred_quality"]) >= 20)
count = SeqIO.write(good_reads, "good_quality.fastq", "fastq")
print("Saved %i reads" % count)

바이오 파이썬에서도 관련 기능을 제공하지만 FastQC 등을 이용해도 된다.

http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc에서 다운받을 수 있는데 리눅스, 윈도우, 맥을 모두 지원하는 그래픽 기반의 QC 프로그램이다. 사용법은 간단하다. 설치 후에 fastqc만 입력하면 된다.

fastqc s_1.reads.fastq s_2.reads.fastq 
fastqc s_1.read.fastq.gz s_2.reads.fastq.gz

와 같이 간략하게 사용할 수 있다. 결과는 html 로 제공되는데 이미지도 포함하고 있다.

다음으로는 R을 자주 사용하는 사용자라면 R package 중에 ShortRead를 사용할 수도 있다.

library("ShortRead")
sr <- readFastq("s1_reads.fastq")
qa <- qa( sr, lane="character" )
report(qa, dest="my_quality_control")

library("ShortRead")
qaSummary <- qa(".", pattern="*.fastq", type="fastq")
report(qaSummary, dest="my_qc_data")

와 같이 사용할 수 있다.

FastQC 통계 #

Total Reads

grep 'Total Sequences' fastqc_data.txt | awk -F"\t" '{ print $2 }'

GC(%)

grep '%GC' fastqc_data.txt | awk -F"\t" '{ print $2"""%" }'

AT(%) : 100% - GC(%)

Q20 : Per sequence quality scores 이후로 나오는 컨텐츠에서 20 이상인 애들의 count를 전체에서 나누면 Q20 (%)를 구할 수 있음.

Total Count

grep -A36 'Per sequence quality scores' fastqc_data.txt | sed 1,2d | awk '{ sum += $2 } END { print sum }'

Q20 이상 count

grep -A36 'Per sequence quality scores' fastqc_data.txt | sed 1,2d | awk '{ if ($1 > 19) {sum += $2} } END { print sum }'

Q20 이상 % 구하기

echo "scale=2; 97555/99842*100" | bc -l | awk -F"." '{ print $1"""%" }'

Q30 이상 count

grep -A36 'Per sequence quality scores' fastqc_data.txt | sed 1,2d | awk '{ if ($1 > 29) {sum += $2} } END { print sum }'

Q30 이상 % 구하기

echo "scale=2; 71091/99842*100" | bc -l | awk -F"." '{ print $1"""%" }'

마지막으로 NGSQCToolkit 은 관련 페이지를 참조하면 된다.

Incoming Links #

Related Data Sciences (DataScience 0) #

• 파이썬/라이브러리/HTSeq

Related Bioinformaticses (Bioinformatics 1) #

• R (프로그래밍 언어)/Bioconductor/ShortRead
• 서열 기반의 항생제 내성 분석