肿瘤康复网 > 动物基因组测序基础分析流程总结（GWAS全流程第一部分：WGS基础流程）

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

时间：2022-06-05 19:15:09

相关推荐

一、前言：

研一进了纯生信课题组之后就开始了自学之路~，当前课题是做动物GWAS。因为自学的时候看了许多教程和代码都是针对人类基因组的，学习（copy）代码的时候各种大佬博主都是以人类基因组或者模式生物举例，如果做其他动物多多少少会有一些差别。

本文是本人学习过程中留下的随手笔记，不断更的话理论上包括三个部分：

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

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

从fastq到vcf过滤

二、从fastq到vcf

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

（一）软件安装

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

1.miniconda

下载地址
/miniconda/Miniconda3-latest-Linux-x86_64.sh

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

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

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

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

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

conda安装使用也可以参考这一篇Conda 安装使用详解，下载速度慢可以自行添加镜像。

2.BWA

/p/19f58a07e6f4?from=groupmessage #该网站有详细指南

3.GATK

conda install -c bioconda gatk

安装后输入gatk

Usage template for all tools (uses --spark-runner LOCAL when used with a Spark tool)gatk AnyTool toolArgsUsage template for Spark tools (will NOT work on non-Spark tools)gatk SparkTool toolArgs [ -- --spark-runner <LOCAL | SPARK | GCS> sparkArgs ]Getting helpgatk --list Print the list of available toolsgatk Tool --help Print help on a particular toolConfiguration File Specification--gatk-config-filePATH/TO/GATK/PROPERTIES/FILEgatk forwards commands to GATK and adds some sugar for submitting spark jobs--spark-runner <target> controls how spark tools are runvalid targets are:LOCAL:run using the in-memory spark runnerSPARK:run using spark-submit on an existing cluster --spark-master must be specified--spark-submit-command may be specified to control the Spark submit commandarguments to spark-submit may optionally be specified after -- GCS: run using Google cloud dataproccommands 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-runmay 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.

4.PLINK

conda install -c bioconda plink

输入plink

PLINK v1.90b6.24 64-bit (6 Jun ) www.cog-/plink/1.9/(C) - Shaun Purcell, Christopher Chang GNU General Public License v3plink <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).

出现如上文本即可。

（二）代码流程

放上的代码都省去了我自己的绝对路径和部分名称，如果要使用代码可以自己添加上自己的文件名和绝对路径。学习过程中也是感受到了读文档的重要性，很多查不到的问题在软件的文档或者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.txtrm -f 39.bamrm -f 39_sorted.bam#4 创建比对索引文件time samtools index 39_markdup.bam && echo "** index done **"

其实看到做人类基因组的教程里面还有一步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

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

因为我做的物种是羊，这一步结束了后我得到了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

没有服务器的话可以这样依次跑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

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

#5 对vcf进行mergels *genotype.vcf.gz > all_genotype.listgatk MergeVcfs -I all_genotype.list -O sheep.raw.vcf.gz

大功告成~

三、过滤

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

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

从测序位点质量上看，在GATK HaplotypeCaller之后，首选的质控方案是GATK VQSR，原理是利用自身数据和已知变异位点集的overlap，通过GMM模型构建一个分类器来对变异数据进行打分，从而评估每个位点的可行度。

虽然经常看到VQSR的教程，但是很可惜目前应该只有人类上可以做（还是需要高质量的已知变异集），所以我们其他物种只能做hard-filtering硬过滤了。

（一）硬过滤流程

以SNP的硬过滤流程为例

1. 首先我们选择出SNP

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

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 根据过滤结果选择SNPgatk SelectVariants 、-R ref.fa \-V sheep.snp.filter.vcf.gz \--exclude-filtered --restrict-alleles-to BIALLELIC \-O sheep.snp.pass.vcf.gz

到这里硬过滤就完成了，接下来是常规化的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

值得注意的是，对于除了人以外的其他动物来说，染色体条数千奇百怪，如果是常见动物比如我这里的羊，那么就需要用

"–sheep"来告诉软件这是羊的基因组，否则就会报错。“–allow-extra-chr"这条命令也是同样的功能，允许额外的染色体输入。”–recode vcf-iid "用以控制输出文件为vcf格式。

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