biostar handbook(十)|如何进行变异检测

2024-06-23 21:08

本文主要是介绍biostar handbook(十)|如何进行变异检测,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

2013053-e65d78595eb8ffff.png
变异检测流程

什么是基因组变异

基因组变异是一个定义比较模糊的概念. 所谓的变异是相对于一个完美的“参考基因组”而言。但是其实完美的“参考基因组”并不存在,因为我们只是选择某一个物种里的其中似乎比较正常的个体进行测序组装,然后基于它进行后续的研究。简单的说,参考和变异是相对而言,变异也可能完全正常。

常见的基因组变异一般可以归为如下几类:

  • SNP, 单核苷酸多态性, 一个碱基的变化
  • INDEL,插入或缺失, 一个碱基的增加或移除
  • SNV, 单核苷酸变异,一个碱基的改变,可以是SNP,也可以是INDEL
  • MNP, 多核苷酸多态性,一个区块中有多个保守的SNP
  • MNV,多核苷酸变异,一个区块中有多个SNP或INDEL
  • short variations, 小于50bp的变异
  • large-scale variation, 大于50bp的大规模变异
  • SV, 结构变异,通常是上千个碱基,甚至是染色体级别上的变异

研究这些变异需要用到不同的手段,其中普通的DNA二代测序在寻找20bp以下的变异比较靠谱,对于大于1kb的结构变异而言,采用光学图谱(Bionanogenomics)可能更加靠谱一点。因此,对于目前最常用的二代测序而言,还是尽量就找SNP和INDEL吧,几个碱基的变化找起来还是相对容易些和靠谱些。

对于SNP的定义这里也要注意下,最初的SNP的定义指的是单个碱基导致的多态性,在群体中广泛存在(1%),可用来作为分子标记来区分不同个体。目前的定义比较粗暴一点,就是那个和“参考基因组”不同的单个位点。值得注意的是,这个概念可能不同人还有不同的定义,当你和别人就某个问题争执的时候,最好问问他是如何定义这个基本概念。由于SNP的广泛存在,并且变异可能会导致疾病,也就是存在某些SNP会导致疾病。Online Mendelian Inheritance in Man,就是一个人类遗传疾病数据库,建议去看下。

最后说下genotype和haplotype。genotype,基因型指的是一个个体的遗传组成。但是对于基因组变异而言,基因型通常指的是个体在某个位点上的等位基因情况。haplotype, 单倍型最初指的是从单个亲本中遗传的一组基因,而在基因组变异背景下,则是指一组变异。

一次简单的变异检测实战

变异检测(variant calling)即通过比较参考序列和比对结果来找到两者的不同并记录,基本上可以分为如下几步:

  • 序列比对
  • 比对后处理(可选)
  • 从联配中确定变异
  • 根据某些标准进行过滤
  • 对过滤的变异注释

这里面的每一个可选的工具都有很多,不同工具组合后的分析流程得到的结果可能会有很大差异。在变异检测这一部分目前就有很多软件,但是常用并且相对比较可靠的工具有如下几个:

  • bcftools: http://www.htslib.org/doc/bcftools.html
  • FreeBayes: https://github.com/ekg/freebayes
  • GATK: https://software.broadinstitute.org/gatk/
  • VarScan2: http://varscan.sourceforge.net/

当然这些工具最初都是用于人类基因组。

以埃博拉基因组为例完成一次简单的Variant Calling,所需工具为efetch, fastq-dump, emboss/seqretbwa, samtoolsFreeBayessnpEff。这些都可以通过conda快速安装。

第一步: 获取参考基因组序列,并建立索引

# 建立文件夹
mkdir -p refs
# 根据Accession下载
ACC=AF086833
REF=refs/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC | seqret -filter -sid $ACC > $REF
bwa index $REF

第二步: 获取需要比对的测序数据, 以前10w条为例

# 仅要前10w条read
SRR=SRR1553500
fastq-dump -X 100000 --split-files $SRR

第三步:序列比对

BAM=$SRR.bam
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM

第四步:使用freebayes或HaplotypeCaller(GATK4)检测变异

