WGS 분석 실습을 위한 FASTQ 파일 구조 이해하기

 

1. 구조 육안 확인 — 4줄 구조 눈으로 익히기

fastq의 구조 확인을 위한 명령어를 실행함.

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | head -8

 

위와 같은 코드가 출력됨.

@SRR062634.1 HWI-EAS110_103327062:6:1:1092:8469 length=100
GGGTTTTCCTGAAAAAGGGATTCAAGAAAGAAAACTTACATGAGGTGATTGTTTAATGTTGCTACCAAAGAAGAGAGAGTTACCTGCCCATTCACTCAGG
+SRR062634.1 HWI-EAS110_103327062:6:1:1092:8469 length=100
@C'@9:BB:?DCCB5CC?5C=?5@CADC?BDB)B@?-A@=:=:@CC'C>5AA+*+2@@'-?>5-?C=@-??)'>>B?D@?*?A#################

 

첫번째 줄 : @

코드에서 출력된 "@SRR062634.1 HWI-EAS110_103327062:6:1:1092:8469 length=100" 에서 @ 뒤에 있는 정보들이 다 의미가 있음.

 

여기서 :(콜론) 4개인데 이를 쪼개면

HWI-EAS110_103327062 : 6  : 1    : 1092 : 8469
기기명_flowcell ID    : lane : tile : x좌표 : y좌표

 

위와 같이 해석할 수 있음.

 

이를 더 세세하게 뜯어서 해석하면,

• 기기명 HWI-EAS110_103327062 : 시퀀싱한 기계 고유 ID. HWI는 HiSeq/ 이전 장비 prefix.
• lane 6 : 기계 안에 flow cell이 있고, flow cell에 레인이 여러 개 있음. 샘플을 레인에 흘려보내서 시퀀싱을 진행함.
  보통 8개 레인이 있고 여러 샘플을 나눠서 넣음.
• tile 1 : 각 레인을 카메라가 찍는데, 한 번에 다 못 찍으니까 타일 단위로 쪼개서 찍음.
  타일 번호가 몇 번째 사진인지를 나타내.
• x, y 1092:8469 : 그 타일 안에서 이 read가 시퀀싱된 정확한 좌표.
  일루미나는 DNA 조각이 flow cell 표면에 붙어서 빛을 내는 방식이라 위치 정보가 있어.

형식은 실제 값은 어떤 기계로 언제 어디서 시퀀싱했냐에 따라 다 다름.

 

두번째 줄: 염기서열

GGGTTTTCCTGAAAAAGGGATTCAAGAAAGAAAACTTACATGAGGTGATTGTTTAATGTTGCTACCAAAGAAGAGAGAGTTACCTGCCCATTCACTCAGG

 

몇 bp인지 세어보기.

# \n 도 개수로 세는 것을 방지하기 위해 tr -d '\n' 도 추가
gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | sed -n '2p' | tr -d '\n' | wc -c

 

실제 염기서열은 100bp로 Read ID에서 봤던 length=100이랑 일치함.

 

염기 종류별로 몇개인지 카운트하기.

 

그리고 염기 종류별로 몇개인지 카운트해서 A, T, G, C, N 비율로 이 read의 품질을 약식으로 확인 가능함.

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | sed -n '2p' | grep -o . | sort | uniq -c

 

A 34, C 16, G 24, T 26 으로 비율이 균등하지만 N은 없어 정상적인 read로 판단할 수 있음.

 

세번째 줄: +

세번째 줄은 항상 +로 시작하고 줄 1 ID를 반복하거나 그냥 +만 있는 경우도 있음.

 

역할은 줄 2랑 줄 4를 구분하는 구분자로 따로 큰 의미는 없음.

 

네번째 줄: Quality score

@C'@9:BB:?DCCB5CC?5C=?5@CADC?BDB)B@?-A@=:=:@CC'C>5AA+*+2@@'-?>5-?C=@-??)'>>B?D@?*?A#################

 

네번째 줄의 각 문자는 두번째 줄의 각 염기서열에 1:1 대응함.

즉 첫 번째 문자 @가 첫 번째 염기 G의 품질, 두 번째 문자 C가 두 번째 염기 G의 품질임.

 

품질 점수 변환 공식은 밑과 같음.

Phred score = ASCII 값 - 33

 

예시)

@의 ASCII 값이 64이므로, Phred score가 31로 Q31로 계산하여 오류율 0.08%로 꽤 좋은 품질이라고 할 수 있음.

끝에 #은 ASCII 값은 35로, Phred score가 2로 Q2면 오류율 63%로 사실상 쓸 수 없는 품질임.

 

줄 4 끝에 ##### 여러 개가 있었는데 read 끝부분 품질이 급격히 떨어지는 것으로 Trimming이 필요한 이유가 되고, Trimming이 필요하다고 판단할 수 있음.

 

궁금증) 그럼 ASCII문자를 보고 현직자들은 어떻게 분석하지?

