2重測序BSA分析項目結(jié)題報告

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

2、.2.3 深度分布統(tǒng)計 123.3 SNP 檢測與注釋 14331樣品與參考基因組間SNP的檢測 14332樣品之間SNP的檢測 173.3.3 SNP結(jié)果注釋193.4 Small In Del 檢測與注釋223.4.1 樣品與參考基因組間Small InDel 的檢測 223.4.2樣品之間 Small InDel檢測 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)計 363.7結(jié)果可視化374 數(shù)據(jù)下載 384.1 結(jié)果文件查看說明 38參考文獻(xiàn) 391重測序BSA項目結(jié)題報告1項目概況1.1合同關(guān)鍵指標(biāo)(1) 完成X個樣品的重測序，共產(chǎn)生XGbp Clean Data,保證Q30達(dá)到80%。(2) 數(shù)據(jù)評估：測序數(shù)據(jù)量，測序數(shù)據(jù)質(zhì)量和 GC含量的統(tǒng)計。(3) 與基因組比對：比對效率，基因組覆蓋度，基因組覆蓋深度統(tǒng)計。變異檢測和注釋：SNP、In Del的檢測和注釋。(5)

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

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

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

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

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

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

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

11、r詳細(xì)信息見如下：Illumina測序標(biāo)識詳細(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通過使用第四行中每個字符對應(yīng)的 ASCII值進(jìn)行計算，即得到對應(yīng)第二行堿基 的測序質(zhì)量值。如果測序錯誤率用e表示，lllunima HiSeq 4000的堿基質(zhì)量值用Qphred 表示，則有下列關(guān)系：Qphred = -10lOg10(e)Mun ima Casava 1.8版本測序錯誤率與測序質(zhì)量值簡明對應(yīng)關(guān)系如下表所

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

28、濾掉）；和相鄰INDEL過濾（兩個In Del距離 小于10bp過濾掉）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)識符的意義解釋， 而標(biāo)題行和數(shù)據(jù)行包含各樣品的變異檢測結(jié)果信息，格式如下所示：重測序BSA項目結(jié)題報告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各列意義說明

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

31、-pro duced-by-the-gatk為了確保樣本SNP的可信性，對樣本檢測的SNP的reads支持?jǐn)?shù)，相鄰SNP的距 離統(tǒng)計累積分布。#重測序BSA項目結(jié)題報告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)，如嘌呤與嘧啶之間的變異。一般來說轉(zhuǎn)換比顛換更 容易發(fā)生，故轉(zhuǎn)換/顛換(Ti/Tv )的比例一般大于1,具體數(shù)值和所測物種有關(guān)。對 于二倍體或者多倍體物種，若同源染色體上的某一SNP位點均為同一種堿基，貝U該SNP位點稱為純合SNP位點；若同源染色體上的SNP位點包含不同類型的堿基，則 該SNP位點稱為雜合SNP位點。純合SNP數(shù)量越多，則樣品與參考基因組之間差異 越大，雜合SNP數(shù)量越多，則樣品的雜合程度越高，具體結(jié)果和樣品的材料選擇有 關(guān)。332樣品之間SNP的檢測根據(jù)樣品與參考基因組的比對結(jié)果，匯

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

