Step 03. Alignment (BWA + Samtools)

트리밍된 페어드엔드 FASTQ를 BWA로 레퍼런스에 정렬하고, samtools로 바로 BAM을 만들어 인덱싱까지 진행합니다. 기본 쓰레드는 20을 사용합니다.

A) Alignment란?

NGS 리드를 레퍼런스 게놈의 위치에 배치하는 과정입니다. 보통 seed-and-extend 전략과 FM-index(Burrows–Wheeler Transform 기반)를 사용해 빠르게 후보 위치를 찾고, 갭/불일치를 허용한 지역 정렬로 최종 매핑을 확정합니다. 결과는 SAM/BAM 형식으로 저장되며, 각 리드에는 MAPQ(매핑 품질)과 CIGAR(매칭/삽입/삭제 패턴) 같은 메타정보가 포함됩니다.

B) BWA란?

BWA (Burrows–Wheeler Aligner)는 BWT/FMI를 활용하는 빠른 리드 매퍼입니다. 로우에러의 짧은~중간 길이 리드(일반적 Illumina)에 최적화되어 있으며, 변이콜링 파이프라인(예: GATK)과 호환성이 매우 좋습니다.

C) BWA 알고리즘 비교 (backtrack / SW / MEM)

정리: 현대 Illumina WGS/Exome는 BWA-MEM이 사실상 표준입니다. backtrack/SW는 레거시/특수 케이스에 한정해 사용하세요.

D) samtools는 무엇을 하나요?

• samtools view: SAM⇄BAM 변환, 태그/리전 필터 (-bS = BAM 출력, SAM 입력 허용).
• samtools sort: (좌표 기준) 정렬된 BAM 생성. -@로 스레드 병렬화.
• samtools index: .bai 인덱스 생성 → 특정 구간만 빠르게 접근.
• samtools flagstat/idxstats: 매핑 요약 통계(QC) 확인.

우리가 쓰는 주요 옵션: -@ <threads>(병렬), view -bS(BAM 출력), sort -o(출력 지정). 필요하면 view -C -T REF로 CRAM 저장도 가능(호환성 확인).

E) 왜 bwa | samtools 스트리밍으로 처리하나?

F) 단일 샘플 예제 (스트리밍 → 정렬 BAM + 인덱스)

Copy
THREADS=20
REF="/home/kang/raw_data/gatk_bundle/v0/Homo_sapiens_assembly38.fasta"
R1="/home/kang/God_Nas/WGS_2nd/Trimmed_fastq/NP-21_R1.trimmed.fastq.gz"
R2="/home/kang/God_Nas/WGS_2nd/Trimmed_fastq/NP-21_R2.trimmed.fastq.gz"
OUT="/home/kang/God_Nas/WGS_2nd/aln"
SAMPLE="NP-21"
mkdir -p "${OUT}"

# Read Group(RG) — GATK 단계에서 필수/중요
RG="@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tLB:${SAMPLE}\tPL:ILLUMINA"

# BWA-MEM 정렬 (SAM을 파일로 쓰지 않고 바로 BAM으로)
bwa mem \
  -t ${THREADS} \
  -M \            # Picard/GATK 호환(보조 정렬을 secondary로 마킹)
  -R "${RG}" \
  "${REF}" "${R1}" "${R2}" \
| samtools view  -@ ${THREADS} -bS - \
| samtools sort  -@ ${THREADS} -o "${OUT}/${SAMPLE}.sorted.bam" -

samtools index -@ ${THREADS} "${OUT}/${SAMPLE}.sorted.bam"

# (선택) QC 요약
samtools flagstat -@ ${THREADS} "${OUT}/${SAMPLE}.sorted.bam" > "${OUT}/${SAMPLE}.flagstat.txt"
samtools idxstats "${OUT}/${SAMPLE}.sorted.bam" > "${OUT}/${SAMPLE}.idxstats.txt"

samtools view -H /path/to/sample.sorted.bam | grep '^@RG'

gatk AddOrReplaceReadGroups \
  -I in.bam -O out.bam \
  -ID FLOWCELL1.L1 -SM SAMPLE1 -LB LIB1 -PL ILLUMINA -PU FLOWCELL1.L1.BC01
samtools index out.bam

G) 배치 실행 스크립트 (run_alignment.sh)

트리밍 디렉토리의 *_R1.trimmed.fastq.gz를 기준으로 페어를 찾아 일괄 정렬합니다.

Copy
#!/usr/bin/env bash
set -euo pipefail

THREADS=20
REF="/home/kang/raw_data/gatk_bundle/v0/Homo_sapiens_assembly38.fasta"
IN="/home/kang/God_Nas/WGS_2nd/Trimmed_fastq"
OUT="/home/kang/God_Nas/WGS_2nd/aln"
mkdir -p "${OUT}"

shopt -s nullglob
for R1 in "${IN}"/*_R1.trimmed.fastq.gz; do
  SAMPLE="$(basename "${R1}" _R1.trimmed.fastq.gz)"
  R2="${IN}/${SAMPLE}_R2.trimmed.fastq.gz"
  [[ -f "${R2}" ]] || { echo "[WARN] ${SAMPLE}: R2 not found — skip"; continue; }

  RG="@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tLB:${SAMPLE}\tPL:ILLUMINA"

  echo "[ALIGN] ${SAMPLE} with BWA-MEM"
  bwa mem -t ${THREADS} -M -R "${RG}" "${REF}" "${R1}" "${R2}" \
    | samtools view -@ ${THREADS} -bS - \
    | samtools sort -@ ${THREADS} -o "${OUT}/${SAMPLE}.sorted.bam" -

  samtools index -@ ${THREADS} "${OUT}/${SAMPLE}.sorted.bam"
  samtools flagstat -@ ${THREADS} "${OUT}/${SAMPLE}.sorted.bam" > "${OUT}/${SAMPLE}.flagstat.txt"
done

echo "All done. BAMs in: ${OUT}"

H) 파일 형식 한눈에

• FASTQ: 원시 리드(염기+품질). 보관/재분석을 위해 원본은 반드시 보존.
• SAM: 텍스트 정렬 포맷(아주 큼). 사람이 보기 쉬우나 실무 저장에는 부적합.
• BAM: SAM의 바이너리 압축. 정렬/인덱스/QC/변이콜링 표준 입력.
• CRAM(선택): 레퍼런스 기반 압축으로 디스크 절약. 뷰어/툴 호환성 확인 필요.