2重測(cè)序BSA分析項(xiàng)目結(jié)題報(bào)告

1、重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告客戶單位：報(bào)告單位：聯(lián)系人：聯(lián)系電話:傳真：報(bào)告日期：項(xiàng)目負(fù)責(zé)人：審核人:重測(cè)序 BSA 項(xiàng)目結(jié)題報(bào)告目錄目錄 11 項(xiàng)目概況 11.1 合同關(guān)鍵指標(biāo) 11.2 項(xiàng)目基本信息 11.3 項(xiàng)目執(zhí)行情況 21.4項(xiàng)目結(jié)果概述 22 項(xiàng)目流程 32.1 實(shí)驗(yàn)流程 32.2 信息分析流程 33 生物信息學(xué)分析 53.1 測(cè)序數(shù)據(jù)質(zhì)控 53.1.1 原始數(shù)據(jù)介紹 53.1.2 堿基測(cè)序質(zhì)量分布 73.1.3堿基類型分布 93.1.4 低質(zhì)量數(shù)據(jù)過(guò)濾 103.1.5測(cè)序數(shù)據(jù)統(tǒng)計(jì) 103.2 與參考基因組比對(duì)統(tǒng)計(jì) 113.2.1 比對(duì)結(jié)果統(tǒng)計(jì) 113.2.2 插入片段分布統(tǒng)計(jì) 113

2、.2.3 深度分布統(tǒng)計(jì) 123.3 SNP 檢測(cè)與注釋 14331樣品與參考基因組間SNP的檢測(cè) 14332樣品之間SNP的檢測(cè) 173.3.3 SNP結(jié)果注釋193.4 Small In Del 檢測(cè)與注釋223.4.1 樣品與參考基因組間Small InDel 的檢測(cè) 223.4.2樣品之間 Small InDel檢測(cè) 22343 Small In Del 的注釋233.5 關(guān)聯(lián)分析 263.5.1高質(zhì)量SNP篩選263.5.2 SNP-index方法關(guān)聯(lián)結(jié)果 263.5.3 ED方法關(guān)聯(lián)結(jié)果 283.5.4候選區(qū)域篩選 293.6 候選區(qū)域的功能注釋 303.6.1候選區(qū)域的SNP注釋3

3、03.6.2 候選區(qū)域的基因注釋 30 候選區(qū)域內(nèi)基因的 GO 富集分析 3候選區(qū)域內(nèi)基因的KEGG富集分析3候選區(qū)域內(nèi)基因COG分類統(tǒng)計(jì) 363.7結(jié)果可視化374 數(shù)據(jù)下載 384.1 結(jié)果文件查看說(shuō)明 38參考文獻(xiàn) 391重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告1項(xiàng)目概況1.1合同關(guān)鍵指標(biāo)(1) 完成X個(gè)樣品的重測(cè)序，共產(chǎn)生XGbp Clean Data,保證Q30達(dá)到80%。(2) 數(shù)據(jù)評(píng)估：測(cè)序數(shù)據(jù)量，測(cè)序數(shù)據(jù)質(zhì)量和 GC含量的統(tǒng)計(jì)。(3) 與基因組比對(duì)：比對(duì)效率，基因組覆蓋度，基因組覆蓋深度統(tǒng)計(jì)。變異檢測(cè)和注釋：SNP、In Del的檢測(cè)和注釋。(5)

4、 關(guān)聯(lián)分析：通過(guò)計(jì)算兩個(gè)混池間等位基因的基因型頻率確定與目標(biāo)性狀關(guān) 聯(lián)的區(qū)域。(6) 候選SNP注釋：對(duì)關(guān)聯(lián)區(qū)域內(nèi)的SNP注釋，包括位置信息和非同義突變信息。(7) 候選基因注釋：對(duì)關(guān)聯(lián)區(qū)域內(nèi)的基因進(jìn)行 GO、KEGG、COG、NR、SwissProt 數(shù)據(jù)庫(kù)注釋。1.2項(xiàng)目基本信息(1)樣品信息：樣品編號(hào)BMK編號(hào)親本1 (父本)P親本2 (母本)M混池1B1混池2B2注：BMK編號(hào)：百邁客對(duì)樣品的統(tǒng)一編號(hào)，實(shí)驗(yàn)建庫(kù)和后續(xù)信息分析均使用該編號(hào)?；斐匾?guī)模：30+30；群體類型：F2群體；研究性狀：水稻千粒重(2) 參考基因組信息：根據(jù)水稻的基因組大小以及 GC含量等信息，最終選取日本晴水稻基因

5、組作為 參考基因組。具體信息如下所示：1. 測(cè)序物種信息：水稻(Oryza sativa)，實(shí)際基因組大小為419.8 Mb，GC含 量為 45.67%;2. 參考物種信息：日本晴水稻(Oryza sativaindica)基因組，組裝出的基因組大小為374.3 Mb，GC含量為43.56%, ScaffoldN50為500Kb,該基因組 組裝到染色體水平，有基因注釋信息，版本號(hào) 為v7.0,下載地址： http:/rapdb.d na.affrc.go.jp/。1.3項(xiàng)目執(zhí)行情況(1) 樣品信息到位時(shí)間為2016年XX月XX日。(2) 樣品檢測(cè)合格時(shí)間為2016年XX月XX日。(3) 項(xiàng)目啟

6、動(dòng)時(shí)間為2016年XX月XX日。 項(xiàng)目分析完成時(shí)間為2016年XX月XX日。1.4項(xiàng)目結(jié)果概述(1) 數(shù)據(jù)質(zhì)控測(cè)序共獲得XXGbp數(shù)據(jù)量，過(guò)濾后得到的Clean Read為XXGbp , Q30達(dá)到80%, 平均每個(gè)樣品測(cè)序深度X。樣品與參考基因組平均比對(duì)效率為 XX%，平均覆蓋深度 為X，基因組覆蓋度為XX% (至少一個(gè)堿基覆蓋)。(2) 變異檢測(cè)SNP檢測(cè)：樣品P、M之間共獲得XX個(gè)SNP，其中非同義突變的SNP共XX個(gè)； 樣品B1、B2之間共獲得XX個(gè)SNP，弓I起非同義突變的SNP共XX個(gè)。In Del檢測(cè)：樣品P、M之間共獲得XX個(gè)Small I nDel;樣品B1、B2之間共獲得

