WES（二）——数据预处理

0 创建目录

# directories目录
ref="/data/supporting_files/hg38/hg38.fa"
known_sites="/data/supporting_files/hg38/Homo_sapiens_assembly38.dbsnp138.vcf"
aligned_reads="/data/aligned_reads"
reads="/data/reads"
results="/data/results"
data="/data/WES/data"

1. 原始数据质控

可以用trimmomatic或FastQC软件进行质控。推荐使用fastp，其使用简单、运速较快，默认可以不设置任何参数，就能够同时进行序列过滤、校正并产生质控报告。

# ----------------- STEP 1: QC - Run fastqcs 数据质控 -------------------------#java -jar /software/Trimmomatic-0.39/trimmomatic-0.39.jar  PE -threads 36 -phred33  ${i}_R1.fastq.gz  ${i}_R2.fastq.gz /data/${i}_1.clean.fq.gz  /data/${i}_unpaired1.fq /data/${i}_2.clean.fq.gz  /data/${i}_unpaired2.fq ILLUMINACLIP:/software/Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36  HEADCROP:10fastqc ${reads}/${i}_1.filt.fastq.gz -o ${reads}/  #正向测序文件
fastqc ${reads}/${i}_2.filt.fastq.gz -o ${reads}/   #反向测序文件
#生成文件为：.html和.zip文件。先检查看一下html文件看是否需要进一步QC

2. 使用BWA进行比对

# ------------- STEP 2: Map to reference using BWA-MEM 唯一比对-----------------# BWA index reference索引参考文件 
bwa index ${ref}# BWA alignment比对
bwa mem -M -t 10 -R "@RG\tID:SRR062634\tPL:ILLUMINA\tS：M:SRR062634" ${ref} ${reads}/${i}_1.filt.fastq.gz ${reads}/${i}_2.filt.fastq.gz > ${aligned_reads}/${i}.sam
samtools view -bS 02.sam/${i}.sam | samtools sort -o 03.bam/${i}.bam#samtools view HCC16.sam | less     #查看文件
#samtools flagstat ${i}.paired.sam   #查看文件总结信息

3. 标记冗余序列并排序

gatk的MarkDuplicates模块只把重复序列在输出的新结果中标记出来，但不删除。

把参数--remove-all-duplicates设置为ture，那么重复序列就被删除掉，不会在结果文件中留存。建议只标记出来不删除这些冗余序列，以便在需要的时候还可以对其做分析。

# ---------- STEP 3: Mark Duplicates and Sort - GATK4 标记去除冗余序列 ------------gatk MarkDuplicatesSpark \  # MarkDuplicatesSpark函数对重复reads标记并对sam文件进行标记
-I ${aligned_reads}/${i}.bam \
-O ${aligned_reads}/${i}_sorted_dedup_reads.bam \
-M  ${aligned_reads}/${i}_sorted_dedup_reads.metrics \
#--remove-all-duplicates false

4. 创建比对索引文件

为sorted_dedup_reads.bam创建索引文件，方便随机访问文件中的任意位置，而且后面的“局部重比对”步骤也要求这个BAM文件一定要有索引。

完成之后，会生成一份sorted_dedup_reads.bam.bai文件。

samtools index ${i}_sorted_dedup_reads.bam

5. Indel局部区域重新比对

在进行这一步骤之前还有一个merge的操作，将同个样本的所有比对结果合并成唯一一个大的BAM文件。之所以会有这种情况，是因为有些样本测得非常深，其测序结果需要经过多次测序（或者分布在多个不同的测序lane中）才全部获得，这个时候我们一般会先分别进行比对并去除重复序列后再使用samtools进行合并。

samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

【Indel局部区域重新比对】非必须，我们后面的变异检测使用GATK的HaplotypeCaller模块，仅当这个时候才可以减少这个Indel局部重比对的步骤。原因是GATK的HaplotypeCaller中，会对潜在的变异区域进行相同的局部重比对！但是其它的变异检测工具或者GATK的其它模块就没有这么干了！所以切记！

6. 碱基质量重校正

#GATK提供了BQSR，这一步以用于矫正测序碱基的“真实质量值”。BQSR会设置参考SNP数据集作为“真实”变异集，这些数据集之外的位点，若测序数据与参考基因组存在差异，则认为是“错误”，以此来判定并矫正测序碱基的质量值。

第一步，BaseRecalibrator建模，计算出所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件(sample_name.recal_data.table)。

第二步，ApplyBQSR运行，这一步利用第一步得到的校准表文件(sample_name.recal_data.table)重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新BAM文件。

# --------------- STEP 4: Base quality Score recalibration 碱基质量校正 ----------------# 1. build the model 建模
gatk BaseRecalibrator \
-I ${aligned_reads}/SRR062634_sorted_dedup_reads.bam \
-R ${ref} \
--known-sites dbsnp_146.hg38.vcf.gz \
--known-sites 1000G_phase1.snps.high_confidence.hg38.vcf.gz \
--known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \  #已知的参考文件
-O ${data}/recal_data.table# 2. Apply the model to adjust the base quality scores应用模型调整碱基质量评分
gatk ApplyBQSR \
-I ${aligned_reads}/SRR062634_sorted_dedup_reads.bam 
-R ${ref} \
--bqsr-recal-file {$data}/recal_data.table \
-O ${aligned_reads}/SRR062634_sorted_dedup_bqsr_reads.bam 

参考

全基因组测序(WGS)流程及实践 - 知乎

全基因组测序简要分析流程 - 沈文龙的博客 | Wenlong Shen's Blog

从零开始完整学习全基因组测序（WGS）数据分析：第4节 构建WGS主流程 | Public Library of Bioinformatics