WES

whole exome-sequencing analysis pipeline

全外显子数据分析

(一) 环境搭建

1、GATK依赖Java 8/JDK 1.8 (Oracle or OpenJDK)

#查看一下环境中是否有java，如果有，版本是否符合要求
java --version

2、下载安装GATK4

wget https://github.com/broadinstitute/gatk/releases/download/4.2.0.0/gatk-4.2.0.0.zip
unzip gatk-4.2.0.0.zip
#将当前路径加到环境变量中
echo 'export PATH=/home/gongyuqi/biosoft/GATK/gatk-4.2.0.0:$PATH' >> ~/.bashrc
source ~/.bashrc

3、下载其他需要的软件

#首先创建WES的conda环境
conda create -n WES
#其次下载其他需要的软件
conda install -y python=3.6.2
conda install -y bwa sra-tools samtools bcftools snpEFF multiqc 
qualimap  
#激活环境，准备开始实战演练
conda activate WES
#创建存放各阶段数据的文件夹
cd /home/gongyuqi/project/WES
mkdir raw qc clean mutaion
cd qc && mkdir raw_qc clean_qc

(二) WES测试数据下载

1、数据来源GSE153707，我们这里从EBI直接下载fastq文件

#将下列数据下载的脚本保存到download.sh文件中
dir=/home/gongyuqi/.aspera/connect/etc/asperaweb_id_dsa.openssh
x=_1
y=_2
for id in {11,12}
do
ascp -QT -l 300m -P33001 -i $dir era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR121/0$id/SRR121359$id/SRR121359$id$x.fastq.gz .
ascp -QT -l 300m -P33001 -i $dir era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR121/0$id/SRR121359$id/SRR121359$id$y.fastq.gz . 
done
#赋予执行权限
chmod +x download.sh
#放后台进行数据下载
nohup ./download.sh &

2、了解文章WES数据处理的相关步骤和参数

(三) GATK准备bam文件用于找变异

1、 比对GATK官网提供的参考基因组