7、XX 個(gè) Small In Del。(3) 關(guān)聯(lián)分析：SNP-index關(guān)聯(lián)算法，共得到XX個(gè)與性狀相關(guān)的侯選區(qū)域，總長(zhǎng) 度為XXbp ; ED關(guān)聯(lián)算法，共得到XX個(gè)與性狀相關(guān)的侯選區(qū)域，總長(zhǎng)度為 XXbp， 兩種方法取交集得到XX個(gè)與性狀相關(guān)的侯選區(qū)域，總長(zhǎng)度為 XXbp。關(guān)聯(lián)區(qū)域內(nèi)包 含非同義突變SNP位點(diǎn)的基因共XX個(gè)，同義突變SNP位點(diǎn)的基因共XX個(gè)。92項(xiàng)目流程2.1實(shí)驗(yàn)流程實(shí)驗(yàn)流程按照l(shuí)llumina公司提供的標(biāo)準(zhǔn)protocol執(zhí)行，包括樣品檢測(cè)、文庫(kù)構(gòu)建、文庫(kù)質(zhì)量檢測(cè)和上機(jī)測(cè)序，具體流程如下Genemle Cluslers片啟選掙和PCR富年文庫(kù)純化去接頭污毀A2 R*gnrd

8、le ClutLert diid Sitqufic* Piir*d End實(shí)驗(yàn)流程圖樣品檢測(cè)合格后，用超聲破碎的方法將 DNA隨機(jī)打斷成350bp的片段，DNA片 段經(jīng)末端修復(fù)、3端加A、加測(cè)序接頭、純化、PCR擴(kuò)增完成測(cè)序文庫(kù)的構(gòu)建。文庫(kù) 經(jīng)質(zhì)檢合格后通過(guò)Illumina HiSeqTM4000進(jìn)行測(cè)序。2.2信息分析流程信息分析的內(nèi)容包括：數(shù)據(jù)質(zhì)控（去除接頭和低質(zhì)量數(shù)據(jù)）、與參考基因組比 對(duì)、變異檢測(cè)與注釋（SNP、In Del）、關(guān)聯(lián)分析、候選SNP及候選基因的注釋。重測(cè)序BSA生物信息分析具體流程如下圖所示:重測(cè)序BSA生物信息分析流程圖3生物信息學(xué)分析3.1測(cè)序數(shù)據(jù)質(zhì)控3.1.1原始

9、數(shù)據(jù)介紹高通量測(cè)序（如lllunima HiSeq 4000等測(cè)序平臺(tái)）得到的原始圖像數(shù)據(jù)文件，經(jīng) 堿基識(shí)別（Base Calling）分析轉(zhuǎn)化為原始測(cè)序序列（Sequeneed Reads，我們稱之 為Raw Data或Raw Reads結(jié)果以FASTQ （簡(jiǎn)稱為fq）文件格式存儲(chǔ)，其中包含測(cè)序 序列（Reads的序列信息以及其對(duì)應(yīng)的測(cè)序質(zhì)量信息。測(cè)序樣品中真實(shí)數(shù)據(jù)隨機(jī)截取結(jié)果如下：HWI-7001455:110: C3B41ACXX:4:1101:1401:2163 1:N:0: TAAGGCCTCTCTCCTATCTTTCCAACCATCTGATAACACCGAACATCCATATTGA

10、GCCCACACTTCTTGATGATCTTTCAATATTTTATGATCCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJHHHHHHHFFFFFFFFEEEEEEFASTQ格式文件中每個(gè)Read由四行描述，其中第一行以“ 開(kāi)頭，隨后為Mu mi na測(cè)序識(shí)別符（Seque nee Ide ntifiers）和描述文字（選擇性部分）；第二行是 堿基序列；第三行以“ +”頭，隨后為Illumina測(cè)序識(shí)別符（選擇性部分）；第四行 是對(duì)應(yīng)序列的測(cè)序質(zhì)量。Illumina測(cè)序識(shí)別符（Sequenee Identifie

11、r詳細(xì)信息見(jiàn)如下：Illumina測(cè)序標(biāo)識(shí)詳細(xì)信息HWI-7001455110Unique instrument nameRun IDFlowcell IDC3B41ACXXFlowcell lane110114012163TAAGGCTile number within the flowcell lanex-coordinate of the cluster within the tile y-coordinate of the cluster within the tileMember of a pair, 1 or 2 (paired-end or mate-pair reads onl

12、y)Y if the read fails filter (read is bad), N otherwise0 when none of the control bits are on, otherwise it is an even numberIndex sequence通過(guò)使用第四行中每個(gè)字符對(duì)應(yīng)的 ASCII值進(jìn)行計(jì)算，即得到對(duì)應(yīng)第二行堿基 的測(cè)序質(zhì)量值。如果測(cè)序錯(cuò)誤率用e表示，lllunima HiSeq 4000的堿基質(zhì)量值用Qphred 表示，則有下列關(guān)系：Qphred = -10lOg10(e)Mun ima Casava 1.8版本測(cè)序錯(cuò)誤率與測(cè)序質(zhì)量值簡(jiǎn)明對(duì)應(yīng)關(guān)系如下表所

13、示：測(cè)序錯(cuò)誤率測(cè)序質(zhì)量值對(duì)應(yīng)字符5%131%2050.1%30?0.01%40I堿基識(shí)別(Base Calling)分析軟件：Illunima Casava 1.8版本測(cè)序參數(shù)：雙端測(cè)序(Paired end, PE)測(cè)序序列讀長(zhǎng)：151bp3.1.2堿基測(cè)序質(zhì)量分布每個(gè)堿基測(cè)序錯(cuò)誤率是通過(guò)測(cè)序Phrec數(shù)值(Phred score, Qphred)得到，而Phrec數(shù)值是在堿基識(shí)別(Base Calling)過(guò)程通過(guò)一種預(yù)測(cè)堿基判別發(fā)生錯(cuò)誤概率模 型計(jì)算得到的，對(duì)應(yīng)關(guān)系如下表所顯示：Phred分值不正確的堿基識(shí)別堿基正確識(shí)別率101/1090%201/10099%301/100099.9%4

