不同碳源培養下解脂耶氏酵母的轉錄組差異分析

常 晨，劉偉豐，郭 陽，包怡紅,*

（1.東北林業大學林學院，黑龍江 哈爾濱 150040；2.中國科學院微生物所，北京 100020）

解脂耶氏酵母（Yarrowia lipolytica）是一種嚴格好氧菌，通常生活在含脂類和蛋白質較多的地方，生長于25～30 ℃、pH 3～8的環境下。Y. lipolytica已經有了很長的工業應用史，它的生理生化特性和基因特性使其成為代謝工程最好的宿主[1-5]。

Y. lipolytica的生理活性和分子生物學特點有很多：1）生長速度快、產量高、安全無致病性，因此可以安全地應用于食品和藥物生產，有研究表明Y. lipolytica在干香腸的發酵過程中生長良好，發揮了重要的發酵作用。同時，該酵母能以凝乳為底物發酵產生高品質的奶酪[6-8]；2）對蛋白質的糖基化修飾程度低，與高等真核生物相似；3）具有非常強的表達異源蛋白的能力，因此也大量的應用在有機酸的生產中[9-10]；4）促進生物燃料的合成[11]；5）可以對因油脂引起的水污染進行生物修復[12]。

高通量測序技術也被稱為下一代測序技術[13-14]，該技術挖掘研究對象基因的轉錄信息迅速，并且十分全面，在生物體轉錄組基因表達分析中普遍采用該技術，精確便捷地挖掘出相關功能基因[15-17]。Y. lipolytica由于其獨特的生理活性和分子生物學特點，近二十年來獲得了研究者越來越多的關注。本研究利用高通量測序技術，探索不同碳源對解脂酵母生長的影響。

1 材料與方法

1.1 材料與試劑

1.1.1 菌株

Y. lipolytica菌株Po1f為實驗室所保存。

1.1.2 培養基

YPD、YNB等培養基均按照文獻[18]配方配制。

RiboMinus Transcriptome Isolation Kit、Ribominus Kit、PureLink PCR Micro Kit、Novex precast gel products 美國英杰（Invitrogen）生命技術有限公司；total RNA-Seq Kit 美國應用生物系統公司；TRlzol Reagent 美國Life Technologies公司；反轉錄試劑盒HiScript 1st Strand cDNA Synthesis Kit 普洛麥格（北京）生物技術有限公司；rTaqDNA聚合酶、DNase I、DNA Marker 寶日醫生物技術（北京）有限公司；所有引物合成于蘇州金唯智有限公司；其他常規生化試劑均為國產分析純。

1.2 儀器與設備

恒溫搖床 上海智誠公司；離心機、低溫離心機德國Sigma公司；瓊脂糖核酸電泳儀 北京百晶生物公司；全自動凝膠成像儀 上海天能科技有限公司；PB-10 pH計、分析天平 上海精密儀器儀表有限公司。

1.3 方法

1.3.1 樣品制備

將Y. lipolytica菌株Po1f劃線于YPD固體平板，24 h后挑取生長良好的單菌落，接種至YPD液體培養基，30 ℃、200 r/min過夜培養，以1%的接種量分別轉接到碳源為葡萄糖和油酸的YNB培養基，30 ℃、200 r/min培養至對數期。4 ℃離心，沉降菌體，棄上清液，加RNase free的磷酸鹽緩沖液重懸，離心，棄上清液，轉移至凍存管中，立即放進液氮中速凍，然后于液氮中保存，后續用于轉錄組RNA樣品抽提，同時制備2 個平行重復樣品。

1.3.2 RNA提取與文庫構建

