2_重测序BSA分析项目结题报告

重测序BSA项目结题报告

客户单位:

报告单位:

联系人:

联系电话:

传真:

报告日期:

项目负责人:审核人:

目录

目录 (1)

1 项目概况 (1)

1.1 合同关键指标 (1)

1.2 项目基本信息 (1)

1.3 项目执行情况 (2)

1.4项目结果概述 (2)

2 项目流程 (3)

2.1 实验流程 (3)

2.2 信息分析流程 (3)

3 生物信息学分析 (5)

3.1 测序数据质控 (5)

3.1.1 原始数据介绍 (5)

3.1.2 碱基测序质量分布 (7)

3.1.3碱基类型分布 (9)

3.1.4 低质量数据过滤 (10)

3.1.5测序数据统计 (10)

3.2 与参考基因组比对统计 (11)

3.2.1 比对结果统计 (11)

3.2.2 插入片段分布统计 (11)

3.2.3 深度分布统计 (12)

3.3 SNP检测与注释 (14)

3.3.1 样品与参考基因组间SNP的检测 (14)

3.3.2 样品之间SNP的检测 (17)

3.3.3 SNP结果注释 (19)

3.4 Small InDel检测与注释 (22)

3.4.1 样品与参考基因组间Small InDel的检测 (22)

3.4.2样品之间Small InDel检测 (22)

3.4.3 Small InDel的注释 (23)

3.5关联分析 (26)

3.5.1 高质量SNP筛选 (26)

3.5.2 SNP-index方法关联结果 (26)

3.5.3 ED方法关联结果 (28)

3.5.4候选区域筛选 (29)

3.6候选区域的功能注释 (30)

3.6.1 候选区域的SNP注释 (30)

3.6.2 候选区域的基因注释 (30)

3.6.2.1候选区域内基因的GO富集分析 (31)

3.6.2.2候选区域内基因的KEGG富集分析 (33)

3.6.2.3候选区域内基因COG分类统计 (36)

3.7结果可视化 (37)

4 数据下载 (38)

4.1结果文件查看说明 (38)

参考文献 (39)

1 项目概况

1.1 合同关键指标

(1) 完成X个样品的重测序，共产生XGbp Clean Data，保证Q30达到80%。

(2) 数据评估：测序数据量，测序数据质量和GC含量的统计。

(3) 与基因组比对：比对效率，基因组覆盖度，基因组覆盖深度统计。

(4) 变异检测和注释：SNP、InDel的检测和注释。

(5) 关联分析：通过计算两个混池间等位基因的基因型频率确定与目标性状关

联的区域。

(6) 候选SNP注释：对关联区域内的SNP注释，包括位置信息和非同义突变信息。

(7) 候选基因注释：对关联区域内的基因进行GO、KEGG、COG、NR、SwissProt

数据库注释。

1.2 项目基本信息

(1)样品信息:

样品编号BMK编号

亲本1（父本）P

亲本2（母本）M

混池1 B1

混池2 B2

注：BMK编号：百迈客对样品的统一编号，实验建库和后续信息分析均使用该编号。

混池规模：30+30；群体类型：F2群体；研究性状：水稻千粒重

(2) 参考基因组信息：

根据水稻的基因组大小以及GC含量等信息，最终选取日本晴水稻基因组作为参考基因组。

具体信息如下所示：

1.测序物种信息：水稻（Oryza sativa），实际基因组大小为419.8 Mb，GC含

量为45.67%；

2.参考物种信息：日本晴水稻（Oryza sativa indica）[1]基因组，组装出的基因

组大小为374.3 Mb，GC含量为43.56%，ScaffoldN50为500Kb，，该基因组

组装到染色体水平，有基因注释信息，版本号为v7.0，下载地址：

http://rapdb.dna.affrc.go.jp/。

1.3 项目执行情况

(1) 样品信息到位时间为2016年XX月XX日。

(2)样品检测合格时间为2016年XX月XX日。

(3) 项目启动时间为2016年XX月XX日。

(4) 项目分析完成时间为2016年XX月XX日。

1.4项目结果概述

(1)数据质控

