FASTQ와 BAM, 명령어로 어떻게 다룰까요? 품질 확인부터 정렬 파일 점검까지

📌 이 글에서 알 수 있는 것

• FASTQ 4줄 구조와 품질 점수(Phred+33) 읽는 법
• FastQC·MultiQC로 무엇을 보고, 어떤 fail은 무시해도 되는지
• fastp 옵션 하나하나의 뜻과 과도한 트리밍의 위험
• samtools view·sort·index·flagstat·markdup·CRAM 변환과 자주 하는 실수

🧬 FASTQ: 리드 하나는 4줄

FASTQ는 시퀀서가 낸 리드를 담는 텍스트 파일입니다. 리드 하나가 4줄입니다.

@read1/1
GATTTGGGGTTCAAAGCAGT
+
IIIIIHHHGG###FFFEE##

• 1줄: @로 시작하는 리드 이름
• 2줄: 염기 서열
• 3줄: + (구분선)
• 4줄: 염기마다 품질 점수 한 글자. 글자의 ASCII 값에서 33을 빼면 Phred 점수입니다(Phred+33)¹. I(73)는 Q40, #(35)는 Q2입니다. Q30은 오류 확률 0.1%입니다.

정렬(alignment) 원리와 변이 호출은 3편에서 다뤘으니, 이 글은 그 앞뒤에서 파일을 다루는 명령어에 집중합니다.

🔍 FastQC·MultiQC: 품질 확인

mkdir -p qc
fastqc --outdir qc -t 4 sample_R1.fastq.gz sample_R2.fastq.gz

• --outdir qc: 보고서를 저장할 폴더. 폴더가 미리 있어야 합니다. 0.12까지는 -o로 줄여 썼지만, 0.13.0 실행 스크립트에는 --outdir만 있습니다².
• -t 4: 사용할 전체 CPU 스레드 수(기본 4). 0.12까지는 ‘동시에 처리할 파일 수’였지만, 0.13.0부터는 파일 하나에 최대 4스레드까지 나눠 씁니다².
• 자주 하는 실수: 폴더를 만들지 않고 실행해 오류가 납니다.

multiqc . -o multiqc_out

• .: 현재 폴더에서 FastQC·fastp·samtools 등이 남긴 로그를 찾아 보고서 하나로 묶습니다. 분석은 돌리지 않으니 마지막 단계 뒤에 다시 돌리세요³.
• -o: 결과 폴더. 기본 파일은 multiqc_report.html과 multiqc_data/입니다.
• 자주 하는 실수: 예전 실행 로그가 섞인 폴더에서 돌리면 이름이 같은 시료는 마지막에 읽힌 결과로 덮어써집니다. 폴더를 나누거나 -d로 폴더 이름을 시료 이름 앞에 붙이세요³.

볼 것은 위치별 염기 품질, 어댑터 비율, 시료 간 GC 분포와 리드 수 차이입니다. FastQC의 pass·warn·fail은 고정 기준값일 뿐이고, 문서도 라이브러리에 따라 특정 모듈의 fail이 흔하다고 설명합니다⁴.

RNA-seq 고유의 QC 지표(rRNA, strandedness 등)는 12편을 보세요.

✂️ fastp: 어댑터와 저품질 리드 정리

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o trim_R1.fastq.gz -O trim_R2.fastq.gz \
      --detect_adapter_for_pe -j sample.fastp.json -h sample.fastp.html -w 8

• -i·-I: read1·read2 입력. -o·-O: read1·read2 출력(소문자가 R1, 대문자가 R2)⁵.
• --detect_adapter_for_pe: paired-end에서 어댑터 서열 자동 탐지를 추가로 켭니다. 기본은 두 리드의 겹침(overlap) 분석만으로 어댑터를 자릅니다⁵.
• -j·-h: JSON·HTML 보고서. MultiQC는 JSON을 읽습니다³.
• -w 8: 작업 스레드 수(기본 3).
• 자주 하는 실수: -o와 -O를 뒤바꾸거나 R2를 빠뜨려 쌍이 깨집니다.