실제 분석에서는 이를 외우지 않고, FastQC 돌리면 그래프로 한눈에 보여줘서 직접 계산할 일이 거의 없음.

현직에서는 실제로 밑과 같은 기준만 알고 있으면 됨.

 

• Q30 이상 → 정상
• Q20 이하 → 의심
• Q10 이하 → Trimming 대상

그래서 FastQC 결과 해석이 중요한 것.

지금은 "내부적으로 어떻게 돌아가는지" 이해하기만하고 , 실제로는 툴이 다 해준다는 것을 인지하고 넘어가면 됨.

 

---

 

2. read 수 / 파일 크기 파악 — 데이터 규모 감각 잡기

read 총 개수 세기

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | wc -l

 

4줄이 1개의 read로 나온 숫자를 4로 나누면 400000 / 4 로 총 read의 개수가 100만 read인 것을 알 수 있음.

전 글이였더 데이터 다운로드에서 --maxSpotId 1000000으로 받은 것과 일치함.

줄 1: Read ID
줄 2: 염기서열
줄 3: 구분자
줄 4: Quality score

파일 크기 파악

ls -lh ~/bioinfo/bi_study/wgs_practice/raw/

 

두 파일 다 각 83MB인 것 확인 가능.

 

# 만약 파일 압축을 풀고 싶다면 이 명령어 실행
## 하지만 굳이 안풀어도 분석 가능하니 안풀어도 됨

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | wc -c

---

 

3. paired-end 구조 이해 — _1이랑 _2 관계 확인

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | sed -n '1p'
gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_2.fastq.gz | sed -n '1p'

 

이 두 파일은 같은 DNA 조각을 양쪽 끝에서 읽은 것이기 때문에 Read ID가 똑같음. 

DNA 조각:  ────────────────────────────────
          →→→→→→→→→→   (_1이 이쪽 끝에서 읽음)
                                ←←←←←←←←←←  (_2가 이쪽 끝에서 읽음)

 

그래서 _1이랑 _2는 항상 쌍으로 움직임. read 순서도 1:1로 대응함

 

염기서열 비교

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | sed -n '2p'
gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_2.fastq.gz | sed -n '2p'

 

# 결과 (상보적 역방향(reverse complement))
GGGTTTTCCTGAAAAAGGGATTCAAGAAAGAAAACTTACATGAGGTGATTGTTTAATGTTGCTACCAAAGAAGAGAGAGTTACCTGCCCATTCACTCAGG
CAGAAGAATGTAAACTCCAGGGGGACAGGAGTGCGTGCCGTGTGANTACCGGCTCNCCCTCTGTTCAATGGGGCGATCACCAGCACCTGTAACTTGCCTC

 

같은 DNA 조각인데 서열이 다름.

 

_1은 왼쪽 끝에서 오른쪽 방향으로 읽고, _2는 오른쪽 끝에서 왼쪽 방향으로 읽어서 방향이 반대로 서열이 다르게 나옴.

 

_2 서열에서 GANTACCGGCTCNCC에 있는 N은 기계가 염기를 판독 못했다는 뜻임.

Quality score가 바닥일 때 N으로 표시됨.

---

 

4. 이상한 read 찾기 — 품질 낮은 read, N 포함된 read 직접 뽑아보기

N 있는 염기서열 줄만 뽑아서 확인하기

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | awk 'NR%4==2 && /N/' | head -5

 

N이 몇 개나 있는지 전체 카운트(N이 하나라도 있는 read 개수 찾기)

# N이 하나라도 read에서 확인되면 카운트
gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | awk 'NR%4==2 && /N/' | wc -l

 

100만 read 중 10,397개에 N이 있음.

약 1% 정도로 정상범위에 있다고 함.

번외) N 개수 세기

gzcat ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz | awk 'NR%4==2' | grep -o 'N' | wc -l

 

33474개 나오는 것을 확인.

---

 

5. FastQC 돌리기 — 위에서 눈으로 본 것들이 시각화로 어떻게 나오는지 비교

conda install -c bioconda fastqc

fastqc ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_1.fastq.gz \
       ~/bioinfo/bi_study/wgs_practice/raw/SRR062634_2.fastq.gz \
       -o ~/bioinfo/bi_study/wgs_practice/qc/

---

FASTQC 결과 확인

위 코드를 통해 FastQC를 돌리면 qc 디렉토리에 html 파일이 2개가 생김.

open ~/bioinfo/bi_study/wgs_practice/qc/SRR062634_1_fastqc.html
open ~/bioinfo/bi_study/wgs_practice/qc/SRR062634_2_fastqc.html

 

위 코드를 통해 html 파일을 열면 브라우저를 통해 fastqc 결과를 확인 가능.

 

위 브라우저에서 결과 및 시각화를 확인할 수 있음.

 

Basic Statistics