测序共获得XXGbp数据量，过滤后得到的Clean Reads为XXGbp，Q30达到80%，平均每个样品测序深度X。样品与参考基因组平均比对效率为XX%，平均覆盖深度为X，基因组覆盖度为XX%（至少一个碱基覆盖）。

(2)变异检测

SNP检测：样品P、M之间共获得XX个SNP，其中非同义突变的SNP共XX个；样品B1、B2之间共获得XX个SNP，引起非同义突变的SNP共XX个。

InDel检测：样品P、M之间共获得XX个Small InDel；样品B1、B2之间共获得XX个Small InDel。

(3)关联分析：SNP-index关联算法，共得到XX个与性状相关的侯选区域，总长度为XXbp；ED关联算法，共得到XX个与性状相关的侯选区域，总长度为XXbp，两种方法取交集得到XX个与性状相关的侯选区域，总长度为XXbp。关联区域内包含非同义突变SNP位点的基因共XX个，同义突变SNP位点的基因共XX个。

2 项目流程

2.1 实验流程

实验流程按照Illumina公司提供的标准protocol执行，包括样品检测、文库构建、文库质量检测和上机测序，具体流程如下:

实验流程图

样品检测合格后，用超声破碎的方法将DNA随机打断成350bp的片段，DNA片段经末端修复、3'端加A、加测序接头、纯化、PCR扩增完成测序文库的构建。文库经质检合格后通过Illumina HiSeq TM4000进行测序。

2.2 信息分析流程

信息分析的内容包括：数据质控（去除接头和低质量数据）、与参考基因组比对、变异检测与注释（SNP、InDel）、关联分析、候选SNP及候选基因的注释。

重测序BSA生物信息分析具体流程如下图所示：

重测序BSA生物信息分析流程图

3 生物信息学分析

3.1 测序数据质控

3.1.1 原始数据介绍

高通量测序（如Illunima HiSeq 4000等测序平台）得到的原始图像数据文件，经碱基识别（Base Calling）分析转化为原始测序序列（Sequenced Reads），我们称之为Raw Data或Raw Reads，结果以FASTQ（简称为fq）文件格式存储，其中包含测序序列（Reads）的序列信息以及其对应的测序质量信息。测序样品中真实数据随机截取结果如下：

@HWI-7001455:110: C3B41ACXX:4:1101:1401:2163 1:N:0: TAAGGC CTCTCTCCTATCTTTCCAACCATCTGATAACACCGAACATCCATATTGAGCCCACACTTCTTGATGATCTTTCAATATTTTATGAT + CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJHHHHHHHFFFFFFFFEEEEEE

FASTQ格式文件中每个Read由四行描述，其中第一行以“@”开头，随后为Illumina测序识别符（Sequence Identifiers）和描述文字（选择性部分）；第二行是碱基序列；第三行以“+”开头，随后为Illumina测序识别符（选择性部分）；第四行是对应序列的测序质量。

Illumina测序识别符（Sequence Identifiers）详细信息见如下：

Illumina 测序标识详细信息

HWI-7001455 Unique instrument name

110 Run ID

C3B41ACXX Flowcell ID

4 Flowcell lane

1101 Tile number within the flowcell lane

1401 'x'-coordinate of the cluster within the tile

2163 'y'-coordinate of the cluster within the tile

1 Member of a pair, 1 or

2 (paired-end or mate-pair reads only)

N Y if the read fails filter (read is bad), N otherwise

0 0 when none of the control bits are on, otherwise it is an even number

TAAGGC Index sequence

通过使用第四行中每个字符对应的ASCII值进行计算，即得到对应第二行碱基

的测序质量值。如果测序错误率用e表示，Illunima HiSeq 4000的碱基质量值用Qphred 表示，则有下列关系：

Q phred = -10log10(e)

Illunima Casava 1.8版本测序错误率与测序质量值简明对应关系如下表所示：测序错误率测序质量值对应字符

5% 13 .

1% 20 5

0.1% 30 ?

0.01% 40 I

碱基识别（Base Calling）分析软件：Illunima Casava 1.8版本

测序参数：双端测序(Paired end，PE)

测序序列读长：151bp

3.1.2 碱基测序质量分布

每个碱基测序错误率是通过测序Phred数值（Phred score，Qphred）得到，而Phred数值是在碱基识别（Base Calling）过程通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表所显示：