기본 설정에서도 어댑터 제거와 품질 필터(Q15 미만 염기가 40%를 넘거나 15 bp보다 짧은 리드 제거)는 켜져 있고, 리드 끝을 품질로 잘라 내는 --cut_front·--cut_tail·--cut_right는 꺼져 있습니다. NovaSeq·NextSeq 데이터는 polyG 꼬리 자르기가 자동으로 켜집니다⁵. 2023·2025년 개정 논문은 중복 평가·제거 모듈과 일괄 처리 기능을 소개했습니다⁶⁷.

과도한 트리밍을 조심하세요. 자르는 옵션을 겹겹이 켜면 리드가 짧아져 여러 곳에 붙는 리드가 늘고, 리드 길이가 시료마다 달라집니다. README도 --cut_right를 더 공격적인 방식으로 소개하고, 앞쪽 자르기가 중복 제거를 방해한다고 적습니다⁵. 기본값에서 시작해, 보고서에서 어댑터가 남는지 확인한 뒤에만 옵션을 더하세요.

📦 SAM·BAM·CRAM과 FLAG·MAPQ

SAM은 정렬 결과를 담는 탭 구분 텍스트, BAM은 그 압축 이진 파일, CRAM은 참조 서열과의 차이만 저장해 더 작게 만든 형식입니다. 정렬 줄마다 QNAME, FLAG, RNAME, POS, MAPQ, CIGAR 등 11개 필수 열이 있고 좌표는 1부터 셉니다¹.

• MAPQ: 정렬 위치가 틀렸을 확률을 Phred로 나타낸 값(255는 정보 없음)입니다¹. 정렬기마다 범위가 달라 BWA-MEM은 최대 60, STAR는 유일하게 붙은 리드에 기본 255를 줍니다(명세의 ‘정보 없음’과 뜻이 다름)⁸⁹. -q 기준은 정렬기에 맞춰 정하세요.
• FLAG: 리드 상태를 비트 합으로 나타냅니다¹.

samtools flags 1028이나 Broad Institute의 Explain SAM Flags 페이지가 숫자를 이름으로 풀어 줍니다¹⁰.

🧰 samtools: 정렬·색인·점검

이하 옵션은 samtools 1.25 매뉴얼 기준입니다¹¹.

bwa index ref.fa
bwa mem -R '@RG\tID:s1\tSM:sample' ref.fa trim_R1.fastq.gz trim_R2.fastq.gz | samtools view -b -o sample.unsorted.bam
samtools sort -@ 8 -o sample.sorted.bam sample.unsorted.bam
samtools index sample.sorted.bam

• bwa index(한 번만) 후 bwa mem으로 트리밍한 FASTQ를 정렬합니다(-R: read group).
• sort: 기본은 좌표(염색체·위치) 순 정렬. -@ 8은 스레드, -o는 출력 파일입니다. -n을 주면 리드 이름순입니다.
• index: .bai 색인을 만들어 영역 질의를 가능하게 합니다. 512 Mbp보다 긴 염색체가 있으면 -c로 CSI 색인을 씁니다.
• 자주 하는 실수: 정렬하지 않았거나 이름순으로 정렬한 BAM에 index를 돌려 실패합니다. 색인은 좌표 정렬 파일만 됩니다.

samtools view -b -q 20 -F 0x904 -o sample.filt.bam sample.sorted.bam chr1:1000000-2000000
samtools view -c -f 4 sample.sorted.bam
samtools view -H sample.sorted.bam | grep '^@SQ' | head -3