使用TRIzol Reagent（Invitrogen）法[19]提取每個樣品的總RNA。通過Agilent 2100 Bioanalyzer和1%瓊脂糖凝膠鑒定樣品RNA質量，通過NanoDrop對RNA定量。取1 μg RNA完整值大于7的總RNA用于文庫的制備。文庫的制備參照NEBNext?UltraTMRNA文庫制備試劑盒操作說明，包括以下步驟：首先使用NEBNext Poly（A）mRNA Magnetic Isolation Module（NEB）進行poly（A）mRNA分離；然后使用NEB Next First Strand Synthesis Reaction Buffer和NEB Next Random Primers使mRNA片段化，以片段化的mRNA為模板，使用ProtoScript II逆轉錄酶合成第一鏈cDNA，并使用Second Strand Synthesis Enzyme Mix合成第二鏈cDNA；之后用AxyPrep Mag PCR Clean-up純化合成的二鏈cDNA；最后用End Prep Enzyme Mix對純化的cDNA進行末端修復、添加A尾，然后進行TA連接以向兩端添加測序接頭；最后使用AxyPrep Mag PCR Cleanup（Axygen）純化PCR產物以得到最終文庫，之后使用Agilent 2100 Bioanalyzer驗證，并通過Qubit定量，使用Q-PCR方法對所構建文庫的有效濃度進行準確定量，之后通過GENEWIZ對序列進行處理和分析。

1.3.3 測序數據過濾

測序期間，接頭序列會測到少量reads，或者測序長度太長，引起reads的3’端bases質量偏低，會對進一步的數據處理與分析產生影響。所以，對初始數據預處理，以過濾掉質量較低的數據，消除污染和接頭序列的影響[20-21]。測序數據過濾所使用的軟件為Trimmomatic（v0.30）。

1.4 測序數據分析

本研究測序數據的處理分析利用高通量Illumina HiSeq測序平臺，在建立了不同碳源培養出的Y. lipolytica的轉錄組數據庫后，可以對數據進行全面的分析與注釋，以葡萄糖為碳源生長的Y. lipolytica為對照組，以油酸為碳源生長的Y. lipolytica為實驗組。本次研究的具體流程主要包括：原始序列質量評估、參考基因比對分析、基因表達水平分析。

2 結果與分析

2.1 測序數據質量評估與分析

影響測序堿基質量有如下幾種原因：測序所用試劑、測序機器以及樣品自身等。一般情況下發生錯誤率比較高的是5’端的前幾個堿基，但3’端堿基錯誤率會隨著測序序列長度的延長逐漸增加，這種現象是高通量測序技術的特點所導致的。一般情況下，每個堿基位置的測序錯誤率應該低于0.50%[22-23]。測序錯誤率用e表示，Illumina HiSeqTM的堿基質量值用Qphred表示，兩者關系如下：

兩個樣品測序原始數據如表1所示，測序數據過濾后如表2所示。

表1 測序初始數據的質量統計Table 1 Quality statistics of original sequencing data

注：Q20、Q30.計算Phred數值時大于20、30的堿基占總堿基的百分比；GC. G、C的堿基數量總和占總堿基數量的百分比；N.每百萬堿基中N的數量。表2同。

表2 過濾后數據質量統計Table 2 Quality statistics of fi ltered data

從表1、2可以看出，對照組的轉錄組包含18 223 407 raw reads，過濾后得到17 923 921 clean reads，Q20值97.46%，Q30值83.94%，序列GC含量51.69%；實驗組的轉錄組共得到22 680 303 raw reads，經過過濾后得到22 656 852 clean reads，Q20值99.62%，Q30值96.70%，序列GC含量55.42%。

本研究用FastQC（v0.10.1）分析測序質量，對照組和實驗組的堿基位置質量分數結果如圖1所示，樣品堿基序列的平均質量分布情況如圖2所示。

圖1 對照組（A）和實驗組（B）clean reads堿基位置質量分數圖Fig. 1 Clean reads base position quality score maps of control group (A)and experimental group (B)

測序質量分數越高，堿基可信程度越高，通常堿基的質量分數為13%時，錯誤率為5%，質量分數為20%時，錯誤率為1%，質量分數為30%時，錯誤率為0.10%。根據圖1顯示對照組和實驗組的堿基位置質量分數均很高，大部分都在30%以上，錯誤率不足0.10%。