Phred分值不正确的碱基识别碱基正确识别率

10 1/10 90%

20 1/100 99%

30 1/1000 99.9%

40 1/10000 99.99%

在Hiseq4000测序系统测序时，首先会对文库进行芯片制备，目的是将文库DNA 模板固定到芯片上，在固定DNA模板的过程中，每个DNA分子会形成一个簇，一个簇就是一个测序位点，在进行固定过程中极少量的簇与簇之间物理位置会发生重叠，在测序时，测序软件通过前4个碱基对这些重叠的点进行分析和识别，将这些重叠点位置分开，保证每个点测到的是一个DNA分子，因此测序序列5′端前几个碱基的错误率相对较高。另外测序错误率会随着测序序列（Sequenced Reads）的长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的。因此在进行碱基测序质量分布分析时，样品的碱基质量分布在前4个碱基和后十几个碱基的质量值会低于中间测序碱基，但其质量值都高于Q30，根据质量值和错误率的关系，我们将质量值转换成错误率，绘制错误率分布图如下：

样品P碱基错误率分布

注：横坐标为reads的碱基位置，纵坐标为单碱基错误率，前151bp 为双端测序序列的第一端测序reads的错误率分布情况，后151bp为另一端测序reads 的错误率分布情况。

3.1.3碱基类型分布

碱基类型分布用于检测有无AT、GC分离现象，这种分离现象可能是建库测序过程中差异扩增引起的，直接影响到后续的分析。高通量测序的序列为基因组随机打断后的DNA片段，位点在整个基因组的分布是近似均匀的，同时根据碱基互补配对的原则，A与T和C与G的含量分别是一致的。由于测序仪器本身的局限性，前几个碱基的A/T和C/G含量可能存在着一定波动。

样品各碱基比例分布如下所示：

样品P各碱基比例分布

注：横坐标为reads的碱基位置，纵坐标为碱基所占的比例；不同颜色代表不同的碱基类型，绿色代表碱基G，蓝色代表碱基C，红色代表碱基A，紫色代表碱基T，灰色代表测序中识别不出的碱基N。前151bp为双端测序序列的第一端测序reads的碱基分布，后151bp为另一端测序reads的碱基分布。每个cycle代表测序的每个碱基，如第一cycle即表示该项目所有测序reads在第一个碱基的A、T、G、C、N的分布情况。

该图的结果显示AT、CG碱基基本不发生分离，且曲线较平缓，说明测序结果正常。

3.1.4 低质量数据过滤

测序得到的原始测序序列（Sequenced Reads）或者Raw Reads，里面含有带接头的、低质量的Reads，为了保证信息分析质量，对Raw Reads进行过滤，得到Clean Reads，用于后续信息分析。

数据过滤的主要步骤如下：

(1) 去除带接头（adapter）的reads。

(2) 若一条reads上N（未能确定出具体的碱基类型）的比例大于10%，则过滤掉该Pair-end reads。

(3) 去除低质量reads（质量值Q≤xx的碱基数占整条read的50%以上）。

数据过滤统计结果见下表：

数据过滤统计表

BMKID Raw_Reads Adapter_Related

（%）

Inferior_percent

(%)

Clean_Reads

P

M

B1

B2

注：BMK ID：百迈客对项目样品的统一编号；Raw_Reads：原始测序reads数；Adapter_Related：含接头被过滤的reads比例；Inferior_percent：N含量超过10%的reads和质量值低于10的碱基超过50%的reads比例；Clean_Reads：过滤后剩余的reads数。

3.1.5测序数据统计

各样品测序产出数据评估结果如下表所示：

样品测序数据评估统计

BMKID Raw_Reads Clean_Reads Clean_Base Q30(%)GC(%) P

M

B1

B2

注：BMK ID：百迈客对项目样品的统一编号；Raw_Reads：原始测序reads数目，以四行为一个单位，统计Pair-end 序列的个数；Clean_Reads：过滤后的reads数，计算方法同Raw Reads；Clean_Bases：过滤后的碱基数，Clean Reads 数乘以序列长度；Q30(%)：质量值大于等于30的碱基占总碱基数的百分比；GC(%)：样品GC含量，即G和C类型的碱基占总碱基的百分比。