14、01/1000099.99%在Hiseq4000測(cè)序系統(tǒng)測(cè)序時(shí)，首先會(huì)對(duì)文庫(kù)進(jìn)行芯片制備，目的是將文庫(kù)DNA 模板固定到芯片上，在固定DNA模板的過(guò)程中，每個(gè)DNA分子會(huì)形成一個(gè)簇，一個(gè) 簇就是一個(gè)測(cè)序位點(diǎn)，在進(jìn)行固定過(guò)程中極少量的簇與簇之間物理位置會(huì)發(fā)生重疊， 在測(cè)序時(shí)，測(cè)序軟件通過(guò)前4個(gè)堿基對(duì)這些重疊的點(diǎn)進(jìn)行分析和識(shí)別，將這些重疊點(diǎn)位置分開(kāi)，保證每個(gè)點(diǎn)測(cè)到的是一個(gè) DNA分子，因此測(cè)序序列5端前幾個(gè)堿基的錯(cuò) 誤率相對(duì)較高。另外測(cè)序錯(cuò)誤率會(huì)隨著測(cè)序序列(Sequeneed Read)的長(zhǎng)度的增加 而升高，這是由于測(cè)序過(guò)程中化學(xué)試劑的消耗而導(dǎo)致的。因此在進(jìn)行堿基測(cè)序質(zhì)量 分布分析時(shí)，樣品的堿基

15、質(zhì)量分布在前4個(gè)堿基和后十幾個(gè)堿基的質(zhì)量值會(huì)低于中間 測(cè)序堿基，但其質(zhì)量值都高于 Q30,根據(jù)質(zhì)量值和錯(cuò)誤率的關(guān)系，我們將質(zhì)量值轉(zhuǎn) 換成錯(cuò)誤率，繪制錯(cuò)誤率分布圖如下：300100 200Positioni long reads樣品p堿基錯(cuò)誤率分布注：橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿基錯(cuò)誤率，前151bp為雙端測(cè)序序列的第一端測(cè)序reads的錯(cuò)誤率分布情況，后151bp為另一端測(cè)序reads的錯(cuò)誤率分布情況。3.1.3堿基類型分布?jí)A基類型分布用于檢測(cè)有無(wú)AT、GC分離現(xiàn)象，這種分離現(xiàn)象可能是建庫(kù)測(cè)序 過(guò)程中差異擴(kuò)增引起的，直接影響到后續(xù)的分析。高通量測(cè)序的序列為基因組隨機(jī) 打斷后的DN

16、A片段，位點(diǎn)在整個(gè)基因組的分布是近似均勻的，同時(shí)根據(jù)堿基互補(bǔ) 配對(duì)的原則，A與T和C與G的含量分別是一致的。由于測(cè)序儀器本身的局限性， 前幾個(gè)堿基的A/T和C/G含量可能存在著一定波動(dòng)。樣品各堿基比例分布如下所示：pQSltion alongi樣品p各堿基比例分布注：橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為堿基所占的比例；不同顏色代表不同的堿基類型，綠色代表堿基G，藍(lán)色代表堿基C,紅色代表堿基A，紫色代表堿基T,灰色代表測(cè)序中識(shí)別不岀的堿基N。前151bp為雙端測(cè)序序列的第一端測(cè)序reads的堿基分布，后151bp為另一端測(cè)序reads的堿基分布。每個(gè)cycle代表測(cè)序的每個(gè)堿基，如 第一 cy

17、cle即表示該項(xiàng)目所有測(cè)序reads在第一個(gè)堿基的A、T、G、C、N的分布情況。該圖的結(jié)果顯示AT、CG堿基基本不發(fā)生分離，且曲線較平緩，說(shuō)明測(cè)序結(jié)果正常。重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告樣品測(cè)序數(shù)據(jù)評(píng)估統(tǒng)計(jì)113.1.4低質(zhì)量數(shù)據(jù)過(guò)濾測(cè)序得到的原始測(cè)序序列（Sequeneed Reads或者Raw Reads里面含有帶接頭的、低質(zhì)量的Reads,為了保證信息分析質(zhì)量，對(duì)Raw Reads行過(guò)濾，得到CleanReads用于后續(xù)信息分析。數(shù)據(jù)過(guò)濾的主要步驟如下:去除帶接頭（adapte）的reads 若一條reads上N （未能確定出具體的堿基類型）的比例大于 10%，則過(guò)濾掉該 Pair-e nd

18、reads 去除低質(zhì)量reads （質(zhì)量值Q（x的堿基數(shù)占整條read的50%以上）。數(shù)據(jù)過(guò)濾統(tǒng)計(jì)結(jié)果見(jiàn)下表:數(shù)據(jù)過(guò)濾統(tǒng)計(jì)表BMKIDAdapter_Related Inferior percentRaw Reads( %) (%)Clean ReadsB1B2注：BMK ID :百邁客對(duì)項(xiàng)目樣品的統(tǒng)一編號(hào)； Raw_Reads：原始測(cè)序reads數(shù); Adapter_Related :含接頭被過(guò)濾 的reads比例；lnferior_percent: N含量超過(guò)10%的reads和質(zhì)量值低于10的堿基超過(guò)50%的reads比例；Clean_Reads： 過(guò)濾后剩余的reads數(shù)3.1.5測(cè)序

19、數(shù)據(jù)統(tǒng)計(jì)各樣品測(cè)序產(chǎn)出數(shù)據(jù)評(píng)估結(jié)果如下表所示:BMKIDRaw_Reads Clean_ReadsClean_BaseQ30(%)GC(%)PMB1B2注：BMK ID :百邁客對(duì)項(xiàng)目樣品的統(tǒng)一編號(hào)； Raw_Reads：原始測(cè)序reads數(shù)目,以四行為一個(gè)單位，統(tǒng)計(jì)Pair-end 序列的個(gè)數(shù)；Clean_Reads：過(guò)濾后的reads數(shù)，計(jì)算方法同Raw Reads； Clean_Bases：過(guò)濾后的堿基數(shù)，Clean Reads 數(shù)乘以序列長(zhǎng)度；Q30(%):質(zhì)量值大于等于30的堿基占總堿基數(shù)的百分比； GC(%):樣品GC含量，即G和C類 型的堿基占總堿基的百分比。3.2與參考基因組比

