使用一個數(shù)據(jù)樣本（EXON_FYY_1.fastq.gz，EXON_FYY_2.fastq.gz）作為示例，樣本多可以寫成循環(huán)。

Step 0：創(chuàng)建目錄

原始測序數(shù)據(jù)放在0.rawdata中，參考文件放在0.ref_38

cd ~/path/exon
mkdir tools 0.rawdata 0.ref_38 1.QC 2.clean_QC 3.bwa_index_38 4.align 5.stat 6.gatk 7.Germline_Mutations 8.Mutect2
mkdir 1.QC/fast_QC  1.QC/multi_QC
mkdir 2.clean_QC/fast_QC  2.clean_QC/multi_QC  2.clean_QC/trim_galore_QC
mkdir 6.gatk/table
mkdir 7.Germline_Mutations/gvcf/
mkdir 8.Mutect2/annotation 8.Mutect2/maf 

Step 1：質(zhì)量控制

利用工具 fastqc、multiqc進行前期質(zhì)量查看，確保測序數(shù)據(jù)基本合格。
html 文件是質(zhì)控報告，zip 文件是 html 文件圖表對應(yīng)的類。樣本多的話可以用 multiqc 對所有樣本的結(jié)果整合。

fastqc ./0.rawdata/*.fastq.gz --outdir ./1.QC/fast_QC --threads 60 --quiet
ls -lh 1.QC/fast_QC
multiqc ./1.QC/fast_QC/*.zip -o ./1.QC/multi_QC
ls -lh 1.QC/multi_QC

Step 2：trim-galore去接頭

如果從質(zhì)控報告看到一些 reads 的質(zhì)量值較差或測到了接頭，需要去除掉這些質(zhì)量較差的 reads 或接頭。

module load /opt/modules/5.0.1/modulefiles/trim-galore/0.6.10
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --path_to_cutadapt /usr/bin/cutadapt --cores 24 \
-o ./2.clean_QC/trim_galore_QC ./0.rawdata/EXON_FYY_1.fastq.gz ./0.rawdata/EXON_FYY_2.fastq.gz \
>> ./2.clean_QC/trim_galore_QC/EXON_FYY_trim.log 2>&1
ls -lh 2.clean_QC/trim_galore_QC

# 再質(zhì)控：
fastqc ./2.clean_QC/trim_galore_QC/*.fq.gz --outdir ./2.clean_QC/fast_QC --threads 60 --quiet
multiqc ./2.clean_QC/fast_QC/*.zip -o ./2.clean_QC/multi_QC
ls -lh 2.clean_QC/multi_QC

*_val_1.fq.gz 表示經(jīng)過質(zhì)量修剪后的有效數(shù)據(jù)，即去除了低質(zhì)量的部分。

Step 3: 構(gòu)建參考基因組的 BWA 索引

這步驟時間較長可以提前運行。生成的 BWA 索引文件說明：
hg38.fasta.64.amb: 存儲了參考基因組序列的一些元信息，如序列長度、名稱等。
hg38.fasta.64.ann: 存儲了一些額外的注釋信息，有助于在比對過程中進行一些元數(shù)據(jù)的操作。
hg38.fasta.64.bwt: 這是Burrows-Wheeler變換的數(shù)據(jù)，是BWA索引的核心部分之一。BWT是一種數(shù)據(jù)轉(zhuǎn)換算法，用于將原始DNA序列進行變換，以便在比對時能夠更高效地進行匹配。
hg38.fasta.64.pac: 這是FM索引的數(shù)據(jù)，也是BWA索引的核心部分之一。FM索引是一種數(shù)據(jù)結(jié)構(gòu)，用于實現(xiàn)快速查找和比對DNA序列。
hg38.fasta.64.sa: 這是后綴數(shù)組的數(shù)據(jù)，也是BWA索引的核心部分之一。后綴數(shù)組是一種用于存儲DNA序列的子串信息的數(shù)據(jù)結(jié)構(gòu)，用于在比對時快速查找匹配。

bwa index -p ./3.bwa_index_38/hg38 ./0.ref_38/hg38.fasta
ls -lh 3.bwa_index_38

Step 4: 將reads比對到參考基因組上 bwa mem (1h/file)

sam 文件是純文本格式，占用存儲空間，可通過 samtools 轉(zhuǎn)換成二進制的 bam 文件，記錄了每一條 reads 比對到參考基因組的結(jié)果。

bwa mem -M -t 24 -R "@RG\tID:EXON_FYY\tSM:EXON_FYY\tLB:WXS\tPL:Illummina" ./3.bwa_index_38/hg38.fasta \
    ./2.clean_QC/trim_galore_QC/EXON_FYY_1_val_1.fq.gz \
    ./2.clean_QC/trim_galore_QC/EXON_FYY_2_val_2.fq.gz \
    > ./4.align/EXON_FYY.sam
# sam排序并轉(zhuǎn)為bam
samtools sort -@ 10 ./4.align/EXON_FYY.sam -o ./4.align/EXON_FYY.bam
ls -lh 4.align

Step 5: 比對結(jié)果進行質(zhì)控

可視化數(shù)據(jù)：快速、全面地對測序數(shù)據(jù)的比對質(zhì)量進行評估和質(zhì)控，確保輸入質(zhì)量是可靠的。
使用samtools stats“收集數(shù)據(jù)”，它生成了一個包含所有原始統(tǒng)計數(shù)據(jù)的文件；使用plot-bamstats自動生成一系列直觀的圖表

samtools stats --reference ./0.ref_38/hg38.fasta ./4.align/EXON_FYY.bam > ./5.stat/EXON_FYY.stat
plot-bamstats -p ./5.stat/EXON_FYY ./5.stat/EXON_FYY.stat

Step 6: 準備 gatk 文件

1. 刪除重復(fù)序列 (need large space & long times) (2h/file)

gatk MarkDuplicatesSpark -I ./4.align/EXON_FYY.bam -O ./6.gatk/EXON_FYY_dedup.bam --tmp-dir /home/chencx/test --remove-all-duplicates true
# 查看刪除前后的數(shù)目變化
samtools view ./4.align/EXON_FYY.bam |wc -l
samtools view ./6.gatk/EXON_FYY_dedup.bam |wc -l

2. 堿基質(zhì)量重矯正--消除測序儀器或者系統(tǒng)帶來的誤差

這部分可以提前下載好

# 2.1. 為人類全基因組構(gòu)建 samtools fasta 索引
samtools faidx ./0.ref_38/hg38.fasta
# 2.2. 為人類全基因組構(gòu)建 dictionary
picard CreateSequenceDictionary R=./0.ref_38/hg38.fasta \
O=./0.ref_38/hg38.dict
# 2.3. 給變異位點創(chuàng)建索引
wget -c https://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/00-All.vcf.gz -P ./0.ref_38
gzip -dk ./0.ref_38/00-All.vcf.gz
# hg38都帶有chr要添加上
awk 'BEGIN{OFS="\t"} /^#/{print; next} {$1="chr"$1; print}' ./0.ref_38/00-All.vcf > ./0.ref_38/00-All_chr.vcf
head -n 80 ./0.ref_38/00-All_chr.vcf

變異檢測之前的一個關(guān)鍵預(yù)處理步驟：校正測序數(shù)據(jù)中非隨機的、系統(tǒng)性的堿基質(zhì)量分數(shù)誤差，為下游的變異檢測提供更準確的數(shù)據(jù)，從而減少假陽性和假陰性。結(jié)果生成一個詳細的“重校準對照表”，指導(dǎo)下一個工具（ApplyBQSR）如何調(diào)整原始BAM文件中每一個堿基的質(zhì)量分數(shù)。

#2.4. 計算需重校正的reads和特征值 (need large space & long times)
gatk BaseRecalibrator -I ./6.gatk/EXON_FYY_dedup.bam \
-R ./0.ref_38/hg38.fasta \
-O ./6.gatk/table/EXON_FYY_recall.table \
--known-sites ./0.ref_38/00-All_chr.vcf \
--tmp-dir /home/chencx/test

實際修改BAM文件中的堿基質(zhì)量分數(shù)。它利用上一步（BaseRecalibrator）生成的“重校準配方”（.table文件），為每一個堿基計算并賦予一個新的、更準確的質(zhì)量分數(shù)，并輸出一個新的BAM文件。

# 2.5. 利用校準表文件重新調(diào)整堿基質(zhì)量，并使用這個新的質(zhì)量值重新輸出一份新的 BAM 文件
gatk ApplyBQSR -I ./6.gatk/EXON_FYY_dedup.bam \
-bqsr ./6.gatk/table/EXON_FYY_recall.table \
-O ./6.gatk/EXON_FYY_bqsr.bam \
-R ./0.ref_38/hg38.fasta 

至此，BQSR流程完成，產(chǎn)生的 EXON_FYY_bqsr.bam 文件就是進行Germline變異檢測（如用 HaplotypeCaller）或Somatic變異檢測（如用 Mutect2）的最佳輸入文件。

歡迎大家評論交流！
（每帖分享：時間，好像很緊張但又無極限）