动物基因组测序基础分析流程总结（GWAS全流程第一部分：WGS基础流程）_gwas分析流程

一、前言：

研一进了纯生信课题组之后就开始了自学之路~，当前课题是做动物GWAS。因为自学的时候看了许多教程和代码都是针对人类基因组的，学习（copy）代码的时候各种大佬博主都是以人类基因组或者模式生物举例，如果做其他动物多多少少会有一些差别。
本文是本人学习过程中留下的随手笔记，不断更的话理论上包括三个部分：

1. 动物全基因组测序（WGS）基础流程
2. GWAS流程和数据挖掘、可视化
3. GWAS后续工作

这篇是关于动物全基因组测序（WGS）基础流程的笔记，简单来说就两件事

1. 从fastq到vcf
2. 过滤

二、从fastq到vcf

主线是三个部分：①原始数据质控 ②数据预处理 ③变异检测。我们从拿到质控好的下机数据（fastq）开始，对动物WGS基础流程进行梳理。

（一）软件安装

从fastq文件到vcf的过程中用到了以下几个软件，其中BWA好像只能到官网下载安装，附上了别人写的指导链接，GATK和PLINK是可以在conda库里找到的，可以进conda虚拟环境后直接用以下命令安装。

1.miniconda

下载地址
https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

下载之后将安装包移动到指定位置安装

#安装
bash Miniconda3-latest-Linux-x86_64.sh
#查看conda版本号
conda --version
1
2
3
4

接下来创建并且激活conda环境

#创建conda虚拟环境 --name或者省略为-n，后接命名即可。如果有python需求，可以创建时指定，如conda create --name your_env_name python=3.7
conda create --name your_env_name
#激活虚拟环境
conda activate your_env_name
#退出虚拟环境
conda deactivate
1
2
3
4
5
6

因为要满足软件的各种奇奇怪怪的需求（版本支持），我会完全将conda环境与系统环境隔离，所以通常没有将conda加入环境变量。我会先cd到miniconda/bin目录下，用source activate命令激活conda中的base环境，然后将这里作为日常使用的环境。此时，从我自定义的虚拟环境中用conda deactivate命令退出时，会默认退出到base环境中。
conda安装使用也可以参考这一篇Conda 安装使用详解，下载速度慢可以自行添加镜像。

2.BWA

https://www.jianshu.com/p/19f58a07e6f4?from=groupmessage #该网站有详细指南

3.GATK

conda install -c bioconda gatk
1

安装后输入gatk

 Usage template for all tools (uses --spark-runner LOCAL when used with a Spark tool)
    gatk AnyTool toolArgs

 Usage template for Spark tools (will NOT work on non-Spark tools)
    gatk SparkTool toolArgs  [ -- --spark-runner <LOCAL | SPARK | GCS> sparkArgs ]

 Getting help
    gatk --list       Print the list of available tools

    gatk Tool --help  Print help on a particular tool

 Configuration File Specification
     --gatk-config-file                PATH/TO/GATK/PROPERTIES/FILE

 gatk forwards commands to GATK and adds some sugar for submitting spark jobs

   --spark-runner <target>    controls how spark tools are run
     valid targets are:
     LOCAL:      run using the in-memory spark runner
     SPARK:      run using spark-submit on an existing cluster 
                 --spark-master must be specified
                 --spark-submit-command may be specified to control the Spark submit command
                 arguments to spark-submit may optionally be specified after -- 
     GCS:        run using Google cloud dataproc
                 commands after the -- will be passed to dataproc
                 --cluster <your-cluster> must be specified after the --
                 spark properties and some common spark-submit parameters will be translated 
                 to dataproc equivalents

   --dry-run      may be specified to output the generated command line without running it
   --java-options 'OPTION1[ OPTION2=Y ... ]'   optional - pass the given string of options to the 
                 java JVM at runtime.  
                 Java options MUST be passed inside a single string with space-separated values.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

4.PLINK

conda install -c bioconda plink
1

输入plink