20、對(duì)統(tǒng)計(jì)重測(cè)序獲得的測(cè)序reads需要重新定位到參考基因組上，才可以進(jìn)行后續(xù)變異分析。bwaF軟件主要用于二代高通量測(cè)序(如Illunima HiSeq 4000等測(cè)序平臺(tái))得到的短序列與參考基因組的比對(duì)。通過(guò)比對(duì)定位Clea n Reads在參考基因組上的位置，統(tǒng)計(jì)各樣品的測(cè)序深度、基因組覆蓋度等信息，并進(jìn)行變異的檢測(cè)。3.2.1比對(duì)結(jié)果統(tǒng)計(jì)樣品的比對(duì)結(jié)果見(jiàn)下表：比對(duì)結(jié)果統(tǒng)計(jì)BMK IDTotal_readsMapped(%)Properly_mapped(%)PMB1B2注：BMK ID :百邁客對(duì)項(xiàng)目樣品的統(tǒng)一編號(hào)； Total_Reads： Total_Reads數(shù)，雙端分別統(tǒng)計(jì)，即re

21、ad1和read2記 為2條reads； Mapped(%):定位到參考基因組的 Clean Reads數(shù)占所有 Clean Reads數(shù)的百分比；Properly mapped : 雙端測(cè)序序列均定位到參考基因組上且距離符合測(cè)序片段的長(zhǎng)度分布。此項(xiàng)目樣品的平均比對(duì)效率均在 XX%以上，說(shuō)明樣品測(cè)序正常。3.2.2插入片段分布統(tǒng)計(jì)通過(guò)檢測(cè)雙端序列在參考基因組上的起止位置，可以得到樣品DNA打斷后得到的測(cè)序片段的實(shí)際大小，即插入片段大小(Insert Size)，是信息分析時(shí)的一個(gè)重要參數(shù)。插入片段大小的分布一般符合正態(tài)分布，且只有一個(gè)單峰，In sert Size分布圖可以展示各個(gè)樣品的插入片

22、段的長(zhǎng)度分布情況。每個(gè)樣品測(cè)序數(shù)據(jù)插入片段大小分布的分析使用picard軟件工具包中Collectin sertSizeMetric.jar 軟件實(shí)現(xiàn)。DewISVInsert sizeInsert size di血ibuhonLqFnu 會(huì) fllx樣品p插入片段分布圖注：橫坐標(biāo)為插入片段長(zhǎng)度，縱坐標(biāo)為其對(duì)應(yīng)的reads數(shù)。由上圖可知，插入片段長(zhǎng)度分布符合正態(tài)分布，說(shuō)明測(cè)序數(shù)據(jù)文庫(kù)構(gòu)建無(wú)異常。 323深度分布統(tǒng)計(jì)Read定位到參考基因組后，可以統(tǒng)計(jì)參考基因組上堿基的覆蓋情況。參考基因組上被reads覆蓋到的堿基數(shù)占基因組的百分比稱為基因組覆蓋度；堿基上覆蓋的 reads數(shù)為覆蓋深度。基因組覆

23、蓋度可以反映參考基因組上變異檢測(cè)的完整性，覆蓋到的區(qū)域越多， 可以檢測(cè)到的變異位點(diǎn)也越多。覆蓋度主要受測(cè)序深度以及樣品與參考基因組親緣 關(guān)系遠(yuǎn)近的影響。基因組的覆蓋深度會(huì)影響變異檢測(cè)的準(zhǔn)確性，在覆蓋深度較高的區(qū)域（非重復(fù) 序列區(qū)），變異檢測(cè)的準(zhǔn)確性也越高。另外，若基因組上堿基的覆蓋深度分布較均 勻，也說(shuō)明測(cè)序隨機(jī)性較好。樣品的堿基覆蓋深度分布曲線和覆蓋度分布曲線見(jiàn)下圖：重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告D0.00 -102034)4050Sequencing depth6o.“器 qoz-e-nEnooc4e0d0.&樣品P的深度分布注：上圖反映了測(cè)序深度分布的基本情況，橫坐標(biāo)為測(cè)序深度，左縱坐標(biāo)為該深

24、度對(duì)應(yīng)的堿基所占百分比，對(duì) 應(yīng)紅色曲線，右縱坐標(biāo)為該深度及以下的堿基所占百分比，對(duì)應(yīng)藍(lán)色曲線。各樣品的平均覆蓋深度和各深度對(duì)應(yīng)的基因組覆蓋比例如下表所示:樣品覆蓋深度和覆蓋度比例統(tǒng)計(jì)BMK IDAve_depthCov_ratio_1X(%)Cov_ratio_5X(%)Cov_ratio_10X(%)MB1B2注: BMK ID :百邁客對(duì)項(xiàng)目樣品的統(tǒng)一編號(hào)；Ave-depth :樣品平均覆蓋深度；Cov_ratio :覆蓋深度在給定深度及以上的堿基數(shù)占參考基因組總堿基數(shù)的比例。由上表可知，此項(xiàng)目基因組平均覆蓋深度約為 X，基因組覆蓋度約為XX% (至少覆蓋1X)根據(jù)染色體各位點(diǎn)的覆蓋深度情

25、況進(jìn)行作圖，若覆蓋深度在染色體上的分布比 較均勻，則可以認(rèn)為測(cè)序隨機(jī)性比較好。樣品的染色體覆蓋深度分布見(jiàn)下圖：15Geriomewlde diflrlbvtlon of crovenigeJsfinrxp.poikoH 弋-mis r rll-H-MM MiSOMb10ah*Chramcmw ipXiwOrUnCM.OwlCM0*3cwOrtCM0*7CWOriCiria6htO1J樣品p染色體覆蓋深度分布圖注：橫坐標(biāo)為染色體位置，縱坐標(biāo)為染色體上對(duì)應(yīng)位置的覆蓋深度取對(duì)數(shù)(Iog2)得到的值。由上圖可以看出基因組被覆蓋的較均勻，說(shuō)明測(cè)序隨機(jī)性較好。圖上深度不均的地方可能是由于重復(fù)序列、PCR偏