3.2 与参考基因组比对统计

重测序获得的测序reads需要重新定位到参考基因组上，才可以进行后续变异分

析。bwa[2]软件主要用于二代高通量测序（如Illunima HiSeq 4000等测序平台）得到

的短序列与参考基因组的比对。通过比对定位Clean Reads在参考基因组上的位置，

统计各样品的测序深度、基因组覆盖度等信息，并进行变异的检测。

3.2.1 比对结果统计

样品的比对结果见下表:

比对结果统计

BMK ID Total_reads Mapped(%)Properly_mapped(%) P

M

B1

B2

注：BMK ID：百迈客对项目样品的统一编号；Total_Reads：Total_Reads数，双端分别统计，即read1和read2记为2条reads；Mapped(%)：定位到参考基因组的Clean Reads数占所有Clean Reads数的百分比；Properly mapped：双端测序序列均定位到参考基因组上且距离符合测序片段的长度分布。

此项目样品的平均比对效率均在XX%以上，说明样品测序正常。

3.2.2 插入片段分布统计

通过检测双端序列在参考基因组上的起止位置，可以得到样品DNA打断后得到

的测序片段的实际大小，即插入片段大小(Insert Size)，是信息分析时的一个重要参

数。插入片段大小的分布一般符合正态分布，且只有一个单峰，Insert Size分布图可

以展示各个样品的插入片段的长度分布情况。

每个样品测序数据插入片段大小分布的分析使用picard软件工具包中CollectInsertSizeMetric.jar软件实现。

样品P插入片段分布图

注：横坐标为插入片段长度，纵坐标为其对应的reads数。

由上图可知，插入片段长度分布符合正态分布，说明测序数据文库构建无异常。

3.2.3 深度分布统计

Reads定位到参考基因组后，可以统计参考基因组上碱基的覆盖情况。参考基因组上被reads覆盖到的碱基数占基因组的百分比称为基因组覆盖度；碱基上覆盖的reads数为覆盖深度。

基因组覆盖度可以反映参考基因组上变异检测的完整性，覆盖到的区域越多，可以检测到的变异位点也越多。覆盖度主要受测序深度以及样品与参考基因组亲缘关系远近的影响。

基因组的覆盖深度会影响变异检测的准确性，在覆盖深度较高的区域（非重复序列区），变异检测的准确性也越高。另外，若基因组上碱基的覆盖深度分布较均匀，也说明测序随机性较好。

样品的碱基覆盖深度分布曲线和覆盖度分布曲线见下图：

样品P的深度分布

注：上图反映了测序深度分布的基本情况，横坐标为测序深度，左纵坐标为该深度对应的碱基所占百分比，对应红色曲线，右纵坐标为该深度及以下的碱基所占百分比，对应蓝色曲线。

各样品的平均覆盖深度和各深度对应的基因组覆盖比例如下表所示：

样品覆盖深度和覆盖度比例统计

BMK ID Ave_depth Cov_ratio_1X(%)Cov_ratio_5X(%)Cov_ratio_10X(%) P

M

B1

B2

注：BMK ID：百迈客对项目样品的统一编号；Ave-depth：样品平均覆盖深度；Cov_ratio：覆盖深度在给定深度及以上的碱基数占参考基因组总碱基数的比例。

由上表可知，此项目基因组平均覆盖深度约为X，基因组覆盖度约为XX%（至少覆盖1X）。

根据染色体各位点的覆盖深度情况进行作图，若覆盖深度在染色体上的分布比较均匀，则可以认为测序随机性比较好。样品的染色体覆盖深度分布见下图：

样品P染色体覆盖深度分布图

注：横坐标为染色体位置，纵坐标为染色体上对应位置的覆盖深度取对数（log2）得到的值。

由上图可以看出基因组被覆盖的较均匀，说明测序随机性较好。图上深度不均一的地方可能是由于重复序列、PCR偏好性引起的。

3.3 SNP检测与注释

3.3.1 样品与参考基因组间SNP的检测

SNP的检测主要使用GATK[3]软件工具包实现。根据Clean Reads在参考基因组的定位结果，使用Picard[4]进行去重复（Mark Duplicates）、GATK进行局部重比对(Local Realignment)、碱基质量值校正（Base Recalibration）等预处理，以保证检测得到的SNP准确性，再使用GATK进行单核苷酸多态性（Single Nucleotide Polymorphism，SNP）的检测，过滤，并得到最终的SNP位点集。