圖2 對照組（A）和實驗組（B）clean reads堿基序列平均質量分布圖Fig. 2 Clean reads mean sequence quality maps of control group (A)and experimental group (B)

當絕大部分堿基序列的平均質量分數的峰值大于30%時，序列質量較好。根據圖2顯示對照組的堿基平均質量分數峰值在36%，實驗組的堿基平均質量分數峰值在37%。由此可見，此次測序數據質量符合標準，可進行下一步分析。

2.2 參考序列比對分析

用已經獲得完整注釋的Y. lipolytica作為參考序列，將過濾后的測序clean data與參考基因組進行對比分析，選擇合適的參考基因組對信息分析成功至關重要，數據比對率可以在一定程度上反映實驗測序樣品與選擇的參考基因組的相似性關系。使用Hisat2（v2.0.1）軟件[24]，默認參數進行短reads的比對。

2.2.1 clean data與參考基因組比對分析

表3 clean reads與參考基因組的比對結果Table 3 Comparison between clean reads and reference genomes

從表3可以看到，2 個樣品定位到基因組上的測序序列數均大于70%，而多重比對測序序列數均小于10%，說明本次研究的測序序列沒有受到污染。

2.2.2 reads在參考基因組不同區域的分布情況

通過對定位到基因組上的測序序列數比對到基因組上的情況統計分析，將區域定位為3 個部分，分別是：外顯子、內含子和基因間隔區域。當測序序列定位到內含子時，通常是有非成熟的mRNA污染，或者注釋不完全的基因組。當測序序列定位到基因間隔區域時，通常是背景噪音，或者注釋不完全的基因組。

圖3分別為對照組和實驗組的reads在參考基因組不同區域的分布情況，可以看出兩個樣品的reads定位到外顯子的區域最多，分別為78.65%和86.90%；定位到內含子區域分別為20.21%和12.10%；定位到基因間隔區域最少，分別為1.13%和1%。

圖3 對照組（A）和實驗組（B）的reads在參考基因組不同區域的分布情況Fig. 3 Distribution of reads in control group (A) and experimental group (B) in different regions of reference genome

2.2.3 reads在染色體上的密度分布情況

通過對定位到基因組上測序序列數比對到基因組上的各個染色體進行統計分析。計算窗口內部比對到堿基位置上的reads數目，并計算其在染色體上的深度分布，并取log2值。一般情況下，定位到該染色體內部的reads總數，會隨著染色體長度的增加而增多。從定位到染色體上的reads數與染色體長度的關系圖中，可以直觀看出測序的均勻度。圖4為對照組和實驗組的reads在染色體上的密度分布圖，縱坐標表示序列在染色體上的深度分布并取log2值，橫坐標表示染色體的長度。

圖4 對照組（A）和實驗組（B）的reads在染色體上的密度分布圖Fig. 4 Density distribution of reads on chromosome in control group (A) and experimental group (B)

2.3 基因表達水平分析

衡量基因的表達水平由其在基因上的豐度判斷，基因豐度越高，基因表達水平越高。本次研究通過使用HtSeq軟件（V 0.6.1）計算基因表達，該軟件使用每百萬reads中來自于某基因每千堿基長度的reads數（reads per kilobase per million mapped reads，RPKM）計算基因表達量[25]。其公式為：

在公式中比對到基因組外顯子上的reads數與比對到全基因組上的reads數的比值表示所有reads數中定位到這個基因的百分比，然后除以外顯子長度。表4分別統計兩個樣品的不同表達水平區間中的基因數。

表4 不同表達水平區間的基因數量統計Table 4 Gene number statistics for different expression level intervals

2.3.1 RNA-seq整體質量評估

2.3.1.1 飽和度檢查

為評估在當前測序深度下的RPKM，通過軟件RSeQC對總的比對reads重采樣，同時采用相對錯誤率衡量評估的RPKM的準確性。

式中：RPKMobs表示每個百分比抽樣下的當前轉錄本的RPKM值；RPKMreal表示總數據量下當前轉錄本的RPKM值。

