一、前言:
研一进了纯生信课题组之后就开始了自学之路~,当前课题是做动物GWAS。因为自学的时候看了许多教程和代码都是针对人类基因组的,学习(copy)代码的时候各种大佬博主都是以人类基因组或者模式生物举例,如果做其他动物多多少少会有一些差别。
本文是本人学习过程中留下的随手笔记,不断更的话理论上包括三个部分:
- 动物全基因组测序(WGS)基础流程
- GWAS流程和数据挖掘、可视化
- GWAS后续工作
这篇是关于动物全基因组测序(WGS)基础流程的笔记,简单来说就两件事
- 从fastq到vcf
- 过滤
二、从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
接下来创建并且激活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
因为要满足软件的各种奇奇怪怪的需求(版本支持),我会完全将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
安装后输入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.
4.PLINK
conda install -c bioconda plink
输入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).
出现如上文本即可。
(二)代码流程
放上的代码都省去了我自己的绝对路径和部分名称,如果要使用代码可以自己添加上自己的文件名和绝对路径。学习过程中也是感受到了读文档的重要性,很多查不到的问题在软件的文档或者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 **"
其实看到做人类基因组的教程里面还有一步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进行merge
ls *genotype.vcf.gz > all_genotype.list
gatk 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 选择出SNP
gatk 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 根据过滤结果选择SNP
gatk 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格式。
关于上文提到的批量运行的流程脚本都是,我有空会放上来~,有任何疑问和错误欢迎留言讨论喔。
四、参考资料
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