全外显子组测序分析 - 流程1

使用一个数据样本（EXON_FYY_1.fastq.gz，EXON_FYY_2.fastq.gz）作为示例，样本多可以写成循环。

Step 0：创建目录

原始测序数据放在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：质量控制

利用工具 fastqc、multiqc进行前期质量查看，确保测序数据基本合格。
html 文件是质控报告，zip 文件是 html 文件图表对应的类。样本多的话可以用 multiqc 对所有样本的结果整合。

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去接头

如果从质控报告看到一些 reads 的质量值较差或测到了接头，需要去除掉这些质量较差的 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

# 再质控：
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 表示经过质量修剪后的有效数据，即去除了低质量的部分。

Step 3: 构建参考基因组的 BWA 索引

这步骤时间较长可以提前运行。生成的 BWA 索引文件说明：
hg38.fasta.64.amb: 存储了参考基因组序列的一些元信息，如序列长度、名称等。
hg38.fasta.64.ann: 存储了一些额外的注释信息，有助于在比对过程中进行一些元数据的操作。
hg38.fasta.64.bwt: 这是Burrows-Wheeler变换的数据，是BWA索引的核心部分之一。BWT是一种数据转换算法，用于将原始DNA序列进行变换，以便在比对时能够更高效地进行匹配。
hg38.fasta.64.pac: 这是FM索引的数据，也是BWA索引的核心部分之一。FM索引是一种数据结构，用于实现快速查找和比对DNA序列。
hg38.fasta.64.sa: 这是后缀数组的数据，也是BWA索引的核心部分之一。后缀数组是一种用于存储DNA序列的子串信息的数据结构，用于在比对时快速查找匹配。

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 转换成二进制的 bam 文件，记录了每一条 reads 比对到参考基因组的结果。

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排序并转为bam
samtools sort -@ 10 ./4.align/EXON_FYY.sam -o ./4.align/EXON_FYY.bam
ls -lh 4.align

Step 5: 比对结果进行质控

可视化数据：快速、全面地对测序数据的比对质量进行评估和质控，确保输入质量是可靠的。
使用samtools stats“收集数据”，它生成了一个包含所有原始统计数据的文件；使用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. 删除重复序列 (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
# 查看删除前后的数目变化
samtools view ./4.align/EXON_FYY.bam |wc -l
samtools view ./6.gatk/EXON_FYY_dedup.bam |wc -l

2. 碱基质量重矫正--消除测序仪器或者系统带来的误差

这部分可以提前下载好

# 2.1. 为人类全基因组构建 samtools fasta 索引
samtools faidx ./0.ref_38/hg38.fasta
# 2.2. 为人类全基因组构建 dictionary
picard CreateSequenceDictionary R=./0.ref_38/hg38.fasta \
O=./0.ref_38/hg38.dict
# 2.3. 给变异位点创建索引
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

变异检测之前的一个关键预处理步骤：校正测序数据中非随机的、系统性的碱基质量分数误差，为下游的变异检测提供更准确的数据，从而减少假阳性和假阴性。结果生成一个详细的“重校准对照表”，指导下一个工具（ApplyBQSR）如何调整原始BAM文件中每一个碱基的质量分数。

#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文件中的碱基质量分数。它利用上一步（BaseRecalibrator）生成的“重校准配方”（.table文件），为每一个碱基计算并赋予一个新的、更准确的质量分数，并输出一个新的BAM文件。

# 2.5. 利用校准表文件重新调整碱基质量，并使用这个新的质量值重新输出一份新的 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流程完成，产生的 EXON_FYY_bqsr.bam 文件就是进行Germline变异检测（如用 HaplotypeCaller）或Somatic变异检测（如用 Mutect2）的最佳输入文件。

欢迎大家评论交流！
（每帖分享：时间，好像很紧张但又无极限）

全外显子组测序分析 - 流程1