34、，符號表如下所示:核苷酸代碼意義核苷酸代碼意義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)計，樣品P和M間共有XX個SNP，樣品B1和B2間共有XX個SNP。333 SNP結(jié)果注釋Sn pEf嚴(yán)是一款用于注釋變異（SNP、Small In Del）和預(yù)測變異影響的軟件。根 據(jù)變異位點在參考基因組上的位置以及參考基因組上的基因位置信息，可以得到變 異位點在基因組發(fā)生的區(qū)域（基因間區(qū)、基因區(qū)或CDS區(qū)等），以及變異產(chǎn)生的影響（同義非同義突變等）。軟件可以使用vcf格式文件作為輸入和輸出。輸出結(jié)果會在 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)識符說明如下：類型意義Effect變異所在的區(qū)域或類型Effect impact變異影響大小(High, Moderate, Low, Modifier )Functional Class功能分類(NONE, SILENT, MISSENSE

37、, NONSENSE )Codon_Change/Distance編碼改變（old_codon/new_codon ）或者變異位點到轉(zhuǎn)錄本的距離（在基因 上下游區(qū)域）Amino_Acid_Change氨基酸編碼改變（原氨基酸類型、位置、改變后氨基酸類型）Amino_Acid_Length氨基酸編碼蛋白的長度（轉(zhuǎn)錄本長度/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警告或錯誤以上的結(jié)果若無法得到，則其對應(yīng)列為空。具體說明可參見SnpEff的說明文檔：http:/s npeff.sourceforge. net/S npEff_ma nu al.html#output23本項目樣品間的SNP注釋具體統(tǒng)計結(jié)果如下所示：SNP注釋結(jié)果統(tǒng)計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為各樣品相對于參考基因組存在的對應(yīng)類型的SNP數(shù)量，P vsM和B1 vs B2為兩個樣品間存在的對應(yīng)類型的 SNP數(shù)量。各行意義說明如下表所示:INTERGENIC基因間區(qū)INTRAGENIC基因內(nèi)（無轉(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文件中基因信息不完整而無法得到準(zhǔn)確的判斷3.4 Small In Del檢測與注釋3.4.1樣品與參考基因組間Small In Del的檢測根據(jù)樣品的Clean Read在參考基因組上的定位結(jié)果，檢測樣品與參考基因組之 間是

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

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

43、區(qū)或CDS區(qū)、是否為移碼突變等。Small InDel的注釋通過SnpEff軟件 實現(xiàn)。發(fā)生移碼突變的In Del可能會導(dǎo)致基因功能的改變，具體注釋結(jié)果見下表：重測序BSA項目結(jié)題報告In Del注釋結(jié)果統(tǒng)計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為各樣品相對于參考基因組存在的對應(yīng)類型的InDeI數(shù)量，P vsM和B1 vs B2為兩個樣品間存在的對應(yīng)類型的 InDel數(shù)量。#重測序BSA項目結(jié)題報告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)（無轉(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ù)倍）整個外顯子被刪除移碼突變（非3的整數(shù)倍插入或刪除）起始密碼子丟失終止密碼子獲得終止密碼子丟失由于gff文件中基因信息不完整而無法得到準(zhǔn)確

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

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

48、U用標(biāo)記在基因組上的位置，可對同一條染色體上標(biāo) 記的 SNPindex值進(jìn)行擬合，本項目并采用DISTANCE方法對 SNP-index進(jìn)行 擬合，取每個SNP左右距離各2M的SNP的厶SNP-index的中值作為該位點擬合后 的關(guān)聯(lián)值。并根據(jù)關(guān)聯(lián)閾值，選擇閾值以上的區(qū)域作為與性狀相關(guān)的區(qū)域。兩個混池分別的SNP-index及厶SNP-index的分布如下圖所示：重測序BSA項目結(jié)題報告SNP-index關(guān)聯(lián)值在染色體上的分布注：橫坐標(biāo)為染色體名稱，彩色的點代表計算出來的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ù)本項目群體的理論分離比，計算關(guān)聯(lián)閾值為XX。根據(jù)計算機(jī)模擬實驗計算結(jié)果，當(dāng)置信度為0.95時，定位區(qū)域為 158439337-160141497( 1.702M)區(qū)間內(nèi)，共得到 XX個區(qū)域，總長度為 XXbp，其 中包含非同義突變SNP位點的基因共XX個，同義突變SNP位點的基因共XX個。理論上，目標(biāo)位點及其附近的連鎖位點應(yīng)趨近于該閾值，因此

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

51、b)Gene nu mberChr0211,688,75512,993,3191.30148Total-1.30148注：Chromosome ID :染色體編號；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重測序BSA項目結(jié)題報告3.5.3 ED萬法關(guān)聯(lián)結(jié)果歐式距離（Euclidean Distanee, ED）算法，是利用測序數(shù)據(jù)尋找混池間存在顯 著差異標(biāo)記，并以此評估與性狀關(guān)聯(lián)區(qū)域的方法 9。理論上，BSA項目構(gòu)建的兩個 混池間除了目標(biāo)性狀相關(guān)位點存在差異，其他位點均趨向于一致，因此非目標(biāo)

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