不同碳源培养下解脂耶氏酵母的转录组差异分析

常 晨，刘伟丰，郭 阳，包怡红,*

（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生长产生的影响提供了科学依据。