比对及去重

比对是生物信息分析最常用的方法，就是把自己的数据拿去跟参考数据进行比较，得到差异的结果。而当所要分析的数据是变异时，差异的结果就是突变的碱基，也可以说是突变的基因。

建立索引

首先需要下载一个参考基因组，正如前面所说，我们使用hg38(GRCh38)。hg38基因组也有多个来源，比如UCSC、GIAB、GATK、NCBI、NCI等等。

将fasta(fna)下载下来后，需要建立索引。这里推荐使用bwa-mem2进行比对，它比标准bwa mem快约3倍，但需要建立专用索引。同时也使用samtools建立索引。

# bwa-mem2索引
bwa-mem2 index GRCh38.d1.vd1.fa

# samtools索引
samtools faidx GRCh38.d1.vd1.fa

bwa-mem2的索引文件为.bwt.2bit.64等，如果这些索引文件不存在，流程会自动退回使用标准bwa mem。

比对

推荐使用bwa-mem2进行比对。它需要CPU支持AVX2或AVX512指令集，至少16线程，以及大于16GB的内存。如果条件不满足，会自动退回使用标准bwa mem。

bwa-mem2 mem \
    -t 16 \
    -M -R "@RG\tID:SRR14724513\tSM:SRR14724513\tPL:illumina\tPU:WES" \
    GRCh38.d1.vd1.fa \
    cleandata/SRR14724513_R1.clean.fq.gz \
    cleandata/SRR14724513_R2.clean.fq.gz | \
samtools sort -@ 4 -m 2G -o bam/SRR14724513.sorted.bam -

samtools index bam/SRR14724513.sorted.bam

比对时务必加入 -R 参数指定Read Group，包含SM（样本名）、ID、PL（测序平台）、PU等信息，不然部分识别SM tag的程序可能会报错。-M 参数标记短split reads为secondary，以便与Picard兼容。

如果使用标准bwa mem，命令格式完全一致，只需将bwa-mem2替换为bwa。

标记重复

需要对重复reads进行标记。这里使用sambamba，它的效率比GATK MarkDuplicates更高。

mkdir -p ./tmp

sambamba markdup \
    --nthreads 16 \
    --overflow-list-size 1000000 \
    --hash-table-size 1000000 \
    --compression-level 1 \
    bam/SRR14724513.sorted.bam \
    bam/SRR14724513.markdup.bam \
    --tmpdir=./tmp \
    2> bam/SRR14724513.sambamba.log

samtools index bam/SRR14724513.markdup.bam

文库复杂度评估

标记重复后，使用samtools flagstat统计reads信息，并通过Lander-Waterman模型估算文库大小：

samtools flagstat bam/SRR14724513.markdup.bam > bam/SRR14724513.flagstat

# 从flagstat中提取总数和重复数
TOTAL_READS=$(grep "in total" bam/SRR14724513.flagstat | awk '{print $1}')
DUPLICATES=$(grep "duplicates" bam/SRR14724513.flagstat | awk '{print $1}')

# Lander-Waterman文库复杂度估算
awk -v N="$TOTAL_READS" -v D="$DUPLICATES" '
BEGIN {
    U = N - D;
    if (N == 0 || U == 0) { C = 0; }
    else if (U >= N) { C = N; }
    else {
        lower = U; upper = N * 1000.0;
        for (i = 0; i < 100; i++) {
            C = (lower + upper) / 2.0;
            U_est = C * (1.0 - exp(-N / C));
            diff = U_est - U;
            if (diff < 0) diff = -diff;
            if (diff < 0.1) break;
            else if (U_est < U) lower = C;
            else upper = C;
        }
    }
    printf "Total Reads:          %d\n", N
    printf "Duplicate Reads:      %d\n", D
    printf "Unique Reads:         %d\n", U
    printf "PERCENT_DUPLICATION:  %.4f\n", (N > 0 ? (D / N) : 0)
    printf "--------------------------------------\n"
    printf "ESTIMATED_LIBRARY_SIZE: %d\n", C
}' > bam/SRR14724513.sambamba.stats

ESTIMATED_LIBRARY_SIZE是文库复杂度的估计值，可以理解为文库中理论上独立的DNA分子数量。这个值越大，说明文库多样性越高。

结果质控

在比对后，需要进行多层次的质控。

覆盖度统计

使用xamdst统计靶向区域的覆盖情况、覆盖深度等信息：

xamdst -1 \
    -p bed/S31285117_Covered.bed \
    -o qc/xamdst \
    --threads 16 \
    --cutoffdepth 1000 \
    bam/SRR14724513.markdup.bam

cp qc/xamdst/coverage.report.json qc/SRR14724513.xamdst.report.json

--cutoffdepth 1000指定关注的深度上限为1000X，-1表示单端模式统计。

比对信息统计

使用gatk来统计比对质量和靶向捕获效率：

gatk CollectAlignmentSummaryMetrics \
    -I bam/SRR14724513.markdup.bam \
    -R GRCh38.d1.vd1.fa \
    -O qc/metrics/SRR14724513.metrics.txt

gatk BedToIntervalList \
    -I bed/S31285117_Covered.bed \
    -O S31285117.interval_list \
    -SD GRCh38.d1.vd1.fa

gatk CollectHsMetrics \
    -BI S31285117.interval_list \
    -TI S31285117.interval_list \
    -I bam/SRR14724513.markdup.bam \
    -O qc/metrics/SRR14724513.hs.txt

CollectHsMetrics会计算On-target rate、Fold enrichment、Fold 80 base penalty等WES特有的质控指标。

样本指纹检测

使用Pengelly 24-SNP面板对样本进行身份鉴定，生成指纹文件：

bcftools mpileup \
    -f GRCh38.d1.vd1.fa \
    -R pengelly_snp.txt \
    --annotate FORMAT/AD,FORMAT/DP \
    bam/SRR14724513.markdup.bam | \
bcftools call -m -Ov -o fingerprint.vcf

这个指纹可以用来跨样本身份追踪，防止样本混淆。

性别评估

需要计算数据的性别，然后与临床信息中的性别进行比较，作为数据的一个初步质控。

使用统计SRY基因覆盖reads数的方案，这是最常用的方案：

# GRCh38
samtools view bam/SRR14724513.markdup.bam Y:2786855-2787682 | wc -l

# GRCh37
samtools view bam/SRR14724513.markdup.bam Y:2654896-2655723 | wc -l

SRY基因位于Y染色体上，男性的测序数据应有一定数量的reads覆盖该区域，女性则几乎没有。一般以30 reads为阈值进行区分。

QC汇总报告

将以上所有质控信息汇总为一个JSON报告：

python3 generate_qc_report.py \
    --sample SRR14724513 \
    --output qc/SRR14724513.qc.json \
    --fastp qc/SRR14724513.fastp_stats.json \
    --xamdst qc/SRR14724513.xamdst.report.json \
    --fingerprint qc/SRR14724513.fingerprint.json \
    --metrics qc/metrics/SRR14724513.metrics.txt \
    --hs qc/metrics/SRR14724513.hs.txt \
    --sry qc/SRR14724513.SRY.count.txt \
    --sry-cutoff 30 \
    --sambamba-stats bam/SRR14724513.sambamba.stats

这个汇总报告包含了fastp质控、覆盖度、比对质量、文库复杂度、性别、指纹等所有质控信息，便于统一查看。