PLINK v1.90b6.24 64-bit (6 Jun 2021)           www.cog-genomics.org/plink/1.9/
(C) 2005-2021 Shaun Purcell, Christopher Chang   GNU General Public License v3

  plink <input flag(s)...> [command flag(s)...] [other flag(s)...]
  plink --help [flag name(s)...]

Commands include --make-bed, --recode, --flip-scan, --merge-list,
--write-snplist, --list-duplicate-vars, --freqx, --missing, --test-mishap,
--hardy, --mendel, --ibc, --impute-sex, --indep-pairphase, --r2, --show-tags,
--blocks, --distance, --genome, --homozyg, --make-rel, --make-grm-gz,
--rel-cutoff, --cluster, --pca, --neighbour, --ibs-test, --regress-distance,
--model, --bd, --gxe, --logistic, --dosage, --lasso, --test-missing,
--make-perm-pheno, --tdt, --qfam, --annotate, --clump, --gene-report,
--meta-analysis, --epistasis, --fast-epistasis, and --score.

"plink --help | more" describes all functions (warning: long).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

出现如上文本即可。

（二）代码流程

放上的代码都省去了我自己的绝对路径和部分名称，如果要使用代码可以自己添加上自己的文件名和绝对路径。学习过程中也是感受到了读文档的重要性，很多查不到的问题在软件的文档或者github上都写得清清楚楚了~。

1. 为Reference 文件建立索引

2. 排序

3. 标记重复

4. 创建比对索引文件

#1 比对 加time是为了看时间，也可以不加（下同）
time bwa mem -t 4 -R '@RG\tID:39_1\tPL:illumina\tSM:39' 39_ref.fa 39_1.fq.gz 39_2.fq.gz | samtools view -Sb - > 39.bam && echo "** bwa mapping done **"
#2 ①用gatk排序
time gatk SortSam -I 39.bam -O 39.sort.bam --SORT_ORDER coordinate --TMP_DIR sort/39
#2 ②也可以用samtools，选一个就行。12.29日发现有同学samtools出来的文件报错，进行了修改
time samtools sort -@ 4 -m 4G -O BAM -o 39_sorted.bam 39.bam && echo "** BAM sort done"rm -f 39.bam 
#3 标记PCR重复
time gatk MarkDuplicates -I 39_sorted.bam -O 39_markdup.bam -M 39_markdup_metrics.txt
rm -f 39.bam
rm -f 39_sorted.bam
#4 创建比对索引文件
time samtools index 39_markdup.bam && echo "** index done **"
1
2
3
4
5
6
7
8
9
10
11
12

其实看到做人类基因组的教程里面还有一步BQSR(Recalibration Base Quality Score)，利用已有的snp数据库，建立相关性模型。可惜我做的物种并没有这么强大的snp数据库~，就跳过啦。接下来直接生成vcf文件

5. 单个体生成vcf文件

如果是单个体，那么到这里就可以直接生成vcf文件了，命令如下：

time gatk HaplotypeCaller -R 39_ref.fa -I 39_markdup.bam -O 39.vcf
1

6. 生成gvcf文件并合并

有多个个体的情况下，首先我们用HaplotypeCaller命令给每一个个体生成gvcf文件，然后用CombineGVCFs按染色体合并gvcf文件。这里我拿第39号个体和第20号染色体做示范，代码如下。

我是用python生成的批量shell脚本然后放到服务器上一起跑的，就会比较快。批量写脚本什么语言都可以，就按照你的需求循环改变名称、个体号、染色体号等各种字符就行。

#1 生成gvcf文件
gatk  HaplotypeCaller -R ref.fa -I 39.rmdup.bam --emit-ref-confidence GVCF --min-base-quality-score 10 -O 39.g.vcf.gz
#2 gvcf文件按染色体合并
ls *.g.vcf.gz > all_gvcf 
gatk CombineGVCFs -R ref.fa -V all_gvcf.list  -L 20 -O chr20.merged.g.vcf.gz

1
2
3
4
5
6

因为我做的物种是羊，这一步结束了后我得到了26条以chr$i.merged.vcf.gz命名的分染色体合并完成的gvcf文件。接下来召唤基因型~

