两分钟教你完成变异检测

在前面的教程中，小果已经给大家分享了质量评估、序列比对等操作，相信小伙伴们也是收获满满。那么今天，小果就带大家迈入新的阶段，揭开变异检测的神秘面纱~像往常一样小果奉行着“先跑流程，后看原理”的流程，这一次先给大家详细地介绍变异检测的全部流程，下一次再给大家分享变异检测的原理。好啦，废话不多说，我们直接开始~

操作代码

像之前一样，小果开门见山先把操作所需要的全部代码先提供给大家，帮助急需操作的小伙伴们解燃眉之急：

# 创建索引
samtools index ./02_align/test.clean.sort.uniq.bam
# 选用freebayes进行变异检测
freebayes -m 30 -q 20 -f ./00_ref/MGI358.SNP.fa -@ ./00_ref/alleles_all.vcf.gz -t ./00_ref/target.358.SE50.subSNP.bed –report-all-haplotype-alleles ./02_align/test.clean.sort.uniq.bam > ./03_SNPCalling/test.clean.SNP.vcf
# 将vcf文件的路径写入vcf.list文件
echo “test ./03_SNPCalling/test.clean.SNP.vcf 1” > ./03_SNPCalling/vcf.list
# 使用vcf2geno_free.pl进行基因型转换
perl vcf2geno_free.pl -in ./03_SNPCalling/vcf.list -dir ./03_SNPCalling/ -std
./00_ref/alleles_all.vcf.gz -summary

以上就是变异检测所需要的全部代码啦，是不是很简单呢~如果小伙伴们平时在生信分析的操作过程中遇到困难，也欢迎大家使用小果开发的生信工具平台http://www.biocloudservice.com/home.html哦~

接下来，小果会给大家详细地介绍用到的工具和每一步代码的参数以及对结果进行分析，小伙伴们千万不要走开哦~

工具介绍

FreeBayes

介绍

FreeBayes是一个用于基因组变异检测和单核苷酸多态性（Single Nucleotide Polymorphism，SNP）分析的开源软件工具，它能够从下一代测序数据中准确地鉴定出SNP、小片段插入/删除（indel）和结构变异等基因组变异。

FreeBayes的工作原理基于贝叶斯定理和单样本比对。它利用已比对的测序数据和参考基因组来计算每个位点上的变异概率，并通过统计模型对突变进行筛选和过滤。与其他常见的SNP鉴定工具相比，FreeBayes在处理高杂合度和复杂变异的情况下表现较好。

FreeBayes具有以下特点：

1. 高灵敏度：FreeBayes在鉴定低频率变异（如突变等）方面表现优秀，能够准确地识别出较罕见的变异。
2. 多样性支持：FreeBayes可以在不同类型的样本（如人类、动植物等）上使用，并且适用于不同种类的测序数据（如Illumina、PacBio、ONT等）。
3. 区域感知：FreeBayes可以针对特定感兴趣区域进行变异检测，从而提高分析的效率和准确性。
4. 并行处理：FreeBayes支持并行计算，可以利用多个处理器或计算节点进行高效的变异检测。
5. 开源免费：FreeBayes是一个开源软件，用户可以自由地使用、修改和分发。

单倍型（单倍体基因型的简称）在遗传学上是指在同一染色体上进行共同遗传的多个基因座（即位点）上等位基因的组合；通俗的说法就是若干个决定同一性状的紧密连锁的基因构成的基因型。 按照某一指定基因座上基因重组发生的数量，单倍型可以指至少两个基因座甚至整个染色体。

FreeBayes也可以称为贝叶斯遗传变异检测器，其主要目标是通过分析单倍型来寻找小型多态性事件（例如SNP（单核苷酸多态性）、indels（插入和删除事件）、MNP（多核苷酸多态性））和复杂事件（复合插入和替换事件）。它通过使用带有编码质量分数（Phred+33，具体信息小果在第一期中分享过）的BAM文件中的reads比对结果，对任意数量的个体和参考基因组（FASTA格式）进行比对，从而确定在参考基因组的每个位置上群体最可能的基因型组合。这种方法可以帮助我们了解不同个体的基因组差异，并揭示遗传变异与疾病等重要生物学现象之间的关联。