#====step 1 首先质空，查看数据是否需要进行过滤
#我们这里的质控结果显示reads的整体质量还不错（GC含量看着不太过分基本不影响后续的分析），可以直接用于比对
ls ../raw/*.gz|xargs fastqc -t 4 -o ./
multiqc *.zip

#====step 2 比对
ls *1.fastq.gz >fq1.txt
ls *2.fastq.gz >fq2.txt
paste fq1.txt fq2.txt > sample.txt
index=/home/gongyuqi/ref/GATK/hg38/v0/Homo_sapiens_assembly38.fasta.64
cat sample.txt | while read id
do nohup
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
sample=${fq1%%_*}
bwa mem -t 8 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" /$index $fq1 $fq2 |samtools sort -@ 8 -o ../align/$sample.sort.bam - &>>../align/$sample.log &
done

# ====step 3 查看比对结果
ls *.bam|while read id;do samtools flagstat $id;done

2、标记PCR重复

# ====step 1 MarkDuplicates
ls *.bam|while read sample
do nohup gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" \
   MarkDuplicates \
   --INPUT $sample \
   --OUTPUT ${sample%%.*}_marked.bam \
   -M ${sample%%.*}.metrics &
done

# ====step 2 FixMateInformation
ls *_marked.bam|while read sample
do nohup gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" \
   FixMateInformation \
   -I ${sample} \
   -O ${sample%%_*}_marked_fixed.bam \
   -SO coordinate &
done
#建立*_fixed.bam文件索引
ls *_fixed.bam|while read sample
do nohup samtools index $sample &
done

3、碱基矫正 BaseRecalibrator

ref=/home/gongyuqi/ref/GATK/hg38/v0/Homo_sapiens_assembly38.fasta
snp=/home/gongyuqi/ref/GATK/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf
indel=/home/gongyuqi/ref/GATK/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz

ls *_marked_fixed.bam|while read sample
do nohup gatk \
         --java-options "-Xmx20G -Djava.io.tmpdir=./" \
         BaseRecalibrator \
         -R $ref \
         -I $sample \
         --known-sites $snp \
         --known-sites $indel \
         -O ${sample%%_*}_recal.table \
         1>${sample%%_*}_log.recal 2>&1 &
done

#BWA比对时设置@RG的重要性在这一步体现出来了，算法通过@RG中的ID来识别各个独立的测序过程
ls *_marked_fixed.bam|while read sample
do nohup gatk \
         --java-options "-Xmx20G -Djava.io.tmpdir=./"   ApplyBQSR \
         -R $ref  \
         -I $sample  \
         -bqsr ${sample%%_*}_recal.table \
         -O ${sample%%_*}_bqsr.bam \
         1>${sample%%_*}_log.ApplyBQSR  2>&1 &
done
# 查看比对结果
ls *bqsr.bam|while read id
do nohup samtools flagstat $id &
done

gatk --java-options "-Xmx20G" AnalyzeCovariates \
     -before recal_data.table1 \
     -after recal_data.table2 \
     -plots AnalyzeCovariates.pdf

4、HaplotypeCaller

这一步是使用GATK的HaplotypeCaller找变异，但现在GATK官网推荐GATK4+Mutect2找变异，在这里我们还是简单运行一下HaplotypeCaller。

ls *bqsr.bam|while read sample
do nohup gatk \
         --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \
         -R $ref  \
         -I $sample \
         --dbsnp $snp \
         -O ../vcf/${sample%%_*}_raw.vcf \
         1>../vcf/${sample%%_*}_log.HC 2>&1 &
done

(三) Mutect2找变异

Mutect2 v4.1.0.0的教程目前已经不被GATK官方强烈推荐。官方更新推荐了Mutect2 v4.1.1.0 and later版本的使用教程，更新过的教程比较推荐用于正常对照样本>40的WES测序数据。如果没有这么多的样本量，Mutect2 v4.1.0.0的教程也同样可以使用，只是有部分参数在新的Mutect2版本中被弃用或是将其功能整合到另外的参数去了。具体问题可以查看官网说明。

1、过滤种系突变

首先需要一个germline variant sites VCF文件，去官网下载af-only-gnomad.hg38.vcf.gz文件。记得一并下载对应的.tbi文件

ref=/home/gongyuqi/ref/GATK/hg38/v0/Homo_sapiens_assembly38.fasta
germine_vcf=/home/gongyuqi/ref/GATK/hg38/af-only-gnomad/af-only-gnomad.hg38.vcf.gz

nohup gatk --java-options "-Xmx20g" Mutect2 \
-R $ref \
-I SRR12135912_bqsr.bam \
-tumor SRR12135912 \
-I SRR12135911_bqsr.bam \
-normal SRR12135911 \
--germline-resource $germine_vcf \
--af-of-alleles-not-in-resource 0.0000025 \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
--bam-output HQ461-untreated-Mutect2.bam \
-O ../Mutect2/HQ461-untreated.Mutect2.vcf &

生成文件如下：

2、考虑样品纯度及假阳性

通常，病人的肿瘤样本是混合有正常组织的，这个时候我们可以选择GetPileupSummaries工具来计算肿瘤样本的污染程度。对于Mutect2流程而言，GATK只考虑了肿瘤样本中正常样本的污染情况，不考虑正常样本中肿瘤样本的污染情况。这里我们的处理组:HQ461,对照组untreated。，虽然我们的样本是肿瘤细胞的单克隆株，我们这里还是跑一下这个流程。

step 1、生成af-only-gnomad.hg38.SNP_biallelic.vcf.gz

ref=/home/gongyuqi/ref/GATK/hg38/v0/Homo_sapiens_assembly38.fasta
germine_vcf=/home/gongyuqi/ref/GATK/hg38/af-only-gnomad/af-only-gnomad.hg38.vcf.gz
nohup gatk SelectVariants \
-R $ref \
-V $germine_vcf \
--select-type-to-include SNP \
--restrict-alleles-to BIALLELIC \
-O /home/gongyuqi/ref/GATK/hg38/af-only-gnomad/af-only-gnomad.hg38.SNP_biallelic.vcf.gz &

step 2、GetPileupSummaries

germine_biallelic_vcf=/home/gongyuqi/ref/GATK/hg38/af-only-gnomad/af-only-gnomad.hg38.SNP_biallelic.vcf.gz
genomic_intervals=/home/gongyuqi/ref/GATK/hg38/v0/wgs_calling_regions.hg38.interval_list

nohup gatk GetPileupSummaries \
-I SRR12135911_bqsr.bam \
-L $genomic_intervals \
-V $germine_biallelic_vcf \
-O ../Mutect2/SRR12135911.pileups.table &

nohup gatk GetPileupSummaries \
-I SRR12135912_bqsr.bam \
-L $genomic_intervals \
-V $germine_biallelic_vcf \
-O ../Mutect2/SRR12135912.pileups.table &

step 3、 CalculateContamination

计算污染比例，用以过滤掉somatic variant中一些可能由污染导致的假阳性突变

nohup gatk CalculateContamination \
-I SRR12135912.pileups.table \
-matched SRR12135911.pileups.table \
-O HQ461-untreated.calculatecontamination.table &

step 3、 FilterMutectCalls

结合CalculateContamination的评估，进行最后的过滤

ref=/home/gongyuqi/ref/GATK/hg38/v0/Homo_sapiens_assembly38.fasta
nohup gatk FilterMutectCalls \
-R $ref \
-V HQ461-untreated.Mutect2.vcf \
--contamination-table HQ461-untreated.calculatecontamination.table \
-O HQ461-untreated.filtered.Mutect2.vcf &

至此，我们找变异的步骤就差不多走完了。你会发现我们并没有使用到CollectSequencingArtifactMetrics和FilterByOrientationBias参数，那是因为我们当前版本的Mutect2已经将FilterByOrientationBias的功能整合进FilterMutectCalls中了，详情见GATK论坛中的一篇帖子。

另外，在做WES分析时，对照组是很重要的。拿肿瘤样本举例，tumor-only的模式会得到很多很多很多的假阳性，详情见GATK论坛中的另一篇帖子。对照组的存在会大大降低假阳性率。

(四) SnpEff & SnpSift注释

以下流程参考SnpEff & SnpSift官网文档

1、SnpEff软件及所需注释数据库的下载

step 1、下载安装SnpEff软件

下载SnpEff，现在SnpEff和SnpSift是绑定在一起的。

cd ~/biosoft/
#下载
wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip &
#安装
unzip snpEff_latest_core.zip

step 2、下载SnpEff软件需要的数据库文件

下载 SnpEff databases: 官网给的命令是java -jar snpEff.jar download GRCh38.76。但是运行这个命令会报错：找不到GRCh38.76

#运行下面的命令查看SnpEff目前支持的database
java -jar snpEff.jar databases | less
#运行下面的命令查看SnpEff支持的GRCh38 database,目前是GRCh38.99
java -jar snpEff.jar databases | grep -i GRCh38
#看来SnpEff网站更新不及时呀，在网站维护和更新这一点上，要重点表扬GATK团队！！！

#于是下载GRCh38.99
nohup java -jar snpEff.jar download GRCh38.99 &
#但是这个网络问题是无解的。我最后还是自己手动下载snpEff_v5_0_GRCh38.99.zip文件
#解压完，最后GRCh38.99文件夹绝对路径/home/gongyuqi/biosoft/snpEff/data/GRCh38.99

GRCh38.99内容如下

2、SnpEff、SnpSift使用

step 1、SnpEff注释

#the "verbose" mode (command line option -v), this makes SnpEff to show a lot of information which can be useful for debugging.
java -Xmx16g -jar snpEff.jar -v GRCh38.99 /home/gongyuqi/project/WES/Mutect2/HQ461-untreated.filtered.Mutect2.vcf > HQ461-untreated.filtered.ann.vcf

step 2、SnpSift过滤

该部分同时参考官网说明文档及某优秀博客用SnpSift过滤VCF文件

#保留Filter字段为'PASS'或缺失值的记录
cat HQ461-untreated.filtered.ann.vcf | java -jar SnpSift.jar filter "( na FILTER ) | (FILTER = 'PASS')" > HQ461-untreated.filtered.ann.snpsift.vcf
#保留INFO中ANN字段的IMPACT为'HIGH'或'MODERATE'的记录
cat HQ461-untreated.filtered.ann.snpsift.vcf | java -jar SnpSift.jar filter "(ANN[0].IMPACT has 'HIGH') | (ANN[0].IMPACT has 'MODERATE')" > HQ461-untreated.filtered.ann.snpsift.second.vcf

#简单看一下三个vcf文件的大小
ls -lh *.vcf
ls *.vcf|while read id; do cat $id|wc -l;done

step 3、sed命令复现SnpSift的功能

这里不难发现，SnpSift其实就是对文本文件的处理。要是shell脚本够扎实，也完全可以不依赖SnpSift。重要的是，在学习NGS组学分析过程中，但凡有锻炼自己shell脚本的地方，就一定要抓住机会动手写一写💪 这里我们可以用一个简单的sed命令解决

#第一个`sed`命令打印出所有标记为PASS的位点
#第二个`sed`命令承接上一个管道输出的内容，保留标记为HIGH或者MODERATE的位点
cat HQ461-untreated.filtered.ann.vcf | grep -v "#"| \
sed -n '/PASS/p' - | \
sed -n '/\(HIGH\|MODERATE\)/p' - \
>test_PASS_HIGH_MODERATE.third.vcf

#检验一下我们`sed`脚本的正确性
#看看位点数目是否一致
cat HQ461-untreated.filtered.ann.snpsift.second.vcf|grep -v "#"|wc -l
cat test_PASS_HIGH_MODERATE.second.vcf|wc -l

#再具体看看内容是否一致
cat HQ461-untreated.filtered.ann.snpsift.second.vcf|grep -v "#"|less -SN
cat test_PASS_HIGH_MODERATE.second.vcf|less -SN

HQ461-untreated.filtered.ann.snpsift.second.vcf

test_PASS_HIGH_MODERATE.second.vcf

(五) 原文数据复现

真正接受检验的时候到了，看看我们的pipeline能不能复现文章的数据。这也是进一步检验我们pipeline的正确性的重要环节（当然也不能迷信文章的结果，但是7+文章的数据分析质量应该还是可以作为一个较好的参考）。

再看看我们过滤后的vcf文件中是否有检测到CDK12的G到A的突变

学习WES心得体会

1、学习某个软件的用法时，主要参考官方文档！！！可以选择性参考相关博客主博客，但切记拿来主义！！！深刻理解用到的每个参数。

2、WES的相关分析还有很多，我们这里只涉及到SNV。同样的，突变的注释也涉及到多种软件，我们这里只涉及到了SnpEff。学海无涯啊~~~从培养起自主学习的思维和习惯开始。

3、拿到最终的注释文件（.vcf），我们还可以用到maftools、ComplexHeatmap等R包将突变情况进行可视化。

4、不急不慢，不慌不忙，坚持下去❤