定量飽和曲線檢查是為了明確基因表達水平定量與數據量之間的關系。當基因的表達量越高時，越容易被準確定量；相反，當基因的表達量低時，越難被準確定量，也就是需要更大的測序數據量才能夠被準確定量。圖5、6為對照組和實驗組的定量飽和曲線檢查分布圖，橫坐標表示定位到基因組上的reads數占總reads數的百分比，縱坐標則是相對誤差百分比。

圖5 對照組的定量飽和曲線檢查分布圖Fig. 5 Quantitative saturation curves of control group

圖6 實驗組的定量飽和曲線檢查分布圖Fig. 6 Quantitative saturation curves of experimental group

從圖5、6可以看出，隨著重采樣比例的增加，相對誤差百分比不斷減少，也就說明定位到基因組上的reads數占總reads數的百分比增加，同時從圖5、6還可以看出轉錄本表達水平越高相對誤差越小。

2.3.1.2 均一性檢查

理想情況下，對于轉錄組測序，測序reads之間為獨立抽樣，所有表達的轉錄本上的reads呈現均一化分布。然而，許多研究表明，很多因素都會影響這種均一化的分布。例如，建庫過程中，片段斷裂和RNA逆轉錄的順序不同，導致RNA-seq最終的數據呈現嚴重的3’偏好性；并且生物體內從5’或者3’的降解過程同樣會導致分布的不均一。

均一性分布的曲線算法是：1）把每個轉錄本從5’到3’的方向平均等分為100 份；2）計算每個等份中的堿基平均測序深度，并用最大值來進行歸一化處理。

圖7表示對照組和實驗組不同長度的轉錄本的reads密度分布圖，不同的顏色表示不同長度范圍的轉錄本。2.3.2 差異表達分析

圖7 對照組（A）和實驗組（B）的均一性分布曲線Fig. 7 Curves of homogeneity distribution in control group (A) and experimental group (B)

不同的樣品，基因差異分析所用軟件也不同。對于有生物學重復的樣品，其基因差異分析應使用Bioconductor軟件包的DESeq2（V1.6.3）進行分析[26-27]，該分析方法所基于的模型是負二項分布，第i個基因在第j個樣本中的read count值為Kij， 則：

對于沒有生物學重復的樣品，其基因差異分析需用Bioconductor軟件包的DESeq（v1.18.0）進行分析，該分析方法基于的模型同樣是負二項分布，第i個基因在第j個樣本中的read count值為Kij，則有：

在一些特殊情況下，需用Bioconductor軟件包的EdgeR（V3.4.6）進行基因差異分析。本次分析所使用的軟件是EdgeR[28]。

2.3.2.1 基因表達水平對比

圖8、9分別為兩個樣品的RPKM的分布圖和盒形圖。在盒型圖中有5 個統計量，分別是最大值、上四分位數、中值、下四分位數和最小值。從圖8、9可以很直觀地看出利用不同碳源生長的Y. lipolytica之間基因表達的差異。

從圖8可以看出，當lgRPKM值小于1.5時對照組的基因密度較小；而當lgRPKM值大于1.5時，對照組的基因密度大于實驗組。從圖9可以看出，不同碳源條件下實驗組和對照組的數據分散程度基本相同，且對照組的中值較大，因此對照組的整體表達量較高。

圖8 RPKM分布圖Fig. 8 RPKM distribution map

圖9 RPKM盒形圖Fig. 9 RPKM box chart

2.3.2.2 差異表達基因篩選

為探究利用不同碳源生長的Y. lipolytica之間差異表達基因，對樣品數據進行差異表達分析。差異表達基因篩選條件為：差異基因表達變化2 倍以上且FDR≤0.05（P值：統計差異顯著性值；FDR：FDR校正P值）。結果顯示，實驗組和對照組之間有536 個顯著性差異表達基因，其中376 個基因上調表達，160 個基因下調表達。圖10為火山圖，表示基因差異表達顯著性變化的分布情況，橫坐標表示基因在不同樣本中表達倍數變化，縱坐標表示基因表達量變化差異的統計學顯著性。