关于小型多态性事件（即基因组中的较小规模的变异情况）的更多详细信息小果也会在后面变异检测的原理详解分享中给大家带来，小伙伴们记得持续关注哦~

调用命令：

freebayes -f ref.fa aln.bam \>var.vcf

常用参数

小果在这里还为大家整理了FreeBayes常用的一些参数及其功能：

1. –fasta-reference：该参数用于指定参考基因组的FASTA文件。参考基因组是进行比对和变异检测的基础，通过指定FASTA文件可以告知FreeBayes使用哪个基因组作为参考进行分析。
2. –bam：该参数用于指定比对测序数据的BAM文件。BAM文件是存储测序数据比对结果的二进制格式文件，包含了每个测序Read的比对位置、质量等信息。通过指定BAM文件，FreeBayes可以从中读取比对信息，并进行变异检测。
3. –min-alternate-count：该参数用于设定最小的非参考等位基因计数阈值。在单个个体中，如果某个位置的非参考等位基因的观察数低于这个阈值，那么该位置的等位基因将被过滤掉。默认情况下，对于双倍体样本，该值为2。
4. –min-alternate-fraction：该参数用于设定最小的非参考等位基因频率阈值。如果某个位置的非参考等位基因在个体中的频率低于这个阈值，那么该位置的等位基因将被过滤掉。
5. –ploidy：该参数用于指定样本的倍性（即染色体的拷贝数）。默认情况下，FreeBayes假设样本为二倍体（diploid），可以根据实际情况进行修改。
6. –genotype-qualities：该参数用于计算基因型的似然比，并将结果输出到VCF文件中。基因型质量分数（genotype qualities）用于衡量基因型推断的可靠性和准确性。
7. –min-mapping-quality：该参数用于设定最小的比对质量阈值。如果某个比对的质量低于这个阈值，那么该比对将被过滤掉。这个参数在只需要进行变异检测而不需要额外需求时，可以使用默认值（对于双倍体生物来说）。
8. –min-base-quality：该参数用于设定最小的碱基质量阈值。如果某个碱基的质量低于这个阈值，那么该碱基将被过滤掉。
9. –min-supporting-mapping-quality：该参数用于设定最小的支持比对质量阈值。如果某个比对的质量低于这个阈值，那么该比对将不会被视为支持变异。
10. –min-supporting-base-quality：该参数用于设定最小的支持碱基质量阈值。如果某个碱基的质量低于这个阈值，那么该碱基将不会被视为支持变异。
11. –min-indel-support：该参数用于设定最小的插入/删除的支持阈值。如果某个插入或删除的支持数低于这个阈值，那么该插入/删除将被过滤掉。
12. –min-strand-bias-pvalue：该参数用于设定最小偏向某一链的偏倚检验P值阈值。如果位点的偏倚检验P值高于这个阈值，那么该位点将被过滤掉。
13. –exclude-unobserved-genotypes：该参数用于排除未观察到的基因型，即在分析结果中不包含未观察到的基因型。
14. –genotype-qualities：该参数用于计算基因型的边际概率，并输出到VCF文件中。基因型的边际概率是基于观察到的测序数据和基因型质量等信息进行推断的结果。

这些仅是FreeBayes可用参数中的一部分，还有其他一些高级参数和过滤选项，用于更精细地控制变异检测过程。在使用FreeBayes时，大家可以通过查看相关文档或运行 ‘freebayes –help’ 命令获取更详细的参数说明和用法示例。

vcf2geno_free.pl

vcf2geno_free.pl是一个用于将VCF文件转换为基因型文件的Perl脚本，它可以将VCF文件中的多个样本的变异位点的基因型信息提取并输出到一个基因型文件中，方便后续的遗传分析和统计分析。

使用vcf2geno_free.pl时，需要指定输入的VCF文件和输出的基因型文件。

例如，以下命令可以将input.vcf文件转换为output.geno文件：

perl vcf2geno_free.pl –input input.vcf –output output.geno

小伙伴们使用之前需要确保已经安装了Perl环境，并且安装了所需的依赖模块哦。

流程详解

上面的代码可能有的小伙伴不是很理解，小果在这里给小伙伴们详细地逐句解释一下：

samtools index ./02_align/test.clean.sort.uniq.bam