26、好性引起的。3.3 SNP檢測(cè)與注釋3.3.1樣品與參考基因組間SNP的檢測(cè)SNP的檢測(cè)主要使用GATK3軟件工具包實(shí)現(xiàn)。根據(jù)Clean Read在參考基因組的定位結(jié)果，使用Picarcf4進(jìn)行去重復(fù)(Mark Duplicates)、GATK進(jìn)行局部重比對(duì)(Local Realignment)堿基質(zhì)量值校正(Base Recalibration)等預(yù)處理，以保證檢測(cè)得到的 SNP準(zhǔn)確性，再使用GATK進(jìn)行單核苷酸多態(tài)性(Single Nucleotide Polymorphism， SNP)的檢測(cè)，過(guò)濾，并得到最終的SNP位點(diǎn)集。主要檢測(cè)過(guò)程如下： 對(duì)于BWA比對(duì)得到的結(jié)果，使用Picard

27、的Mark Duplicate工具去除重復(fù)，屏 蔽 PCR-duplicatio n 的影響。 使用GATK進(jìn)行In Del Realig nment,即對(duì)存在插入缺失比對(duì)結(jié)果附近的位點(diǎn) 進(jìn)行局部重新比對(duì)，校正由于插入缺失引起的比對(duì)結(jié)果錯(cuò)誤。(3) 使用GATK進(jìn)行堿基質(zhì)量值再校準(zhǔn)(Base Recalibration),對(duì)堿基的質(zhì)量值 進(jìn)行校正。（4）使用GATK進(jìn)行變異檢測(cè)（varia nt calli ng）,主要包括SNP和In Del。對(duì)SNP進(jìn)行嚴(yán)格過(guò)濾：snp clusteM濾（5bp內(nèi)如果有2個(gè) SNP則過(guò)濾掉），In del 附近SNP過(guò)濾（In Del附近5bp內(nèi)的SNP過(guò)

28、濾掉）；和相鄰INDEL過(guò)濾（兩個(gè)In Del距離 小于10bp過(guò)濾掉）5。具體流程可參考GATK官方網(wǎng)站的BestPractice/gatk/guide/best-practices?bpm=DNAseq#varia nt-dis covery-ovw變異結(jié)果使用vcf文件格式展示。vcf文件包括注釋行、標(biāo)題行和數(shù)據(jù)行三部分。 其中注釋行包含文件數(shù)據(jù)行的INFO和FORMAT列中使用的各種標(biāo)識(shí)符的意義解釋， 而標(biāo)題行和數(shù)據(jù)行包含各樣品的變異檢測(cè)結(jié)果信息，格式如下所示：重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告SNP變異結(jié)果信息表#CHR0MPOSIDR

29、JEFALTFIITEFORMATPQUALRINFOChr15634GA140.84 PASSANNOTATIONSGT:AD:DP:GQ:PL1/1:0,6:6:18:169,18,0Chr130071AG141.84 PASSANNOTATIONSGT:AD:DP:GQ:PL1/1:0,6:6:18:170,18,0Chr130478CT95.9PASSANNOTATIONSGT:AD:DP:GQ:PL1/1:0,5:5:15:124,15,0Chr132667AG91.03 PASSANNOTATIONSGT:AD:DP:GQ:PL1/1:0,4:4:12:119,12,0各列意義說(shuō)明

30、如下：1CHROMChr1參考序列的染色體名稱2POS5634參考序列位點(diǎn)坐標(biāo)3ID標(biāo)識(shí)符4REFG參考序列對(duì)應(yīng)位置堿基5ALTASNP位點(diǎn)對(duì)應(yīng)的另外類型的堿基6QUAL140.84變異位點(diǎn)質(zhì)量值7FILTERPASS過(guò)濾狀態(tài)8INFOANNOTATIONS位點(diǎn)注釋信息9FORMATGT:AD:DP:GQ:PL基因型信息格式10R011/1:0,6:6:18:169,18,0樣品的基因型信息VCf文件的詳細(xì)說(shuō)明信息見(jiàn)網(wǎng)頁(yè)：http:/gatkforums.broadi /discussi on/1268/how-should-i-i nterpret-vcf-files

31、-pro duced-by-the-gatk為了確保樣本SNP的可信性，對(duì)樣本檢測(cè)的SNP的reads支持?jǐn)?shù)，相鄰SNP的距 離統(tǒng)計(jì)累積分布。#重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告Cwuisllv* SNP depth discrtkjtian41 MXtkHIenhidoDiplMwf p Fiaa RCrf-Mi 碎g FKKg 3OD 00 4OC:USNuhh Tinqi fiNP dnumcflitinSNP質(zhì)量分布圖注：左邊為 SNP reads支持?jǐn)?shù)目累積圖，右邊為相鄰SNP之間的距離累積圖SNP類型的變異分為轉(zhuǎn)換和顛換兩種，同種類型堿基之間突變稱為轉(zhuǎn)換(Tran sitio n),如嘌呤與

32、嘌呤之間、嘧啶與嘧啶之間的變異，不同類型堿基之間的 突變稱為顛換(Transversion)，如嘌呤與嘧啶之間的變異。一般來(lái)說(shuō)轉(zhuǎn)換比顛換更 容易發(fā)生，故轉(zhuǎn)換/顛換(Ti/Tv )的比例一般大于1,具體數(shù)值和所測(cè)物種有關(guān)。對(duì) 于二倍體或者多倍體物種，若同源染色體上的某一SNP位點(diǎn)均為同一種堿基，貝U該SNP位點(diǎn)稱為純合SNP位點(diǎn)；若同源染色體上的SNP位點(diǎn)包含不同類型的堿基，則 該SNP位點(diǎn)稱為雜合SNP位點(diǎn)。純合SNP數(shù)量越多，則樣品與參考基因組之間差異 越大，雜合SNP數(shù)量越多，則樣品的雜合程度越高，具體結(jié)果和樣品的材料選擇有 關(guān)。332樣品之間SNP的檢測(cè)根據(jù)樣品與參考基因組的比對(duì)結(jié)果，匯