freebayes -f $REF $BAM > ${SRR}_freebayes.vcf
gatk HaplotypeCaller -R $REF-I $BAM -stand-call-conf 30 \-bamout bamout.bam--genotyping-mode DISCOVERY-O ${SRR}_haplotypecaller.vcf

用IGV可视化的效果如下:

2013053-e372b9241c285a31.jpg
IGV对BAM和VCF文件进行可视化

这是最简单的变异检测流程,对于找到的变异还可以进一步过滤,这一部分内容见call variant中关于snp筛选的一些思考

第五步:变异标准化(可选)

由于VCF文件的灵活性,同一种变异可以通过不同的形式表示, 如下图

2013053-2c0545f3d78aa314.jpg
出自 "Unified representation of genetic variants"

变异标准化按照如下规则对变异位点表示进行简化

  • 尽可能以少的字符表示变异
  • 无等位基因可以标识为长度为0
  • 变异位点必须左对齐

看起来很复杂,其实操作起来很简单

bcftools norm -f $REF SRR1553500_freebayes.vcf  > SRR1553500_freebayes_norm.vcf
# Lines   total/split/realigned/skipped:    493/0/0/0

大部分软件,如GATK, freebayes已经是标准化的结果。

变异检测那么简单吗?

经过简单的实战之后,似乎变异检测是一件非常容易的事情,只要敲几行命令就行了。当然最开始我也是想的,毕竟无知者无畏,但是了解的越多,你就会发现事情并没有那么简单。**大部分基因组相关的DNA序列有一些特性是人类的直觉所不能理解的,因为这需要考虑一些背景。

  1. DNA序列可以非常的长
  2. A/T/G/C能够构建任意组合的DNA序列,因此在完全随机情况下,即使随机分配也能产生各种各样的模式。
  3. DNA序列只有部分会受到随机影响,基本上这部分序列都是有功能的。
  4. 不同物种的不同的DNA序列受到不同的随机性影响
  5. 我们按照实验流程将大片段DNA破碎成小的部分,并尝试通过和参考基因组比对找到原来的位置。只有它依旧和原来的位置非常靠近,才能进一步寻找变异。

因此即便序列和基因组某个序列非常接近,从算法的角度是正确比对,但其实偏离了原来正确的位置,那么从这部分找到的变异也是错误的。那我们有办法解决这个问题嘛?基本上不可能,除非技术进步后,我们可以一次性通读所有序列。当然目前比较常用的方法是找到最优变异,并且那个变异能更好的解释问题,且和每条read中的变异都是一致的。这就是目前变异检测软件常用策略:realignment或probabilistic alignment。

变异注释

变异注释意味着猜测遗传变异(SNP, INDEL, CNY, SV)对基因功能,转录本和蛋白序列以及调控序列的影响。为了对变异进行预测,预测软件需要你提供基因组注释信息,并且注释信息的完善程度决定了预测的准确性。变异预测一般会提供如下结果

  • 变异位点所在基因组注释的位置,是转录本上游,还是编码区,还是非编码RNA
  • 列举出收到影响的转录本和基因
  • 确定变异在蛋白序列上的影响, 如stop_gained(终止密码子提前), missense(错义), stop_lost(终止密码子缺失), frameshift(移码)等
  • 对于人类,还可以和已知的位点进行匹配

这些效应的定义可以在序列本体论查询.

变异注释常用软件有:VEP, snpEFF, AnnoVar, VAAST2. 其中VEP是网页工具http://asia.ensembl.org/Tools/VEP, 使用很方便,可惜支持的物种有限。snpEFF可以说是支持物种最多的工具,这里使用它。

snpEff databases > listing.txt
# 确认物种名
grep -i ebola databases
# 下载
snpEff download ebola_zaire

下载之后还需要检查一下snpEFF提供的注释是否和我们所使用参考基因组一致.

snpEff dump ebola_zaire | less

很不幸的是,snpEFF提供的是基于KJ660346构建的数据库,而我们使用的是AF086833。因此需要重新下载对应的参考基因组重新比对,进行注释。对饮参考基因组的地址为https://www.ncbi.nlm.nih.gov/nuccore/KJ660346.1?report=fasta