• samtools 是一个常用的处理 SAM和 BAM文件的工具，它有许多用于对测序数据进行操作和分析的命令。
• index 是 samtools 的一个子命令，用于为 BAM 文件创建索引，索引文件（BAI文件）是一种辅助文件，可以加速对 BAM 文件的访问。
• ./02_align/test.clean.sort.uniq.bam 是待索引的 BAM 文件的路径和文件名。

总结：这段代码是为经过坐标排序的.BAM文件创建索引，生成以.bai为后缀的索引文件，以便快速的访问bam文件

freebayes -m 30 -q 20 -f ./00_ref/MGI358.SNP.fa -@ ./00_ref/alleles_all.vcf.gz -t ./00_ref/target.358.SE50.subSNP.bed –report-all-haplotype-alleles ./02_align/test.clean.sort.uniq.bam > ./03_SNPCalling/test.clean.SNP.vcf

代码中的”-m 30″表示设置最小的读取深度为30，即只考虑那些至少有30个支持序列的位点。

“-q 20″表示设置最低质量阈值为20，即只考虑序列质量得分大于20的碱基。

“-f ./00_ref/MGI358.SNP.fa”指定参考基因组文件为”./00_ref/MGI358.SNP.fa”，即使用该文件作为比对的参考序列。

“-@ ./00_ref/alleles_all.vcf.gz”表示使用”./00_ref/alleles_all.vcf.gz”作为已知的SNP位点注释文件，可能用于帮助准确调用SNP。

“-t ./00_ref/target.358.SE50.subSNP.bed”指定一个目标区域文件，限制分析范围为”./00_ref/target.358.SE50.subSNP.bed”中定义的区域。

“–report-all-haplotype-alleles”选项可能是用于报告所有haplotype的等位基因信息。

“./02_align/test.clean.sort.uniq.bam”是待处理的BAM文件路径，这里将其作为输入数据进行SNP calling。

“> ./03_SNPCalling/test.clean.SNP.vcf”将程序的输出结果重定向到名为”./03_SNPCalling/test.clean.SNP.vcf”的文件中，即将SNP calling 结果保存到该文件中。

总结：这段代码的目的是运行”freebayes”程序进行SNP calling，它会基于给定的参考基因组、已知的SNP位点信息和目标区域文件，对输入的BAM文件中的序列进行分析，调用和报告在目标区域内的SNP位点，并将结果保存到指定的VCF文件中。

echo “test ./03_SNPCalling/test.clean.SNP.vcf 1” > ./03_SNPCalling/vcf.list

该行文本的含义是在当前目录下的”03_SNPCalling”文件夹中创建一个名为”vcf.list”的列表文件，并写入”test ./03_SNPCalling/test.clean.SNP.vcf 1″这行内容，vcf.list文件是一个描述VCF文件的列表文件，列表中包含了要处理的VCF文件的路径和其他相关信息。

perl vcf2geno_free.pl -in ./03_SNPCalling/vcf.list -dir ./03_SNPCalling/ -std
./00_ref/alleles_all.vcf.gz -summary

• perl vcf2geno_free.pl 是调用该 Perl 脚本。
• -in ./03_SNPCalling/vcf.list 指定了输入文件列表的路径和文件名。
• -dir ./03_SNPCalling/ 指定了输出目录的路径，基因型数据将被写入该目录。
• -std ./00_ref/alleles_all.vcf.gz 指定了参考标准文件的路径和文件名。这个文件可能是包含参考基因组的 VCF 文件或其他形式的参考文件。
• -summary 是一个选项，表示生成基因型数据的摘要统计信息。

怎么样，小伙伴们是不是已经对变异检测的流程理解深刻了呢？如果小伙伴们平时在生信分析的操作过程中遇到困难，欢迎大家使用小果开发的生信工具平台http://www.biocloudservice.com/home.html。大家在新接触一个知识的时候，与其先花费大量时间死啃知识点，不如先利用好工具先自己上手接触流程，在跑完一遍全流程后再返回去理解知识点，相信可以更好更快地理解，达到事半功倍的效果！

为了让小伙伴们对变异检测更加精通，成为真正的大神，小果在下两期还会继续分享如何分析变异检测得到的VCF文件和genotype文件以及变异检测的原理详解，小伙伴们记得持续关注哦~那么我们下次再见啦，拜拜┏(＾0＾)┛~