33、總樣品之間所有有差異的變異位點(diǎn)，各 樣品的SNP列表文件格式如下所示：各樣品SNP列表示意#ChrPosRefPMB1B2chromosome_1240CCTCCchromosome_1248GGAGGchromosome_1422AATAAchromosome_1463CNTCCchromosome_1483TNCTTchromosome_1631CCTCCchromosome_1651TTCTT注：Chr: SNP位點(diǎn)所在的染色體名稱；Pos：SNP在參考序列的位置；Ref：參考序列的堿基類型；P、 M、 B1、B2:各樣品在該SNP位點(diǎn)對(duì)應(yīng)的堿基類型。SNP基因型的編碼采用標(biāo)準(zhǔn)核苷酸符號(hào)

34、，符號(hào)表如下所示:核苷酸代碼意義核苷酸代碼意義AAdenosineMA C (aMino group)CCytosineSG C (Strong interaction)GGuanineWA T (Weak interaction)TThymidineBG T C (not A) (B comes after A)UUracilDG A T (not C) (D comes after C)RG A (puRine)HA C T (not G) (H comes after G)YT C (pYrimidine)VG C A (not T, not U) (V comes after U)KG

35、 T (Ketone)NA G C T (aNy)據(jù)統(tǒng)計(jì)，樣品P和M間共有XX個(gè)SNP，樣品B1和B2間共有XX個(gè)SNP。333 SNP結(jié)果注釋Sn pEf嚴(yán)是一款用于注釋變異（SNP、Small In Del）和預(yù)測(cè)變異影響的軟件。根 據(jù)變異位點(diǎn)在參考基因組上的位置以及參考基因組上的基因位置信息，可以得到變 異位點(diǎn)在基因組發(fā)生的區(qū)域（基因間區(qū)、基因區(qū)或CDS區(qū)等），以及變異產(chǎn)生的影響（同義非同義突變等）。軟件可以使用vcf格式文件作為輸入和輸出。輸出結(jié)果會(huì)在 vcf文件的INFO列添 加以下字段：EFF= Effect （ Effect_lmpact | Fun ctio nal_Class

36、 | Codo n_Cha nge | Amino_Acid_Change| Amino_Acid_Length | Gene_Name | Transcript_BioType | Gen e_Codi ng | Tran script_ID | Exon_Ra nk | Ge no type_Number | ERRORS | WARNINGS ）各標(biāo)識(shí)符說(shuō)明如下：類型意義Effect變異所在的區(qū)域或類型Effect impact變異影響大小(High, Moderate, Low, Modifier )Functional Class功能分類(NONE, SILENT, MISSENSE

37、, NONSENSE )Codon_Change/Distance編碼改變（old_codon/new_codon ）或者變異位點(diǎn)到轉(zhuǎn)錄本的距離（在基因 上下游區(qū)域）Amino_Acid_Change氨基酸編碼改變（原氨基酸類型、位置、改變后氨基酸類型）Amino_Acid_Length氨基酸編碼蛋白的長(zhǎng)度（轉(zhuǎn)錄本長(zhǎng)度/3）Gene_Name基因名Transcript_BioType轉(zhuǎn)錄本功能Gene_Coding編碼蛋白(CODING | NON_CODING)Transcript_ID轉(zhuǎn)錄本IDExon/Intron Rank外顯子或內(nèi)含子位次Genotype_Number變異的基因型位次

38、Warnings/Errors警告或錯(cuò)誤以上的結(jié)果若無(wú)法得到，則其對(duì)應(yīng)列為空。具體說(shuō)明可參見(jiàn)SnpEff的說(shuō)明文檔：http:/s npeff.sourceforge. net/S npEff_ma nu al.html#output23本項(xiàng)目樣品間的SNP注釋具體統(tǒng)計(jì)結(jié)果如下所示：SNP注釋結(jié)果統(tǒng)計(jì)TypeP vs MB1 vs B2INTERGENICINTRAGENICINTRONUPSTREAMDOWNSTREAMSPLICE_SITE_ACCEPTORSPLICE_SITE_DONORSTART_LOSTNON_SYNONYMOUS_STARTSYNONYMOUS_CODINGCDS

39、 NON_SYNONYMOUS_CODINGSYNONYMOUS_STOPSTOP_GAINEDSTOP_LOSTOther注：Type : SNP所在區(qū)域或類型；P、M、B1、B2為各樣品相對(duì)于參考基因組存在的對(duì)應(yīng)類型的SNP數(shù)量，P vsM和B1 vs B2為兩個(gè)樣品間存在的對(duì)應(yīng)類型的 SNP數(shù)量。各行意義說(shuō)明如下表所示:INTERGENIC基因間區(qū)INTRAGENIC基因內(nèi)（無(wú)轉(zhuǎn)錄本信息）INTRON內(nèi)含子UPSTREAM基因上游區(qū)域（5K以內(nèi)）DOWNSTREAM基因下游區(qū)域（5K以內(nèi)）SPLICE_SITE_ACCEPTOR剪切供體突變（exon前2bp內(nèi)）SPLICE_SITE_D

40、ONOR剪切受體突變（exon后2bp內(nèi)）NON_SYNONYMOUS_CODING非同義編碼突變NON_SYNONYMOUS_START非同義的起始密碼子突變START_LOST起始密碼子丟失STOP_GAINED終止密碼子獲得STOP_LOST終止密碼子丟失SYNONYMOUS_CODING同義編碼突變SYNONYMOUS_STOP同義終止密碼子突變Other由于gff文件中基因信息不完整而無(wú)法得到準(zhǔn)確的判斷3.4 Small In Del檢測(cè)與注釋3.4.1樣品與參考基因組間Small In Del的檢測(cè)根據(jù)樣品的Clean Read在參考基因組上的定位結(jié)果，檢測(cè)樣品與參考基因組之 間是