snpEff ebola_zaire ${SRR}_freebayes.vcf > ${SRR}_annotated.vcf

最终会生成注释后的VCF文件以及变异位点的描述性报告。

这篇关于biostar handbook(十)|如何进行变异检测的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1088276

相关文章

C#使用SQLite进行大数据量高效处理的代码示例

《C#使用SQLite进行大数据量高效处理的代码示例》在软件开发中,高效处理大数据量是一个常见且具有挑战性的任务,SQLite因其零配置、嵌入式、跨平台的特性,成为许多开发者的首选数据库,本文将深入探... 目录前言准备工作数据实体核心技术批量插入:从乌龟到猎豹的蜕变分页查询:加载百万数据异步处理:拒绝界面

Python使用自带的base64库进行base64编码和解码

《Python使用自带的base64库进行base64编码和解码》在Python中,处理数据的编码和解码是数据传输和存储中非常普遍的需求,其中,Base64是一种常用的编码方案,本文我将详细介绍如何使... 目录引言使用python的base64库进行编码和解码编码函数解码函数Base64编码的应用场景注意

Java进行文件格式校验的方案详解

《Java进行文件格式校验的方案详解》这篇文章主要为大家详细介绍了Java中进行文件格式校验的相关方案,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录一、背景异常现象原因排查用户的无心之过二、解决方案Magandroidic Number判断主流检测库对比Tika的使用区分zip

Java使用Curator进行ZooKeeper操作的详细教程

《Java使用Curator进行ZooKeeper操作的详细教程》ApacheCurator是一个基于ZooKeeper的Java客户端库,它极大地简化了使用ZooKeeper的开发工作,在分布式系统... 目录1、简述2、核心功能2.1 CuratorFramework2.2 Recipes3、示例实践3

基于Flask框架添加多个AI模型的API并进行交互

《基于Flask框架添加多个AI模型的API并进行交互》:本文主要介绍如何基于Flask框架开发AI模型API管理系统,允许用户添加、删除不同AI模型的API密钥,感兴趣的可以了解下... 目录1. 概述2. 后端代码说明2.1 依赖库导入2.2 应用初始化2.3 API 存储字典2.4 路由函数2.5 应

Python使用date模块进行日期处理的终极指南

《Python使用date模块进行日期处理的终极指南》在处理与时间相关的数据时,Python的date模块是开发者最趁手的工具之一,本文将用通俗的语言,结合真实案例,带您掌握date模块的六大核心功能... 目录引言一、date模块的核心功能1.1 日期表示1.2 日期计算1.3 日期比较二、六大常用方法详

Python使用DrissionPage中ChromiumPage进行自动化网页操作

《Python使用DrissionPage中ChromiumPage进行自动化网页操作》DrissionPage作为一款轻量级且功能强大的浏览器自动化库,为开发者提供了丰富的功能支持,本文将使用Dri... 目录前言一、ChromiumPage基础操作1.初始化Drission 和 ChromiumPage

Jackson库进行JSON 序列化时遇到了无限递归(Infinite Recursion)的问题及解决方案

《Jackson库进行JSON序列化时遇到了无限递归(InfiniteRecursion)的问题及解决方案》使用Jackson库进行JSON序列化时遇到了无限递归(InfiniteRecursi... 目录解决方案‌1. 使用 @jsonIgnore 忽略一个方向的引用2. 使用 @JsonManagedR

使用Folium在Python中进行地图可视化的操作指南

《使用Folium在Python中进行地图可视化的操作指南》在数据分析和可视化领域,地图可视化是一项非常重要的技能,它能够帮助我们更直观地理解和展示地理空间数据,Folium是一个基于Python的地... 目录引言一、Folium简介与安装1. Folium简介2. 安装Folium二、基础使用1. 创建

Nginx如何进行流量按比例转发

《Nginx如何进行流量按比例转发》Nginx可以借助split_clients指令或通过weight参数以及Lua脚本实现流量按比例转发,下面小编就为大家介绍一下两种方式具体的操作步骤吧... 目录方式一:借助split_clients指令1. 配置split_clients2. 配置后端服务器组3. 配