主要检测过程如下：

(1) 对于BWA比对得到的结果，使用Picard的Mark Duplicate工具去除重复，屏蔽PCR-duplication的影响。

(2) 使用GATK进行InDel Realignment，即对存在插入缺失比对结果附近的位点进行局部重新比对，校正由于插入缺失引起的比对结果错误。

(3)使用GATK进行碱基质量值再校准（Base Recalibration），对碱基的质量值

进行校正。

(4)使用GATK进行变异检测（variant calling），主要包括SNP和InDel。

(5)对SNP进行严格过滤：snp cluster过滤（5bp内如果有2个SNP则过滤掉），Indel 附近SNP过滤（InDel附近5bp内的SNP过滤掉）；和相邻INDEL过滤（两个InDel距离小于10bp过滤掉）[5]。具体流程可参考GATK官方网站的BestPractice：https://https://www.sodocs.net/doc/db18992400.html,/gatk/guide/best-practices?bpm=DNAseq#variant-dis covery-ovw

变异结果使用vcf文件格式展示。vcf文件包括注释行、标题行和数据行三部分。其中注释行包含文件数据行的INFO和FORMAT列中使用的各种标识符的意义解释，而标题行和数据行包含各样品的变异检测结果信息，格式如下所示：

SNP变异结果信息表

#CHROM POS ID REF ALT QUAL FILTE

R

INFO FORMAT P

Chr1 5634 . G A 140.84 PASS [ANNOTATIONS] GT:AD:DP:GQ:PL 1/1:0,6:6:18:169,18,

Chr1 30071 . A G 141.84 PASS [ANNOTATIONS] GT:AD:DP:GQ:PL 1/1:0,6:6:18:170,18,

Chr1 30478 . C T 95.9 PASS [ANNOTATIONS] GT:AD:DP:GQ:PL 1/1:0,5:5:15:124,15,

Chr1 32667 . A G 91.03 PASS [ANNOTATIONS] GT:AD:DP:GQ:PL 1/1:0,4:4:12:119,12,

各列意义说明如下：

1 CHROM Chr1 参考序列的染色体名称

2 POS 5634 参考序列位点坐标

3 ID . 标识符

4 REF G 参考序列对应位置碱基

5 ALT A SNP位点对应的另外类型的碱基

6 QUAL 140.84 变异位点质量值

7 FILTER PASS 过滤状态

8 INFO [ANNOTATIONS] 位点注释信息

9 FORMAT GT:AD:DP:GQ:PL 基因型信息格式

10 R01 1/1:0,6:6:18:169,18,0 样品的基因型信息

vcf文件的详细说明信息见网页：

https://www.sodocs.net/doc/db18992400.html,/discussion/1268/how-should-i-interpret-vcf-files-pro duced-by-the-gatk

为了确保样本SNP的可信性，对样本检测的SNP的reads支持数，相邻SNP的距离统计累积分布。

SNP质量分布图

注：左边为SNP reads 支持数目累积图，右边为相邻SNP之间的距离累积图。

SNP类型的变异分为转换和颠换两种，同种类型碱基之间突变称为转换（Transition），如嘌呤与嘌呤之间、嘧啶与嘧啶之间的变异，不同类型碱基之间的突变称为颠换（Transversion），如嘌呤与嘧啶之间的变异。一般来说转换比颠换更容易发生，故转换/颠换（Ti/Tv）的比例一般大于1，具体数值和所测物种有关。对于二倍体或者多倍体物种，若同源染色体上的某一SNP位点均为同一种碱基，则该SNP位点称为纯合SNP位点；若同源染色体上的SNP位点包含不同类型的碱基，则该SNP位点称为杂合SNP位点。纯合SNP数量越多，则样品与参考基因组之间差异越大，杂合SNP数量越多，则样品的杂合程度越高，具体结果和样品的材料选择有关。

3.3.2 样品之间SNP的检测

根据样品与参考基因组的比对结果，汇总样品之间所有有差异的变异位点，各样品的SNP列表文件格式如下所示：