‘壹’ TBtools | 基因组重测序数据分析常用上游软件“三兄弟”，聚齐了~
是的，最后，我选择了“一条路走到黑”~ 咱们先不管一些人如何看，东西先整出来。大体如下。
没错。从“原始测序数据”到“SNP Calling”，也就是 从 .fastq 到 .vcf，一套聚齐了。其中涉及到三个目前最常用的软件：
有了这三个插件，所有的“TBtools”用户，可以在 Windows 或者 MacOS 下，点点点，分析自己想分析的数据。当然，时间成本前面有提到，最好是按照至少一天跑一个样品来算。做做 BSAseq 或者少数两三个样品还是可以的。大的群体，那么就可能需要足够好配置的机器....这个暂时不是咱们要考虑的事情。
前两个插件都介绍过了。今天介绍“BCFtools GUI Wrapper”，整体如下。
一次包括两个功能，变异检测和变异过滤。
从界面来看，比较简单。只需要将 MarkDup.PositionSorted.bam 都拖拽放置进去就可以了。
逻辑上需要等一些时间，可以得到初步的 bcf 文件。
“SNP Calling”之后，几乎所有情况下，都会做下一步过滤。具体过滤参数会对后续分析产生一定的影响。此处，默认给了一个相对严格的过滤参数，尽量保留“真阳性”SNP。具体使用方法如下。
简单来说...就是把文件放进去，然后设置一下输出目录，最后点击“Start”即可。
基本完事。三个“玩具”，欢迎“试玩”。我觉得，BSAseq真的可以试试了。因为有了 VCF 之后，剩下的其实就是一个 R包 的事儿。这个就用 R-Plugin 了。而如果你的 PC 真的强，你甚至还可以用搞 GWAS，因为有了 VCF，剩下的不就是？ R 包吗？
‘贰’ 如何在ncbi上查找重测序数据库
以基因的名称（AP1）和物种的名称（如杨树）去检索，如果只需要cDNA序列，数据库选择unigene。 如果需要该基因的基因组序列，选择genome.
‘叁’ 全基因组重测序（whole genome sequencing,WGS）技术简介及常见问题
全基因组重测序（whole genome sequencing,WGS）对全基因组进行测序扫描，能在全基因组范围之内解读所有常见和稀有变异信息，全面检测单核苷酸多态位点（Single Nucleotide Polymorphisms, SNP）、插入缺失位点（Insertion Deletion, InDel）、结构变异（Structure Variation，SV）以及拷贝数变异（Copy Number Variation, CNV）。
按照Illumina公司提供的标准 protocol 执行，包括样品质量检测、文库构建、文库质量检测和文库测序等流程。
全基因组重测序生物信息分析流程
染色体覆盖深度分布图
SNP突变类型分布图
各类型变异在染色体的分布
遗传病致病机制、疾病发生发展、肿瘤分子分型、疾病风险筛查。
Cancer Sci. 2018;109(3):513-522. Nakagawa H, Fujita M. Whole genome sequencing analysis for cancer genomics and precision medicine.
APMIS. 2019;127(5):303-315 Rossing M, Sørensen CS, Ejlertsen B, Nielsen FC. Whole genome sequencing of breast cancer
Leukemia, 2018;32(2):332-342 Burns A, Alsolami R, Becq J, et al. Whole-genome sequencing of chronic lymphocytic leukaemia reveals distinct differences in the mutational landscape between IgHVmut and IgHVunmut subgroups
BMC Genomics. 2019;20(1):231. Published 2019 Mar 20 Li T, Unger ER, Rajeevan MS. Universal human papillomavirus typing by whole genome sequencing following target enrichment: evaluation of assay reprocibility and limit of detection.
‘肆’ 全基因组重测序的概述
全基因组重测序是对已知基因组序列的物种进行不同个体的基因组测序，并在此基础上对个体或群体进行差异性分析。全基因组重测序的个体，通过序列比对，可以找到大量的单核苷酸多态性位点（SNP），插入缺失位点（InDel，Insertion/Deletion）、结构变异位点（SV，Structure Variation）位点和拷贝数变异位点（CNV， number variation)。SBC可以协助客户，通过生物信息手段，分析不同个体基因组间的结构差异， 同时完成注释。
‘伍’ 全基因组重测序的介绍
全基因组重测序是对已知基因组序列的物种进行不同个体的基因组测序，并在此基础上对个体或群体进行差异性分析。SBC将不同梯度插入片段（Insert-Size）的测序文库结合短序列（Short-Reads）、双末端（Paired-End）进行测序，帮助客户在全基因组水平上扫描并检测与重要性状相关的基因序列差异和结构变异，实现遗传进化分析及重要性状候选基因预测。
‘陆’ 全基因组测序，有必要做吗
要不要做，先得了解全基因组测序是在做什么。普通大众是不会关心高端复杂的基因技术和遗传原理的，但是大家可以关心一下，做了全基因测序之后，你得到的是什么。作为已经做过的，我以国内某公司的基因检测报告来跟大家解释一下。
一、祖源分析
不得不说，这个很有用。不同的药物，不同的人服用后会在效果和副作用上有不同的反应。有时候，对于同一病症，有多种治疗药物，我们需要选择对于自身副作用最小、效果最好的。
哇擦，吓死我了，以后再也不喝感冒灵了。
以上就是对于全基因测序主要项目类型的说明，如果经济条件允许，测一下还有可以有的，理论上越早测越好。
‘柒’ 动植物重测序--BSA
针对目标性状，选择 表型极端差异的亲本 构建家系，对该家系 目标性状表型极端的子代 分别混合得到的两个样本池进行全基因组重测序，检测到的两池间DNA差异片段即为候选区域，可进一步定位到目标性状相关的基因或标记。利用该方法，可以快速的对单一性状进行精细定位，加速目标性状的功能基因组研究进程。
对，BSA最主要的要求就是研究的这个性状在亲本和子代中一定是极端差异，这样结果才能更准确。
样本要求：2个亲本+2个子代池（每个子代池≥20个样本）（200mg组织够提DNA即可）
文库类型：WGS文库（有参）
测序数据量：亲本10x/个，子代池≥20x/池（若子代池混样30个，就测30x）
子代池选择：家系群体（F2，RILS，BC，DH等）
1）参考基因组：参考基因组组装质量好坏以及注释是否完全对BSA结果有一定的影响。建议尽量选择组装至染色体水平，且基因组注释较完全的参考基因组。
2）分离群体：理论上任何一对具有相对性状的亲本杂交后产生的分离后代都适用于集群分离分析法。但常用的分离群体为：F2代群体、回交群体、重组自交系、双单倍体群体。
3）亲本选择：亲本间目标性状要求差异显着，尽量选用纯度高的亲本，并进一步通过自交进行纯化；并且，除了关注的1个目标性状外，其他性状尽量一致，无其他性状发生分离。
4）混池个体数：混池的样本数通常选取极端性状（群体总样本数的5～10%内）的个体进行混池。混池中若极端性状个体数太少，可能造成样品数不够，不具备代表性；若样本数过高可能会引入杂合个体，产生干扰。
5）表型统计：表型统计不准确，或表型由多个微效基因控制，可能会引起定位效果不佳。
1、群体完整度对定位效果的影响？
下表参考链接中有经验统计。2 亲本+2 子代池是最优选择，次之为⼀个亲本+2 个子代池；再次为 2 个亲本+1 个子代池； 1 个亲本+1 个子代池；最次为只有子代池，如果只有亲本，那么建议⽼师去做变异检测。若是 mutmap （突变体与野生型杂交），则可提供1个野生型亲本+1个突变型子代池，或者两个子代池。不同群体的index算法原理上稍有不同。
2、如果该物种没有参考基因组怎么办？或者如果样品是两个品种的杂交后代，⽽两个品种都有参考基因组，该选择哪个参考基因组？
如果该物种没有参考基因组，⽽近缘物种有参考基因组，可以选择近缘物种的参考基因组进⾏⽐对分析。如果两个品种都有参考基因组，则要根据样品的表型和其他数据选择⼀个与样品差异较⼩的品种作为参考基因组。
3、如何构建混池？混池规模多⼤能够满⾜要求？
通过表型调查，筛选具有极端表型的个体构建两个混池，质量性状（1： 2： 1） 具有主效位点的数量性状（正态分布）。表型统计不准确，或由多个微效基因控制，会影响定位效果。混池个体性状不极端，存在较多中间型，同样会导致定位结果不佳。⼀般选择群体的 5%到 10%个体构建极端混池，建议极端样品在 30+30 以上。通常在样本量等情况允许的条件下越⼤越好，能够更进⼀步的缩⼩候选区域。具体需根据物种、样本情况等设计方案。
4、简化基因组测序适合做 BSA 性状定位吗？
针对微效多基因控制的数量性状，与⽬标性状相关的位点很多，简化基因组或许有幸捕获到少数位点，但绝⼤部分会被遗漏，这对后续的研究很不利。针对质量性状或由少数主效基因控制的数量性状，简化基因组 BSA 的⽅法风险更⾼。诺⽲致源不做简化基因组的 BSA性状定位。
5、能否保证定位到的候选区域的⼤⼩和候选基因的个数？
候选区域的⼤⼩和候选基因个数的多少，与群体的⼤⼩、亲本材料差异的⼤⼩、⽬标性状特点、测序深度的多少、分析物种的基因组⽔平等多项因素相关，因此不能给出统⼀的标准，但可以参考已完成的项⽬经验进⾏估计。
6、BSA 能否检测 SV、 CNV？
不建议⽤ BSA 的结果检测 SV、 CNV，第⼀、测序深度不够；第⼆、 BSA ⼦代是混池样本，混池样本做⼤结构变异检测误差较⼤，且没有分析意义； 想做⼤结构变异检测，建议⽤新技术做， Bionano、三代、 10X Genomics 或者测 30X ⼆代做分析。
该方法参考于2019年7月发表在《自然-通讯》（ Nature Communications ）杂志上发表了题为 Dissecting a heterotic gene through GradedPool-Seq mapping informs a rice-improvement strategy 的研究论文。
文中将F2群体不同表型均进行了混池测序，100-150个个体/池，个体平均测序深度0.5x，测序数据比对，变异检测，SNP标记筛选，统计检验-背景降噪-定位。这个降噪后分析结果看上去要比SNP-index的结果简单很多，精度也会高。如果样本足够，经费阔绰，建议可以用这个分析方法尝试一下。
参考学习：
1、 BSA你真得懂了吗?
2、 BSA的原理
3、 GPS的论文 报道
‘捌’ gwas和全基因组重测序的区别
基于第二代高通量测序技术，对于有参考序列的物种，针对不同的真菌菌株，可通过全基因组重测序的方法获得全基因组范围内完整的变异信息，讨论群体的遗传结构、影响群体遗传平衡的因素以及物种形成的机制，定位重要性状位点，为后续分子育种打下坚实基础。同时，通过全基因组大样本重测序对真菌重要菌株进行全基因组的基因型鉴定，并与关注的表型数据进行全基因组关联分析（GWAS），找出与关注表型相关的SNP位点，定位性状相关基因。随着测序成本降低和拥有参考基因组序列的物种增多，基因组重测序也成为育种研究中迅速有效的方法之一，在全基因组水平扫描并检测出与重要性状相关的变异位点，具有重大的科研价值和产业价值。
近日，Nature Genetics发表的一篇文章就充分利用了微生物基因组测序与以全基因组重测序为基础的全基因组关联分析结合的方法，揭示了裂殖酵母遗传与表型多样性之间的联系。研究者选取裂殖酵母Schizosaccharomyces pombe作为研究对象，在全球20个国家范围内收集了时间跨度为100年的161个野生株系的S.Pombe，进行了全基因组测序，推测裂殖酵母在公元前340年开始广泛大量出现，祖先种到达美洲的时间为公园1623年。后续研究者又选取223个菌种进行全基因组关联分析，发现至少89个性状表现出一个关联。每个性状最显着的检测到的变异可以解释平均22%的表型差异，且indel的影响比SNP更大。
‘玖’ 基因组重测序SNP_calling
基因组重测序是对已知基因组序列的物种进行不同个体的基因组测序，并在此基础上对个体或群体进行差异性分析。
基于全基因组重测序技术，人们期望快速进行资源普查筛选，寻找到大量遗传变异，实现遗传进化分析及重要性状候选基因的预测。随着测序成本降低和拥有参考基因组序列物种增多，全基因组重测序成为动植物育种和群体进化研究迅速有效的方法。
• SNP 强调群体在这个位点有多样性。
• Genotype 强调个体在这个位点是哪种基因型。SNV属于genotype的范畴，表示这个genotype 是个突变，是第一无二的。
• 首先利用 bwa , samtools , picard 对参考基因组ref.fa建立索引，随后使用基因组比对软件 bwa 将fastq文件比对回参考基因组，生成sam文件(这里使用例子数据由于小于10M，所以使用bwa的-is参数，当参考基因组很大的时候，使用bwtsw参数)
• 使用 picard 对sam文件排序后生成bam文件(这里如果使用samtools的话会删除重复位点，而picard只会标记重复并不会删除)。之后使用 samtools 对bam文件建立索引，利用已经注释过的indel.vcf文件对indel周围进行realignment，这个操作有两个步骤，第一步生成需要进行realignment的位置的信息，第二步 才是对这些位置进行realignment，最后生成一个重排好的BAM。
• 接下来对mapping得到的bam文件进行矫正，消除偏态，去除假阳性，这里主要针对测序质量不是特别好的数据，可以根据fastqc结果对测序结果进行解读，如果比较好的话可以跳过。
• 使用 HaplotypeCaller 发现变异并生成gvcf文件。随后使用 CombineGVCFs 将不同文件整合， GenotypeGVCFs 函数输出vcf文件。
• 下一步使用 VariantRecalibrator 校准，并根据可视化结果确定最佳参数，确定参数之后再运行一次 VariantRecalibrator 得出结果
• 如何确定参数：(1) an注释信息选择:从图中可以看出，这个模拟结果可以很好的将真实的变异位点和假阳性变异位点分开，形成了明显的界限， 也就是说，如果一个变异位点的这两个注释值，只要有一个落在了界限之外，就会被过滤掉。最主要的是要 看右边两个图片，只要能很好的区分开novel和known以及filtered和retained就可以。其实在如何选择注释值 存在一定得主观性，因此，在做VariantRecalibrator时可以做两次，第一次尽可能的多的选择这些注释值， 第一遍跑完之后，选择几个区分好的，再做一次VariantRecalibrator，然后再做ApplyRecalibration。 (2) tranche值的设定:如果这个值设定的比较高的话，那么最后留下来的变异位点就会多，但同时假阳性的位点也会相应增加； run两遍VariantRecalibrator，第一遍的时候多写几个阈值，看哪个阈值结果好，选择一个最好的阈值，再run一遍 VariantRecalibrator。 区分标准： 1. 看结果中已知变异位点与新发现变异位点之间的比例，这个比例不要太大，大多数新发现的变异都是假阳性，如果太多的 话，可能假阳性的比例就比较大； 2. 看保留的变异数目，这个就要根据具体的需求进行选择了。 – 3. 看TI/TV值，对于人类全基因组，这个值应该在2.15左右，对于外显子组，这个值应该在3.2左右，不要太小或太大，越接近 这个数值越好，这个值如果太小，说明可能存在比较多的假阳性。 千人中所选择的tranche值是99，仅供参考。
• 最后使用 annovar 对vcf进行注释并可视化结果
对于基因组
对于转录组
基因组使用bwa
转录组使用STAR 2-pass，第一步比对
根据比对结果重新建立索引
进行第二次比对
基因组使用 SortSam.jar
转录组使用 AddOrReplaceReadGroups.jar
由于 STAR 软件使用的 MAPQ 标准与 GATK 不同，而且比对时会有 reads 的片段落到内含子区间，需要进行一步 MAPQ 同步和 reads 剪切，使用 GATK 专为 RNA-seq 应用开发的工具 SplitNCigarReads 进行操纵，它会将落在内含子区间的 reads 片段直接切除，并对 MAPQ 进行调整。DNA 测序的重测序应用中也有序列比对软件的 MAPQ 与 GATK 无法直接对接的情况，需要进行调整。
之后就是可选的 Indel Realignment，对已知的 indel 区域附近的 reads 重新比对，可以稍微提高 indel 检测的真阳性率，如果对于某些物种没有已知indel区域的话官网推荐可以跳过这步直接生成vcf文件，之后使用第一次生成的vcf文件作为已知indel再走一遍流程。
主要是针对测序质量不是特别高的数据
随后可以使用 annovar 对vcf文件进行注释，可以使用 convert2annovar.pl 进行格式转换。比如将 VCF 和 pileup 格式的文件转换为 annovar 的输入格式。
‘拾’ 全基因组重测序的生物信息分析内容
1．数据量产出
总碱基数量、Totally mapped reads、Uniquely mapped reads统计，测序深度分析。
2．一致性序列组装
与参考基因组序列（Reference genome sequence）的比对分析，利用贝叶斯统计模型检测出每个碱基位点的最大可能性基因型，并组装出该个体基因组的一致序列。
3．SNP检测及在基因组中的分布
提取全基因组中所有多态性位点，结合质量值、测序深度、重复性等因素作进一步的过滤筛选，最终得到可信度高的SNP数据集。并根据参考基因组序列对检测到的变异进行注释。
4．InDel检测及在基因组的分布
在进行mapping的过程中，进行容Gap的比对并检测可信的Short InDel。在检测过程中，Gap的长度为1~5个碱基。
5．Structure Variation检测及在基因组中的分布
目前SBC能够检测到的结构变异类型主要有：插入、缺失、复制、倒位、易位等。根据测序个体序列与参考基因组序列比对分析结果，检测全基因组水平的结构变异并对检测到的变异进进行注释。
全基因组重测序生物信息学分析流程