圖10 差異基因火山圖Fig. 10 Volcano map of differentially expressed genes

由圖10可以看出，很多上調基因的-lgFDR很大，說明這些上調基因的表達量變化差異很顯著。

2.3.2.3 差異基因GO富集分析

對兩個樣品的數據進行整理和分析，把篩選得到的差異基因進行GO富集分析，可以明確實驗樣品的差異基因在哪些基因功能上得到體現。GO富集分析涵蓋了3 個方面，分別描述基因的分子功能、細胞組分、參與的生物過程，因此，GO功能顯著性富集分析能給出差異表達基因與哪些生物學功能顯著相關。

本研究通過對差異表達基因進行GO分析，得到了差異基因GO富集柱狀圖（圖11）。在GO富集柱狀圖中，縱坐標為富集的GO term，橫坐標為該term中差異基因的個數，同時為了區分生物過程、細胞組分和分子功能3 個模塊，采用不同顏色分別表示。本研究挑選了富集最顯著的25 個GO term在圖中展示。

圖11 差異基因GO富集柱狀圖Fig. 11 GO enrichment histogram of differentially expressed genes

如圖11所示，在分子功能分組中涉及催化活性、結合、轉運活性的差異基因比較多；在細胞組分的功能分組中涉及細胞膜組分和細胞組分的差異基因比較多；在生物學過程功能分組中，涉及代謝過程的差異表達基因最多其次為細胞過程和單組織過程。

2.3.2.4 差異基因KEGG富集分析

在生物體內，不同基因相互協調發揮其生物學功能，通過通路顯著性富集分析，以KEGG通路為單位，應用超幾何檢驗，能發現差異表達基因參與的最主要生化代謝途徑和信號轉導途徑，同時找出差異基因相對于所有有注釋的基因顯著富集的通路。

本研究利用散點圖的形式展示KEGG富集分析得到的結果。圖12橫軸表示差異表達的基因中位于該通路條目的基因數目與所有被注釋基因中位于該通路條目的基因總數的比值。在散點圖的右側Q-value的取值范圍是[0,1]。在散點圖中衡量KEGG富集程度有3 個指標，分別是比值、Q-value和富集到此通路上的基因個數。當比值越大，Q-value越接近于零則表示富集越明顯。

從圖12可以看出，兩個樣品的差異基因KEGG代謝通路中顯著富集的有α-亞麻酸代謝、PPAR信號通路、脂肪酸降解、不飽和脂肪酸的生物合成、纈氨酸、亮氨酸和異亮氨酸降解、過氧化物酶體、丙酸代謝、氯代烷烴和氯代烯烴降解等。

圖12 差異基因KEGG富集散點圖Fig. 12 KEGG enrichment scatter plot of differentially expressed genes

3 結 論

本研究采用Illumina HiSeq測序平臺，對以不同碳源生長的Y. lipolytica的轉錄組進行測序分析，測序質量良好。其中對照組是以葡萄糖為碳源的Y.lipolytica，共獲得18 223 407 raw reads，經過過濾后得到17 923 921 clean reads；實驗組是以油酸為碳源的Y. lipolytica，共獲得22 680 303 raw reads，經過過濾后得到22 656 852 clean reads。按照差異基因表達變化2 倍以上且FDR≤0.05的原則，兩個樣品之間共篩選出536 個顯著性差異表達基因，其中376 個基因上調表達，160 個基因下調表達。通過差異基因GO富集分析結果可以明確葡萄糖和油酸在細胞組分、生物學過程和分子功能三大分類中表達基因的差異。最后通過差異基因KEGG富集分析確定差異表達基因參與的最主要生化代謝途徑和信號轉導途徑，以散點圖的形式對分析結果進行展示，結果顯示兩個樣品的差異基因KEGG代謝通路中顯著富集的通路有α-亞麻酸代謝、PPAR信號通路、脂肪酸降解等。研究結果為進一步分析不同碳源對Y. lipolytica生長產生的影響提供了科學依據。