• -b: BAM으로 출력. -h는 SAM 출력에 헤더를 포함, -H는 헤더만 봅니다.
• -q 20: MAPQ 20 미만 제외.
• -F 0x904: 4+256+2048 비트 중 하나라도 있으면 제외(비정렬·secondary·supplementary 제거). -f는 반대로 비트를 모두 가진 것만 남깁니다.
• chr1:1000000-2000000: 영역 질의. 좌표 정렬과 색인이 필요합니다.
• -c -f 4: 정렬 안 된 리드 수만 셉니다.
• 자주 하는 실수: 참조 이름이 chr1인지 1인지 확인하지 않고 영역을 넣어 결과가 비어 나옵니다. @SQ 줄로 먼저 확인하세요.

samtools flagstat sample.sorted.bam > sample.flagstat.txt
samtools idxstats sample.sorted.bam > sample.idxstats.txt
samtools coverage -q 20 sample.sorted.bam
samtools depth -a -b targets.bed sample.sorted.bam > depth.txt

• flagstat: 전체·primary·중복·매핑·properly paired 리드 수를 QC 통과/실패로 나눠 셉니다.
• idxstats: 염색체별 길이·매핑 리드·비매핑 리드 수. 성별(chrX·chrY)이나 참조에 넣은 오염 서열(EBV 등) 확인용입니다.
• coverage: 염색체(또는 -r 영역)별 리드 수, 덮인 염기 비율, 평균 깊이·MAPQ. -q는 최소 MAPQ입니다.
• depth -a -b: BED 영역의 위치별 깊이. -a가 없으면 깊이 0인 위치가 빠집니다. 중복·secondary 리드는 기본으로 제외합니다.
• 자주 하는 실수: depth에서 -q(염기 품질)와 -Q(MAPQ)는 coverage와 뜻이 반대입니다.

🔁 중복 표시와 CRAM 변환

samtools collate -@ 8 -o sample.collate.bam sample.unsorted.bam
samtools fixmate -m sample.collate.bam sample.fixmate.bam
samtools sort -@ 8 -o sample.pos.bam sample.fixmate.bam
samtools markdup -s sample.pos.bam sample.markdup.bam

• collate: 같은 이름의 리드 쌍을 붙여 놓습니다(bwa 출력은 이미 묶여 있어 생략 가능).
• fixmate -m: mate 정보와 markdup에 필요한 MC·ms 태그를 채웁니다.
• sort: markdup은 좌표 정렬을 요구합니다.
• markdup: 중복에 FLAG 1024를 붙입니다. -s는 통계 출력, -r은 아예 지웁니다¹². 같은 위치에서 시작하는 앰플리콘 데이터에는 쓰지 않습니다.
• 중복을 표시하려면 앞 절 sort 대신 이 경로를 쓰고, samtools index sample.markdup.bam 후 flagstat·coverage를 이 파일로 다시 돌리세요.
• 자주 하는 실수: 좌표 정렬된 BAM에 바로 fixmate를 돌리는 것입니다. Picard MarkDuplicates는 좌표 정렬·이름별 입력을 받고 METRICS_FILE이 필요합니다¹³.

samtools view -C -T ref.fa -o sample.cram sample.markdup.bam

• -C: CRAM 출력. -T ref.fa: 정렬에 쓴 그 참조 FASTA.
• samtools 1.22부터 기본 출력이 CRAM 3.1이고¹⁴, Illumina 데이터에서 BAM보다 50~70% 작았습니다¹⁵. GATK 4.6·Picard 3.4 이하(htsjdk 4.3 미만)는 3.1을 못 읽으니 -O cram,version=3.0으로 만드세요¹⁶.
• 자주 하는 실수: 참조 파일을 지우거나 다른 판본으로 바꾸는 것입니다. 문서는 참조를 잃으면 리드 서열을 모두 잃는 것과 같을 수 있다고 경고합니다¹⁷.

⚠️ 한계와 주의

