动物基因组测序基础分析流程总结（GWAS全流程第一部分：WGS基础流程）
一、前言：
研一进了纯生信课题组之后就开始了自学之路~，当前课题是做动物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 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-file PATH/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-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 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.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
动物基因组测序基础分析流程总结（GWAS全流程第一部分：WGS基础流程）相关推荐
- 学习全基因组测序数据分析1：测序技术
本文转载自微信公众号解螺旋的矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教科书级别的生物信息学习材料.虽然本平台只关注宏基因组领域,但此系列文章知识体系完善.干货满满,是值得每位专 ...
- 从零开始完整学习全基因组测序数据分析：第1节 测序技术
欢迎订阅我们的微信公众号:基因学苑 本文转载自微信公众号解螺旋矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教材级别的生物信息学习材料.我们将陆续转载给大家.大家也可以关注公众号解螺 ...
- 2020.10.21【转载】丨GWAS全基因组关联分析流程
感谢CSDN用户 追梦生信人 梳理了GWAS全基因组关联分析的整个流程,并提供了基本的命令,用到的软件包括BWA.samtools.gatk.Plink.Admixture.Tassel等,在此分享出 ...
- GWAS全基因组关联分析流程（BWA+samtools+gatk+Plink+Admixture+Tassel）
我梳理了GWAS全基因组关联分析的整个流程,并提供了基本的命令,用到的软件包括BWA.samtools.gatk.Plink.Admixture.Tassel等,在此分享出来给大家提供参考. 一.BW ...
- 全基因组重测序基础分析
本文仅为个人记录 流程说明 从零开始完整学习全基因组测序(WGS)数据分析:第3节 数据质控 2015a 2016a 全基因组测序数据获取后应该怎么分析 samtools samtools常用命令详解 ...
- 大熊猫源致病大肠杆菌CCHTP全基因组测序及耐药和毒力基因分析
大熊猫源致病大肠杆菌CCHTP全基因组测序及耐药和毒力基因分析 邓雯文1,李才武1,赵思越2,李仁贵1,何永果1,吴代福1,杨盛智2,黄炎1,张和民1,邹立扣2 1. 中国大熊猫保护研究中心,大熊猫国 ...
- 【3月30日直播】新冠病毒全基因组测序——Midnight试剂盒及整体解决方案
识别上方二维码 或点击「阅读原文」 免费报名参加 新冠疫情肆虐全球,基于Nanopore测序技术和数据分析在全球感染性疾病防控中的优势充分显现出来.该平台使用灵活.操作简便.产出快速.分析实时等特征为 ...
- 学习全基因组测序数据分析2：FASTA和FASTQ
本文转载自微信公众号解螺旋的矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教科书级别的生物信息学习材料.虽然本平台只关注宏基因组领域,但此系列文章知识体系完善.干货满满,是值得每位专 ...
- 《全基因组测序WGS数据分析——1.DNA测序技术》学习笔记
WGS(Whole Genome Sequencing) 指将物种细胞里面完整的基因组序列全部DNA,检测并排列,此技术几乎能够鉴定出基因组上任何类型的突变. 对于人类来说,全基因组测序的价值是极大的 ...
- 《全基因组测序WGS数据分析——4.构建WGS主流程》学习笔记
流程的具体形式其实是次要的,WGS本质上只是一个技术手段,重要的是,要明白自己所要解决的问题是什么,所希望获取的结果是什么,然后再选择合适的技术. 这是WGS数据分析的流程图.流程的目的是准确检测出每 ...
最新文章
- 常用的深度学习的linux代码(1.实时监测GPU情况2.当前正常使用的GPU情况3.杀掉特定某个进程4.杀掉特定某个进程)
- C++:类-多态的学习和使用
- SQL Server 2019安装教程
- vue组件调用(用npm安装)
- python反向代理服务器_主机、服务器，代理服务器，反向代理服务器理解（自用）...
- Nacos集群（二）阿里自研弱一致性Distro协议核心实现
- java se下载完怎么启动_【Java SE】如何安装JDK以及配置Java运行环境
- WPS office 下载
- MyBatis整合Spring的实现（16）
- 机票大讲堂之机票的秘密
- 9 椭圆曲线密码体制
- 大鱼吃小鱼c语言编程,scratch大鱼吃小鱼设计思路
- 【MySQL必知必会--理论】
- 【03】Linux笔记
- 大文件打包压缩成的几个小文件怎么解压？
- 《数据库系统概论》复习笔记
- UltraISO PE 绿色版9.1.2.2463
- Linux-selinux
- SpringMVC：通配符的匹配很全面, 但无法找到元素 ‘context:component-scan‘ 的声明
- 数据质量分析之校验规则模板