WES，WGS等DNA测序数据找变异流程服务

2021-10-21 17:10:57 浏览数 (2)

肿瘤或者家系的WES，WGS等DNA测序样品的fastq数据，需要比对到参考基因组并且找变异并且注释，我们仅仅是收取一个计算机资源的费用，800-8000元人民币（根据样品数量不同收费不一样）即可，并且提供全套代码。不管是公共数据集还是你自己的实验测序数据，一样的费用！我们会代替你跑如下所示的流程：

date: 2021-05-14 13:37:00

whole exome-sequencing analysis pipeline

全外显子数据分析

(一) 环境搭建

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

代码语言：javascript复制
#查看一下环境中是否有java，如果有，版本是否符合要求
java --version

2、下载安装GATK4

代码语言：javascript复制
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、下载其他需要的软件

代码语言：javascript复制
#首先创建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文件

代码语言：javascript复制
#将下列数据下载的脚本保存到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官网提供的参考基因组

代码语言：javascript复制
#====step 1 首先质空，查看数据是否需要进行过滤
#我们这里的质控结果显示reads的整体质量还不错（GC含量看着不太过分基本不影响后续的分析），可以直接用于比对
ls ../raw/*.gz|xargs fastqc -t 4 -o ./
multiqc *.zip
代码语言：javascript复制
#====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 "@RGtID:$sampletSM:$sampletLB:WGStPL:Illumina" /$index $fq1 $fq2 |samtools sort -@ 8 -o ../align/$sample.sort.bam - &>>../align/$sample.log &
done
代码语言：javascript复制
# ====step 3 查看比对结果
ls *.bam|while read id;do samtools flagstat $id;done

2、标记PCR重复

代码语言：javascript复制
# ====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

代码语言：javascript复制
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。

代码语言：javascript复制
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文件

代码语言：javascript复制
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
代码语言：javascript复制
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
代码语言：javascript复制
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中一些可能由污染导致的假阳性突变

代码语言：javascript复制
nohup gatk CalculateContamination 
-I SRR12135912.pileups.table 
-matched SRR12135911.pileups.table 
-O HQ461-untreated.calculatecontamination.table &
step 3、 FilterMutectCalls

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

代码语言：javascript复制
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是绑定在一起的。

代码语言：javascript复制
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

代码语言：javascript复制
#运行下面的命令查看SnpEff目前支持的database
java -jar snpEff.jar databases | less
#运行下面的命令查看SnpEff支持的GRCh38 database,目前是GRCh38.99
java -jar snpEff.jar databases | grep -i GRCh38
#看来SnpEff网站更新不及时呀，在网站维护和更新这一点上，要重点表扬GATK团队！！！
代码语言：javascript复制
#于是下载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注释
代码语言：javascript复制
#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文件

代码语言：javascript复制
#保留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
代码语言：javascript复制
#简单看一下三个vcf文件的大小
ls -lh *.vcf
ls *.vcf|while read id; do cat $id|wc -l;done
step 3、sed命令复现SnpSift的功能

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

jar shell pipeline sed

0 人点赞