BWA + SAMtools比对示例

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam