|
http://www.abbs.info e-mail:abbs@sibs.ac.cn ISSN
0582-9879
ACTA BIOCHIMICA et
BIOPHYSICA SINICA 2003, 35(8): 707–716
CN 31-1300/Q |
03138
Estimating Coarse Gene Networks from Yeast Gene Expression Time Series
XU
Xiao-Jiang#, WANG Lian-Shui#, DING Da-Fu*
( Key Laboratory of
Proteomics, Institute of Biochemistry and Cell Biology, Shanghai Institutes for
Biological Sciences, the Chinese Academy of Sciences, Shanghai 200031, China )
Abstract Gene
networks is the collection of gene-gene regulatory relations at the expression
level. In this study, a combined approach of the linear transcriptional
modeling, identification of promoter elements and gene co-expression clustering
is developed to decipher yeast gene networks from expression time series. The
cell must reorganize the genomic expression to programs required for growth and
survival in each environment. The expression of many genes is regulated by
environmental stress. The products of many genes that induced in the
environmental stress are involved in metabolism of carbohydrates, structural
repairs and even sporulation. Interestingly, it is identified that
transcription factors Mcm1 and Dal82 matched their binding sites
TT-CC---T--GGAAA and TGAAAAWTTT in cell cycle progression and environmental
stress response, respectively. These conclusions agree with the known
observations. The results indicate that the approach may be useful for modeling
gene networks from microarray data.
Key
words gene regulatory
network; genome time serial expression; cell cycle; environmental stress
________________________________________
Received:
April 28, 2003 Accepted:
May 16, 2003
This
work was supported by the grants from the National High Technology Research and
Development Program of China (863 Program) (No. 2002AA234021), Knowledge Innovation
Program of the Chinese Academy of Sciences (No. KJCX1-08), and Shanghai Science
and Technology Commission (No. 00JC14018)
# These authors contributed
equally to this work
*Corresponding
author: Tel, 86-21-54921254; Fax, 86-21-54921011; e-mail, dingdafu@server.shcnc.ac.cn
从酵母表达时间序列估计基因调控网络
徐肖江# 王连水# 丁达夫*
( 中国科学院上海生命科学研究院生物化学与细胞生物学研究所蛋白质组学重点实验室, 上海 200031 )
摘要 基因调控网络是生命功能在基因表达层面上的展现。
用组合线性调控模型、 调控元件识别和基因聚类等方法, 从基因组表达谱解读酵母在细胞周期与环境胁迫中的基因调控网络。 结果表明,
细胞在不同环境条件下会调整基因调控网络。 在适应的环境下,
起主要作用的是细胞生长和增殖有关的基因调控网络; 而在响应环境胁迫时, 细胞会再规划调控网络,
抑制细胞生长和增殖相关的基因, 诱导跟适应性糖类代谢与结构修复相关的基因, 还可能启动减数分裂产生孢子。 分别从细胞周期和环境胁迫响应相关基因中, 搜索到转录因子Mcm1结合位点TT-CC---T--GGAAA, 和Dal82在尿囊素代谢途径相关基因上的结合位点TGAAAAWTTT。 从而,
从酵母表达时间序列估计基因调控网络是可行的, 与至今已知的实验观察相当吻合。
关键词 基因调控网络; 基因组表达时间序列; 细胞周期;
环境胁迫
细胞生化网络是生命系统展现功能的基础, 至少有三种层次:基因表达、 代谢和信号转导。 基因表达不仅直接由转录因子调控而且还跟代谢网络与信号转导偶联; 酶与蛋白质又是基因的产物,
因此当投射到基因空间时, 生化网络就约化为基因网络[1]。 目前已有一些方法从cDNA微阵列测量估计基因网络,
例如布尔网络、 贝叶斯网络、 微分方程和线性模型(包括SVD奇异值分解)[2, 3]。 这些方法还有若干不确定的因素, 例如： 目前基因表达阵列测量中基因数(数千以上)远大于时间抽样点数(20左右),
因而无法唯一地确定网络参数; 这些网络模型仅描写基因之间的视在相互作用, 不区分直接的(因果的)相互作用还是间接相互作用。 前者,
指一个基因与编码转录因子的基因之间的关系; 后者, 例如某代谢物影响某基因表达, 而另一基因编码调节此代谢物的酶。 目前,
转录因子与顺式元件之间的相互作用可以利用“转录因子／顺式元件”数据库或者全基因组转录因子结合位点资料[4]给予辨认。 基因组表达的时间抽样不足的问题除了今后充分增加样点数外, 基因按共表达聚类仍能给出基因调控网络的有用估计。 这也显示,
疾病基因调节关系的解读也有助于药物目标的识别和药物设计, 从而为复杂疾病的诊断与治疗提供有价值的线索。
本文组合线性调控模型、
调控元件识别和基因聚类等方法来解读酵母在细胞周期与环境胁迫中的基因调控网络。
1 材料和方法(Materials
and Methods)
1.1 材料
cDNA微阵列(microarray)可测量基因组转录表达数据; 若测量遍及一个时段,
即在若干时间点上测量, 就形成基因组表达时间序列, 而每个基因的表达时间序列又称基因表达谱(expression profile), 用此表达时间序列估计基因调控网络。
作为研究案例, 本文采用三套酿酒酵母(Saccharomyces cerevisiae)基因组表达时间序列。 一是细胞周期(ALPHA,
α因子同步化)表达时间序列, 含有6178个基因, 18个时间采样点,
取自Spellman等[5]的实验;
二是细胞响应环境胁迫的表达时间序列, 含有6152个基因,
包括： 从25
℃到37 ℃的热激, 有8个时间点; 细胞经过氧化氢处理,
有10个时间点; 后二套取自Gasch
等[6]的实验。 酵母基因的序列信息及其功能提示采用数据库SGD
(http://genome-www.stanford.edu/Saccharomyces/)。
1.2 方法
1.2.1 数据预处理 在基因表达时间序列中, 若某基因的缺失值的比例超过20%, 则被剔除[7]。
若某基因的表达谱的绝对值都不超过2(阈值), 亦被剔除,
因为表明它不参与调控过程[3]。
对细胞周期情况, 此阈值取2; 对于环境胁迫情况,
此阈值取3。
1.2.2 线性网络模型线性 网络模型假定基因之间的相互作用是线性的和非瞬时的, 即基因i在时刻tk+1的表达水平是前一时刻tk所有基因j(j=1,
2,…,n)表达水平的加权相加[公式(1)]。
yi(tk+1)=
, tk=t0+kΔt , k=0, 1, 2, …, T–1
(1)
此处Rij表示所有基因j对基因i的调控强度, Δt表示相互作用的平均传递时间[3]。 当调控网络参数给定, 式(1)描写基因网络的动力学, 这是工程学中常见的;
当基因表达时间序列yi(t)给定, 要估计调控网络参数Rij,
现称逆向工程问题(reverse engineering)。
本文的目的是后者。 现已有的几篇论文用“奇异值分解”方法求解过此逆向工程问题[8]。 然而, 目前的测量均限于n>>T,
故他们的解集是不惟一的, 因而难有生物学意义。
本文试图先将众多基因聚类, 寻求类类之间的调控网络, 而类内的基因相互作用用其他知识或信息来辨认。
目前, 常把具有相似表达谱的共表达基因归属一类, 同类中的平均基因表达谱定义为代表此类的原型基因(prototype)的表达谱。 原型基因的调控网络形式上仍为式(1), 但参与网络的基因数不是n, 而是原型基因的总数P(即聚类总数)。 逆工程问题就是要寻求网络参数Rij, 使得用线性网络模型拟合原型基因表达时间序列的残差总量[公式(2)]达到最小。
min∑T-1〖〗k=0∑P〖〗i=1[yi(tk+1)-∑P〖〗j=1Rijyj(tk)]2(2)
当上述时间采样点是满值或无缺值, 则式(2)有解析解,
可参见文献[9]。
当时间采样点中有缺值, 像测量点不等间距, 例如tk+1与tk相隔两个Δt,
只要将式(2)相关项改成公式(3), 其他项不变。
[yi(tk+1)-∑P〖〗j=1(∑P〖〗a=1RiaRaj)yj(tk)](3)
此时涉及Rij的非线性组合, 式(2)没有解析解, 但可用最优化方法,
例如遗传算法[10], 求得数值解。
1.2.3 基因聚类对预处理后的基因表达谱进行归一化, 即表达谱减去其均值再除以其标准差; 并用线性插值法填补缺失数据。 以基因表达谱的皮尔逊(Pearson)相关系数作为共表达聚类的量度, 其基因的分级聚类、
显示和分析通过软件包“Cluster”和“TreeView”完成[7]。
1.2.4 调控元件与转录因子的搜索原型基因所含的基因是共表达的, 常常也是共调节的, 可用软件包AlignACE[11],
搜索基因上游的启动子序列, 寻求所有可能的顺式调控元件或转录模体(motif)。
再用酵母启动子数据库SCPD核查已知转录模体和相应的转录因子[12]。
还可用TRANSFAC 6.0数据库寻觅转录因子功能信息[13]。
2 结果(Results)
2.1 细胞周期调控网络
我们找到563个在细胞周期过程显著表达的基因, 在此分为6类[图1(A)], 线性模建得到6个原型基因之间的调控网络[图1(B)]。 由表1可知,
原型基因1在M/前G1期有表达峰值, 主要是与细胞交配相关的基因。
原型基因2在M/G1期显著表达,
主要是糖类代谢、 氨基酸合成与分解、
能量代谢相关基因, 核糖体蛋白和rRNA合成相关基因;
从这些基因中, 我们搜索到转录因子Gal4、
QBP、 Bas2、 Pho4、
Dal82、 Gcn4 (调节糖类、 氨基酸和核酸代谢相关基因的转录)的DNA结合位点(见表2)。 原型基因3在G1期有表达峰值, 主要是细胞周期控制、
DNA合成和芽孢生长相关基因。 原型基因4表达峰值位于S(组蛋白基因)和G2期及M前期, 其功能主要是完成染色体复制并促进细胞进入M期。 原型基因5表达峰值位于M后期, 主要是细胞营养物质合成相关基因。
原型基因6主要是M后期和G1前期表达, 主要参与染色体预复制复合物形成及M/G1期转换相关基因[5, 14]。
Fig.1 Graphical representation of the
linear model for 6 clusters on data set ALPHA
(A) The column of plots on the left hand
side depicts 563 genes that remain in the data set ALPHA after preprocessing,
grouped in 6 clusters. The column of plots on the right hand side depicts the
associated prototypes. (B) A network of 6 prototypes. The circles represent
prototypes. Lines show interaction, with red arrows standing for activation,
and blue bars for inhibition. Solid lines represent strong effect, and dashed
represent weak effect. a, mating; b, matabolism; c, budding & DNA
synthesis; d, mitosis control; e, nutrion; f, pre-replication.
Table 1 Members of prototypes and their description of gene
functions on the data set ALPHA
|
Cluster |
Description of function |
Member of Prototypes |
|
1 |
Cell growth, division and DNA synthesis |
SST2, FIG1, FIG2, AGA1, AGA2, GFA1, FUS1, FUS2, FUS3, KAR3, KAR4, KAR5,
AFR1, MID2, SLT2, CIK1, MID2 |
|
Signal pathway: |
SST2, FUS1, FUS3, CMK2 |
|
|
Transcription: |
HAP4, IME4, KAR4, TYE7 |
|
|
2 |
Cell growth, division and DNA synthesis: |
CDC48, BUB2, CIN5, TFS1, TOR2, STH1,
UTH1, PAM1, DIG2, SNF6, MFA1 |
|
Carbohydrate, Amino acid, Nucleotide, Energy metabolism: |
ICL2, KNH1, GLC3, TSL1, GSY2, PYK2, GLK1, GPH1, HXK1, PGM2, GPD1,
MUC1, ATH1; SNF6; MET6，MET28, MET31,ARG4, GAP1, ADE2, IRA2, MSU1, CEM1, STF2, ATP17, ALD6 |
|
|
Transcription: |
PHD1, CIN5, MET31, MET28, SNF6, STH1,
MDR1, SRB7 |
|
|
3 |
Cell cycle control: |
SWI4, CLN1, CLN2, CLB6, MCD1, RAD53,
PCL5, ZIP1, SMC3 |
|
DNA synthesis, repair and replication: |
CLN2, CLB6, RAD27, RAD51, RAD53, RAD54, RFA1, RFA2, RNR1, RNR3, POL12,
POL30, MSH2, MSH6, CTF4, DPB2, CDC9, PRI2, DUN1, TFB1 |
|
|
Transcription: |
SWI4, ASF1, MIG1, HCM1, SPT21, GAT1,
TAF25, TFB1, YOX1, HFM1，ABF1 |
|
|
4 |
Cell growth, division and DNA synthesis |
CIS3, BUD3, CHS2, MYO1, BUD4, SWI5,
ACE2, ECM33, CDC5, CLB1, CLB2, ASE1, CYK2; |
|
Carbohydrate, Amino acid, Energy metabolism |
MUP1, CAR2, MET17, MET14, ARO9, DIP5, MUP1, ALG7, PHO4, PMI40,
ECM15, ALG5, PSA1, PMT2, CHS2, ATF2, IDH1, ACE2, ERG4, ELO1, ERG3 |
|
|
Transcription |
DOT4, PHO4, HTA2, HTB2, HHT2, HHF2, HTB1, HTA1, HHT1, HHF1, HST3,
SWI5, ACE2，FKH1, FKH2, NDD1 |
|
|
5 |
Phosphate metabolism |
MIS1, CHA1, FUN63, URA7, PHO3, PHO5,
PHO11, PHO12, EXG1, MSN1, AUR1, FAA4 |
|
Cell cycle control |
DBF2, CDC47, YTM1, DIP2, MAK16, RME1,
SSF1, EXG1, IQG1 |
|
|
Transcription |
DBF2, MSN1, RME1, NSR1, NOP2, HCA4,
MRT4, SSM4, SRP40, GCD14 |
|
|
6 |
Cell cycle control and mitosis |
PCL2, FAR1, CLN3, PCL6, CDC54, CDC6,
PCL9, EGT2, SIC1, TPD3 |
|
Cell growth, division and DNA synthesis |
EGT2, TPD3, FAR1, HSP150, BUD9, CHS1, MFA2, STE2, GPA1, NUT2,
ASH1, XRS2, CDC54, MCM3, MCM2, CDC6, CTS1 |
|
|
Transcription |
NUT2, ASH1, PCL9, GCR3, TPD |
Table 2 Transcriptional motifs from 6 prototypes on data set
ALPHA
|
Motif |
Cluster |
Notes |
|
GGGGC-SSS-S----G |
C2 |
Gal4->GAL7; Gal4->GAL10; QBP->GAL10; Pho4->PHO8; IRE->IME1; Gal4->GAL80 |
|
G-----RR----GG-KGG-G |
C2 |
Gal4->GAL10; Dal82->DAL7; ARC->ARG1 |
|
AAA-G-GGSGC--S |
C2 |
Dal82->DAL4 |
|
GS-S-SSCC-----GGG |
C2 |
QBP->GAL10; heat shock not HSE->DDR2; BUF->CAR1 |
|
GS-GS-S--G---KGG |
C2 |
QBP->GAL1; URSSGA->SGA1 |
|
-G--G-G--SC--G---GG----C |
C2 |
QBP->GAL1; URSSGA->SGA1 |
|
GAAAAGCCC-T |
C2 |
Leu3->LEU2; URS1ERG11->ERG11 |
|
WRACGCGT-R |
C3 |
Mcm1->CLN3 |
|
CGCG-AA-A-AAA |
C3 |
Dal82->DAL4 |
|
CSGC---GG------SM-C--R |
C4 |
QBP->GAL1; Hap1->CYC1 |
|
GCGATGAG-TG |
C5 |
BUF->HO |
|
TT-CC---T--GGAAA |
C6 |
Mcm1->CLN3，CDC47，SWI4，FAR1，HSP150，CDC6，STE6，CCP1，CDC46 |
|
RARCCAGCMR |
C6 |