41、否存在小片段的插入與缺失（Small In Del: 1-5bp）。樣品的插入缺失使用GATK 檢測(cè)。Small InDel變異一般比SNP變異少，同樣反映了樣品與參考基因組之間的差 異，并且編碼區(qū)的In Del會(huì)引起移碼突變，導(dǎo)致基因功能上的變化。3.4.2樣品之間Small In Del檢測(cè)根據(jù)樣品與參考基因組的Small InDei檢測(cè)結(jié)果，提取樣品之間有差異的變異位 點(diǎn)，即樣品之間的Small InDei變異位點(diǎn)。部分結(jié)果如下表所示：樣品測(cè)序數(shù)據(jù)統(tǒng)計(jì)#CHROMPOSREFALTWMut各列意義說(shuō)明如下:列數(shù)標(biāo)題示例意義說(shuō)明1CHROMLG1參考序列的染色體名稱2POS14890321

42、參考序列位點(diǎn)坐標(biāo)3REFC參考序列對(duì)應(yīng)位置堿基序列4ALTCATInDel位點(diǎn)對(duì)應(yīng)的另外類型的堿基序列56W、Mut0/0、0/1、1/1、./.各樣品對(duì)應(yīng)的InDel類型（0/0 :純合且與參考基因組一 . 致；0/1：雜合類型；1/1：純合且與參考基因組不一 致；./.:不能確定）據(jù)統(tǒng)計(jì)，樣品P和M間共有XX個(gè)In Del，樣品B1和B2間共有XX個(gè)In Del。3.4.3 Small In Del 的注釋根據(jù)樣品檢測(cè)得到的Small In Del位點(diǎn)在參考基因組上的位置信息，對(duì)比參考基 因組的基因、CDS位置等信息（一般在gff文件中），可以注釋In Del位點(diǎn)是否發(fā)生在基 因間區(qū)、基因

43、區(qū)或CDS區(qū)、是否為移碼突變等。Small InDel的注釋通過(guò)SnpEff軟件 實(shí)現(xiàn)。發(fā)生移碼突變的In Del可能會(huì)導(dǎo)致基因功能的改變，具體注釋結(jié)果見(jiàn)下表：重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告In Del注釋結(jié)果統(tǒng)計(jì)TypeP vs MB1vs B2INTERGENICINTRAGENICINTRONUPSTREAMDOWNSTREAMSPLICE_SITE_ACCEPTORSPLICE_SITE_DONORSTART_LOSTFRAME_SHIFTCODON_DELETIONCODON_INSERTIONCDSEXON_DELETIONCODON_CHANGE_PLUS_CODON_DELETION

44、CODON_CHANGE_PLUS_CODON_INSERTIONSTOP_GAINEDSTOP_LOSTOther注：Type: InDel所在區(qū)域或類型；P、M、B1、B2為各樣品相對(duì)于參考基因組存在的對(duì)應(yīng)類型的InDeI數(shù)量，P vsM和B1 vs B2為兩個(gè)樣品間存在的對(duì)應(yīng)類型的 InDel數(shù)量。#重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告INTERGENICINTRAGENICINTRONUPSTREAMDOWNSTREAMSPLICE_SITE_ACCEPTORSPLICE_SITE_DONORCODON_CHANGE_PLUS_CODON_DELE TIONCODON_CHANGE_PLUS_CO

45、DON_INSERTIONCODON_DELETIONCODON_INSERTIONEXON_DELETEDFRAME_SHIFTSTART_LOSTSTOP_GAINEDSTOP_LOSTOther基因間區(qū)基因內(nèi)（無(wú)轉(zhuǎn)錄本信息）內(nèi)含子基因上游區(qū)域（5K以內(nèi)）基因下游區(qū)域（5K以內(nèi)）剪切供體突變（exon前2bp內(nèi)）剪切受體突變（exon后2bp內(nèi)）非密碼子邊界上的3的整數(shù)倍的刪除非密碼子邊界上的3的整數(shù)倍的插入密碼子刪除（3的整數(shù)倍）密碼子插入（3的整數(shù)倍）整個(gè)外顯子被刪除移碼突變（非3的整數(shù)倍插入或刪除）起始密碼子丟失終止密碼子獲得終止密碼子丟失由于gff文件中基因信息不完整而無(wú)法得到準(zhǔn)確

46、的判斷253.5關(guān)聯(lián)分析3.5.1高質(zhì)量SNP篩選根據(jù)SNP檢測(cè)結(jié)果，樣品P和M共篩選到XX個(gè)SNP位點(diǎn)，樣品B1和B2 共XX個(gè)SNP，在關(guān)聯(lián)分析前，首先對(duì)樣品 B1和B2間的XX個(gè)SNP進(jìn)行過(guò)濾， 過(guò)濾標(biāo)準(zhǔn)如下：首先過(guò)濾掉有多個(gè)基因型的SNP位點(diǎn)，其次過(guò)濾掉read支持度小于4的SNP位點(diǎn)，再次根據(jù)通過(guò)親本的SNP信息過(guò)濾掉與同表型親本不同的位點(diǎn)， 最終得到高質(zhì)量的可信SNP位點(diǎn)XX個(gè)。SNP過(guò)濾統(tǒng)計(jì)Total SNP多個(gè)等位基因的位點(diǎn)Read支持度小于4的位點(diǎn)利用親本過(guò)濾的位點(diǎn)總過(guò)濾位點(diǎn)數(shù)高質(zhì)量SNPB1 vs B2107042155032254280053.5.2 SNP-index方