가장 상위의 Basic Statistics는 위에서 확인했던 것들이 정리되어 있음.

 

• Total Sequences: 1,000,000 — 아까 wc -l 해서 4,000,000 ÷ 4로 확인한 거랑 일치
• Total Bases: 100 Mbp — 100만 read × 100bp
• Sequence length: 100 — Read ID에서 봤던 length=100
• %GC: 40 — 전체 염기 중 G+C 비율
• Sequences flagged as poor quality: 0 — 전체적으로 심각하게 나쁜 read는 없다는 뜻

Per base sequence quality

Per base sequence quality는 각 위치(bp)별로 전체 100만 read의 품질 분포를 박스플롯으로 보여줌.

 

박스 하나가 그 위치의 품질 분포인데, 박스 위아래 폭이 넓을수록 read마다 품질 편차가 크다는 뜻. 파란 선은 평균.

 

패턴을 보면 앞부분은 Q36 근처로 안정적이고, 70bp 넘어가면서 급격히 떨어짐.

100bp 끝에서는 최솟값이 Q2까지 내려오는데 이게 아까 직접 봤던 #####이 그래프로 표현된 것.

 

빨간 X가 뜨는 기준은 어떤 위치의 하위 10% 품질이 Q10 아래로 떨어지면 실패.

이건 일루미나 데이터에서 워낙 흔한 패턴이라 무조건 나쁜 데이터라는 뜻이 아니고, Trimming으로 뒤쪽 잘라내면 해결됨.

 

Per tile sequence quality

 

y축이 타일 번호, x축이 read 위치.

색이 파란색이면 품질 좋은 거고 빨간색이면 특정 타일에서 품질이 나쁜 것.

위에서 Read ID에서 봤던 타일 번호가 있었는데 그 타일별로 품질 이상이 있는지 보는 것.

전체 파란색이니까 특정 타일에 문제 없다는 뜻으로 해석하면 됨.

Per sequence quality scores

read 전체의 평균 품질 분포를 나타내는 그래프.

아까 Per base는 위치별 품질이었고, 이 그래프는 read 하나당 평균 품질이 몇인지 분포로 보여줌.

대부분 read가 Q35~36에 몰려있으니까 전반적으로 품질 좋다는 뜻으로 해석됨.

Per sequence GC content

이 그래프는 각 read의 GC 비율 분포를 보여줌. read 하나당 GC가 몇 % 인지 계산해서 분포로 그림.

 

두 선의 의미는 빨간 선은 실제 데이터를 뜻하고, 파란 선은 이 게놈의 GC 비율로 계산한 이론적 기댓값임.

 

해석 기준은 

• 두 선이 거의 일치 : 정상
• 봉우리가 두 개 : 다른 종의 DNA 오염 의심
• 한쪽으로 치우침 : 라이브러리 준비 과정 문제

그래프 해석

중심이 둘 다 ~40%로 일치하고 모양도 비슷한데 실제 데이터가 살짝 뾰족함.

이 차이로 경고를 나타냈지만 실질적인 문제는 없음.

Basic Statistics에서 봤던 %GC 40이랑도 일치함.

 

Per base sequence content

각 위치별로 A/T/G/C 비율이 얼마나 되는지 보여주는 그래프.

정상이면 4개 선이 평행하게 유지되어야 함.

앞쪽 1~9bp에서 살짝 흔들리는 건 일루미나 특성상 흔한 거고, 이후엔 안정적으로 평행하니까 정상으로 해석 가능.

 

Per base N content

각 위치별 N 비율을 나타내는 그래프.

아까 직접 세어봤던 N이 여기서 시각화됨.

그래프가 0에 딱 붙어있으니까 N이 거의 없다는 뜻이고 정상으로 해석 가능.

Sequence Length Distribution

read 길이 분포를 나타내는 그래프.

100bp에 모든 read가 몰려있으니까 전부 동일한 길이라는 뜻이고 정상임.

Trimming 후에는 이 분포가 달라짐.

Sequence Duplication Levels

 

중복 read 비율을 나타내는 그래프.

똑같은 서열이 몇 번씩 반복되는지 보여주는 것.

대부분 1번만 나오고 2번 이상 중복이 거의 없으니까 정상.

WGS는 중복이 적은 게 정상이고, 중복이 많으면 PCR 과증폭 문제를 의심해야 함.

Overrepresented sequences

특정 서열이 비정상적으로 많이 반복되는지 보는 것.

어댑터나 오염 서열이 있으면 여기 목록으로 나옴.

"No overrepresented sequences"니까 깨끗한 것.

Adapter Content

 

어댑터 오염 여부를 나타내는 그래프.

시퀀싱할 때 DNA에 붙이는 어댑터 서열이 read 안에 남아있는지 보는 것.

그래프가 0에 딱 붙어있으니까 어댑터 오염 없다는 뜻이고 정상.

만약 여기서 튀어오르면 Trimming 때 반드시 제거해야 함.