#3 生成基因型文件（这一步变成vcf了）
gatk GenotypeGVCFs -R ref.fa -V chr20.merged.g.vcf.gz -O chr20.genotype.vcf.gz
1
2

没有服务器的话可以这样依次跑26条染色体

for i in {1..26}; do gatk GenotypeGVCFs -R ref.fa -V chr$i.merged.g.vcf.gz -O chr$i.genotype.vcf.gz;done
1

有服务器的话还是建议批量写26个脚本一起提交，这一步大概一两天能完成（前几条染色体比较大，需要的时间是后面染色体的几倍）。

#5 对vcf进行merge
ls *genotype.vcf.gz > all_genotype.list
gatk MergeVcfs -I all_genotype.list -O sheep.raw.vcf.gz
1
2
3

大功告成~

三、过滤

上一步得到的是rawdata，我们还需要对原始数据做变异质控。质控是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。

SNP过滤有两种情况，一种是仅根据位点质量信息（测序深度，回帖质量等）对SNP进行粗过滤。另一种是考虑除了测序质量以外的信息进行的过滤。 接下来我会分别介绍到这两种过滤。

从测序位点质量上看，在GATK HaplotypeCaller之后，首选的质控方案是GATK VQSR， 原理是利用自身数据和已知变异位点集的overlap，通过GMM模型构建一个分类器来对变异数据进行打分，从而评估每个位点的可行度。
虽然经常看到VQSR的教程，但是很可惜目前应该只有人类上可以做（还是需要高质量的已知变异集），所以我们其他物种只能做hard-filtering硬过滤了。

（一）硬过滤流程

以SNP的硬过滤流程为例

1. 首先我们选择出SNP

#1 选择出SNP
gatk SelectVariants -select-type SNP -V sheep.raw.vcf.gz -O sheep.raw.snp.vcf.gz
1
2

2. 硬过滤条件

过滤条件可以在GATK的官网找到，这个链接是官方文档对于过滤条件选择的解释，每一条都有。GATK硬过滤条件
在评论区找到了一个总结，截图如下表。

根据这个条件我们进行硬过滤

#2 执行硬过滤
gatk VariantFiltration \
-R ref.fa \
-V sheep.raw.snp.vcf.gz \
--filter-expression "QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name LowQua \
-O sheep.snp.filter.vcf.gz
#3 根据过滤结果选择SNP
gatk SelectVariants 、
-R ref.fa \
-V sheep.snp.filter.vcf.gz \
--exclude-filtered --restrict-alleles-to BIALLELIC \
-O  sheep.snp.pass.vcf.gz
1
2
3
4
5
6
7
8
9
10
11
12
13

到这里硬过滤就完成了，接下来是常规化的MAF、缺失率等SNP过滤。

（二）根据MAF、缺失率进行过滤

这部分有很多软件都可以完成，我用plink做个示范，指定MAF（次要等位基因）过滤阈值为0.05，基因型缺失的过滤阈值为0.1。

#1 MAF、缺失率过滤
plink \
--vcf sheep.snp.pass.vcf.gz \
--sheep \  
--keep-allele-order \
--allow-extra-chr \
--geno 0.1 --maf 0.05\
--recode vcf-iid \
--out sheep.mis0.1.maf0.05.vcf
1
2
3
4
5
6
7
8
9

值得注意的是，对于除了人以外的其他动物来说，染色体条数千奇百怪，如果是常见动物比如我这里的羊，那么就需要用
"–sheep"来告诉软件这是羊的基因组，否则就会报错。“–allow-extra-chr"这条命令也是同样的功能，允许额外的染色体输入。”–recode vcf-iid "用以控制输出文件为vcf格式。

关于上文提到的批量运行的流程脚本都是，我有空会放上来~，有任何疑问和错误欢迎留言讨论喔。

四、参考资料

GATK 人类基因组标准流程 https://zhuanlan.zhihu.com/p/69726572?from_voters_page=true
GWAS（from 橙子牛奶糖）https://www.cnblogs.com/chenwenyan/p/11803311.html
全基因组测序(WGS)流程及实践 https://zhuanlan.zhihu.com/p/335770966