• 참조 불일치: chr1/1 표기, GRCh37/GRCh38 판본이 BAM 헤더(@SQ)와 BED·참조 FASTA 사이에서 다르면 조용히 결과가 비거나 틀립니다.
• read group 손실: 정렬할 때 -R로 넣지 않으면 시료 정보가 없어, 여러 BAM을 합친 뒤 시료를 구분할 수 없고 read group을 요구하는 하위 도구(예: GATK)에서 멈춥니다. 나중에는 samtools addreplacerg로 넣습니다¹¹¹⁷.
• 이름순 vs 좌표순: fixmate는 이름별 입력(GO:query·SO:queryname), index·markdup·영역 질의는 좌표순(SO:coordinate)이 필요합니다. @HD 줄로 확인합니다¹.

🏛️ 지금 어디까지 왔나

• 도구 판본(2026-11-06): FastQC 0.13.0, MultiQC 1.35, fastp 1.4.0, samtools 1.25, Picard 3.5.0입니다.
• 형식 표준: SAM·BAM·CRAM 명세는 hts-specs 저장소에서 관리되며, 이 글은 2025-08-12 판 SAM 명세를 확인했습니다¹.

🧪 연구자 노트

• FastQC 결과로 어댑터와 품질 문제를 먼저 확인한 뒤 다음 단계로 넘어가세요. 그래야 다시 돌리는 시간을 아낄 수 있습니다.
• BAM은 정렬률만 보지 말고 중복률과 삽입 크기도 함께 보세요. 정렬률이 높아도 라이브러리 문제가 숨어 있을 수 있습니다.

🙋 그래서 나에게는? — 실습 전 체크리스트

• 도구 버전을 --version으로 기록했다(FastQC 0.13은 -o 대신 --outdir)
• MultiQC로 시료 전체의 품질·어댑터·리드 수를 한 번에 봤다
• FastQC fail이 라이브러리 종류 때문인지 먼저 따졌다
• fastp는 기본값으로 시작했고, 자르기 옵션을 더한 이유를 기록했다
• BAM 헤더 @SQ 이름·판본이 참조 FASTA·BED와 같다
• 좌표 정렬 후 색인했고, flagstat·idxstats를 시료별로 남겼다
• 정렬기에 맞는 MAPQ 기준과 FLAG 필터를 정했다
• read group(ID·SM)이 들어 있다
• CRAM으로 바꿨다면 같은 판본의 참조 FASTA를 함께 보관한다

⚠️ 이 글은 의학적 조언이 아니며, 특정 제품을 권하지 않습니다. 실제 옵션은 사용하는 도구 버전의 매뉴얼과 기관 SOP, 자체 검증 결과를 우선하세요.

❓ 자주 묻는 질문

Q. FastQC에서 빨간 fail이 여러 개 나왔는데 다시 시퀀싱해야 하나요? A. 먼저 라이브러리 종류를 보세요. RNA-seq의 염기 조성이나 중복 fail은 흔한 정상 패턴이고, 문제는 품질 급락, 남은 어댑터, 다른 시료와 크게 다른 분포입니다.

Q. 중복 리드는 지워야 하나요, 표시만 해야 하나요? A. 보통은 FLAG 1024로 표시만 하고, 하위 도구가 알아서 제외하게 둡니다. RNA-seq은 고발현 유전자에서 중복이 자연히 생겨 보통 지우지 않습니다(12편).

Q. BAM을 CRAM으로 바꿔도 되나요? A. 저장 공간은 크게 줄지만 같은 참조 FASTA가 영구히 필요합니다. 참조를 함께 보관할 수 있을 때만 바꾸세요.

📚 오믹스 기술 시리즈

• 1편 WGS 플랫폼 비교
• 2편 WGS 커버리지·깊이 설계
• 3편 WGS 데이터 처리와 QC
• 4편 WGS vs WES vs 패널
• 5편 리눅스 명령줄 기본기
• 6편 분석 환경 만들기
• 7편 FASTQ·BAM 실전 (이 글)
• 8편 VCF·BED 실전
• 9편 워크플로 관리
• 10편 R과 Python
• 11편 RNA-seq 실험 설계
• 12편 RNA-seq 데이터 분석
• 13편 small RNA-seq 실험 설계
• 14편 small RNA-seq 데이터 분석

참고