47、法關(guān)聯(lián)結(jié)果SNP-index是近年來(lái)發(fā)表的一種通過(guò)混池間的基因型頻率差異進(jìn)行標(biāo)記關(guān)聯(lián)分 析的方法，主要是尋找混池之間基因型頻率的顯著差異，用 (SNPindex)統(tǒng)計(jì)。標(biāo)記SNP與性狀關(guān)聯(lián)度越強(qiáng)， (SNPndex)越接近于1。計(jì)算方法簡(jiǎn)述如下：Maa表示aa池來(lái)源于母本的深度；Paa表示aa池來(lái)源于父本的深度；Mab表示ab池來(lái)源于母本的深度；Pab表示ab池來(lái)源于父本的深度；SNP-index(ab)=Mab/(Pab+Mab); SNP-index(aa)=Maa/(Paa+Maa) (SNP ndex)=SNP-i ndex(aa)-SNP-i ndex(ab。為了消除假陽(yáng)性的位點(diǎn)，禾

48、U用標(biāo)記在基因組上的位置，可對(duì)同一條染色體上標(biāo) 記的 SNPindex值進(jìn)行擬合，本項(xiàng)目并采用DISTANCE方法對(duì) SNP-index進(jìn)行 擬合，取每個(gè)SNP左右距離各2M的SNP的厶SNP-index的中值作為該位點(diǎn)擬合后 的關(guān)聯(lián)值。并根據(jù)關(guān)聯(lián)閾值，選擇閾值以上的區(qū)域作為與性狀相關(guān)的區(qū)域。兩個(gè)混池分別的SNP-index及厶SNP-index的分布如下圖所示：重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告SNP-index關(guān)聯(lián)值在染色體上的分布注：橫坐標(biāo)為染色體名稱，彩色的點(diǎn)代表計(jì)算出來(lái)的SNP-index (或 SNP-index )值，黑色的線為擬合后的SNP-index (或 SNP-index )值。上

49、圖是 B1混池的SNP-index值的分布圖；中圖是 B2混池的SNP-index值的分布圖；下圖是 SNP-index值的分布圖，其中紅色的線代表置信度為0.99的閾值線，藍(lán)色的線代表置信度為0.95的閾值線，綠色的線代表置信度為 0.90的閾值線。根據(jù)本項(xiàng)目群體的理論分離比，計(jì)算關(guān)聯(lián)閾值為XX。根據(jù)計(jì)算機(jī)模擬實(shí)驗(yàn)計(jì)算結(jié)果，當(dāng)置信度為0.95時(shí)，定位區(qū)域?yàn)?158439337-160141497( 1.702M)區(qū)間內(nèi)，共得到 XX個(gè)區(qū)域，總長(zhǎng)度為 XXbp，其 中包含非同義突變SNP位點(diǎn)的基因共XX個(gè)，同義突變SNP位點(diǎn)的基因共XX個(gè)。理論上，目標(biāo)位點(diǎn)及其附近的連鎖位點(diǎn)應(yīng)趨近于該閾值，因此

50、顯著關(guān)聯(lián)的區(qū)域附近應(yīng)該出現(xiàn)一個(gè)較高的峰值。但從結(jié)果上看，沒(méi)有超過(guò)理論閾值的區(qū)域，說(shuō)明本 實(shí)驗(yàn)中沒(méi)有發(fā)現(xiàn)顯著的定位結(jié)果。為了充分利用數(shù)據(jù)，將閾值降低以尋找比較可能 的定位區(qū)域，利用擬合后 SNP-index的99百分位數(shù)，即XXX，最終得到可能的 定位區(qū)域?yàn)?58439337-160141497( 1.702M)，共得到XX個(gè)區(qū)域，總長(zhǎng)度為XXbp， 其中包含非同義突變SNP位點(diǎn)的基因共XX個(gè)，同義突變SNP位點(diǎn)的基因共XX 個(gè)，移碼突變的基因共 XX個(gè)。然而由于未達(dá)到理論閾值，這個(gè)區(qū)域很可能是假陽(yáng) 性區(qū)域，需要進(jìn)一步驗(yàn)證。關(guān)聯(lián)區(qū)域信息統(tǒng)計(jì)表Chromosome IDStartEndSize(M

51、b)Gene nu mberChr0211,688,75512,993,3191.30148Total-1.30148注：Chromosome ID :染色體編號(hào)；Start：關(guān)聯(lián)區(qū)域起始位置；End:關(guān)聯(lián)區(qū)域終止位置；Size:關(guān)聯(lián)區(qū)域大小, 以Mb為單位；Gene number：關(guān)聯(lián)區(qū)域內(nèi)的基因數(shù)量。2重測(cè)序BSA項(xiàng)目結(jié)題報(bào)告3.5.3 ED萬(wàn)法關(guān)聯(lián)結(jié)果歐式距離（Euclidean Distanee, ED）算法，是利用測(cè)序數(shù)據(jù)尋找混池間存在顯 著差異標(biāo)記，并以此評(píng)估與性狀關(guān)聯(lián)區(qū)域的方法 9。理論上，BSA項(xiàng)目構(gòu)建的兩個(gè) 混池間除了目標(biāo)性狀相關(guān)位點(diǎn)存在差異，其他位點(diǎn)均趨向于一致，因此非目標(biāo)

52、位點(diǎn) 的ED值應(yīng)趨向于0。ED方法的計(jì)算公式如下所示，ED值越大表明該標(biāo)記在兩混 池間的差異越大。JlAmz A財(cái)尸 十（匚加+ （Gn世Gm），+ 血 一 T艸/其中：35Amut為A堿基在突變混池中的頻率,Cmut為C堿基在突變混池中的頻率,Gmut為G堿基在突變混池中的頻率,Tmut為T堿基在突變混池中的頻率，Awt為A堿基在野生型混池中的頻率;Cwt為C堿基在野生型混池中的頻率；Gwt為G堿基在野生型混池中的頻率;Twt為T堿基在野生型混池中的頻率。在進(jìn)行分析時(shí)，利用兩混池間基因型存在差異的 SNP位點(diǎn)，統(tǒng)計(jì)各個(gè)堿基在不 同混池中的深度，并計(jì)算每個(gè)位點(diǎn) ED值，為消除背景噪音，對(duì)原始 ED值進(jìn)行乘 方處理9，本項(xiàng)目取原始ED的4次方作為關(guān)聯(lián)值以達(dá)到消除背景噪音的功能，然 后采用局部線性回歸LOESS方法對(duì)ED值進(jìn)行擬合，關(guān)聯(lián)值分布如下圖所示：3EH9UED關(guān)聯(lián)值在染色體上的分布注：橫坐標(biāo)為染色體名稱，彩色的點(diǎn)代表每個(gè)SNP位點(diǎn)的ED值，黑色的線為擬合后的 ED值，紅色的虛線代表顯著性關(guān)聯(lián)閾值，ED值越高，代表該點(diǎn)關(guān)聯(lián)效果越好。取所有位點(diǎn)擬合值的median+3SD作為分析的關(guān)聯(lián)閾值9，計(jì)算得XX。根據(jù)關(guān)聯(lián)閾值判定，定位區(qū)域?yàn)?158439337-160141497( 1.702