Mark Duplicates

# -----------------------------------------
# STEP 3: Mark Duplicates and Sort - GATK4
# -----------------------------------------

echo "STEP 3: Mark Duplicates and Sort - GATK4"
gatk MarkDuplicatesSpark \
    -I ${aligned_reads}/ERR015533.paired.sam \
    -O ${aligned_reads}/ERR015533_sorted_dedup_reads.bam

這一步的主要目的是處理重複的讀取並對讀取進行排序。是資料預處理的重要步驟，其實字面意思蠻好懂的，但為什麼會有重複序列呢？

重複序列的產生方式：

1. PCR 重複序列：當在PCR反應中進行DNA擴增時，來自同一DNA模板的多個拷貝可能會產生多個相同的DNA片段。這些相同的DNA片段在後續的测序中將被記錄為PCR重複序列。

2. 光學重複序列：在高密度定序流程中，使用illumina定序儀時，因定序反應在flowcell上進行，有時會出現多次讀取相同的DNA片段的情況，這些被多次讀取的DNA片段稱為光學重複序列。

當我們對DNA進行测序時，會生成許多讀取（reads），有時候相同的DNA片段會被多次擴增和测序，這就會導致重複的讀取。這些重複的讀取會在後續的變異體呼叫中引入偏誤，因此我們需要將它們標記並刪除，以確保後續步驟中的準確性和可靠性，經過這一步，我們得到了已標記並排序的BAM文件，其中不再包含重複的讀取。

Base Quality Recalibration

# ----------------------------------
# STEP 4: Base Quality Recalibration
# ----------------------------------

echo "STEP 4: Base Quality Recalibration"

# 1. Build the model
gatk BaseRecalibrator \
    -I ${aligned_reads}/ERR015533_sorted_dedup_reads.bam \
    -R ${ref} \
    --known-sites ${known_sites} \
    -O ${data}/recal_data.table

# 2. Apply the model to adjust the base quality scores
gatk ApplyBQSR \
    -I ${aligned_reads}/ERR015533_sorted_dedup_reads.bam \
    -R ${ref} --bqsr-recal-file ${data}/recal_data.table \
    -O ${aligned_reads}/ERR015533_sorted_dedup_bqsr_reads.bam

這一步的主要目的是校準基質質量分數，以提高變異體呼叫的精確性。
在定序過程中，不同位置的鹼基可能有不同的錯誤率。 為了更準確地呼叫變異體，我們需要調整每個位置的基質質量分數，因此要建立模型，根據訓練變異體，確定每個位置的實際錯誤率，然後應用這些校準資訊來調整基質質量分數，減少假陽性，經過這一步，我們得到了校準後的BAM文件，其中基質質量分數已經調整，以更好地反映每個位置的錯誤率。

為什麼需要 BQSR？

BQSR 的需求主要來自於以下幾個因素：

• 系統性誤差：在测序過程中，可能會引入系統性誤差，例如定序儀器的校準偏差或試劑的效能差異，這些誤差會導致鹼基質量值的不準確性。

• 序列上下文效應：鹼基質量值可能受到其上下文環境的影響，不同的鹼基上下文可能會導致不同的質量分數分佈。

• 變異體檢測需求：在變異體檢測中，我們需要高精度的鹼基質量值，以區分真正的變異體（例如 SNPs）和测序誤差。

PS.感興趣可以再來了解他的算法

Collect Alignment & Insert Size Metrics

# -----------------------------------------------
# STEP 5: Collect Alignment & Insert Size Metrics
# -----------------------------------------------

echo "STEP 5: Collect Alignment & Insert Size Metrics"

gatk CollectAlignmentSummaryMetrics R=${ref} I=${aligned_reads}/ERR015533_sorted_dedup_bqsr_reads.bam O=${aligned_reads}/alignment_metrics.txt

gatk CollectInsertSizeMetrics INPUT=${aligned_reads}/ERR015533_sorted_dedup_bqsr_reads.bam OUTPUT=${aligned_reads}/insert_size_metrics.txt HISTOGRAM_FILE=${aligned_reads}/insert_size_histogram.pdf

這一步的主要目的是收集比對和插入大小的度量信息。
我們需要評估比對的質量以及插入大小的分佈，這將有助於後續變異體呼叫和分析，以便進行質量控制和分析，這些度量信息可以幫助我們評估测序數據的質量，檢測任何潛在的問題，經過這一步，我們得到了比對和插入大小度量的結果，可以用於後續質量控制和分析。

Call Variants - GATK HaplotypeCaller

# ----------------------------------------------
# STEP 6: Call Variants - GATK HaplotypeCaller
# ----------------------------------------------

echo "STEP 6: Call Variants - GATK HaplotypeCaller"

gatk HaplotypeCaller \
    -R ${ref} \
    -I ${aligned_reads}/ERR015533_sorted_dedup_bqsr_reads.bam \
    -O ${results}/raw_variants.vcf

# Extract SNPs & INDELs
gatk SelectVariants \
    -R ${ref} \
    -V ${results}/raw_variants.vcf --select-type SNP \
    -O ${results}/raw_snps.vcf

gatk SelectVariants \
    -R ${ref} \
    -V ${results}/raw_variants.vcf --select-type INDEL \
    -O ${results}/raw_indels.vcf

這一步的主要目的是使用HaplotypeCaller工具進行變異體呼叫，這是變異體呼叫的關鍵步驟，我們需要確定DNA序列中的變異位置，經過這一步，我們得到了原始的變異體呼叫集，其中包括SNP和Indel，用於後續分析。

Reference

A Study on Optimizing MarkDuplicate in Genome Sequencing Pipeline
Legacy GATK Forum