基因表达数据分析实验指导

1. 实验基本情况 2. 实验方法：

2.1 表达谱数据的下载

2.2 将表达谱数据导入matlab软件 2.3 补缺失值 2.4 数据标准化

2.5 差异表达基因筛选 2.6 选择差异表达的基因

2.7对差异表达基因送入功能注释

附 -- Matlab的Microarray Data Analysis

1. 实验基本情况

实验目的：

掌握和了解常用的基因表达分析过程，包括数据下载、数据预处理、差异表达分析和基因功能注释。了解GEO、SMD、Matlab软件和WebGestalt数据库的使用。 实验方法：

详见下面的描述。 实验作业：

每位同学从GEO或SMD数据库上下载一套表达谱数据，进行数据预处理，差异表达基因分析或聚类分析等数据分析过程（依据具体问题操作，arraytool或matlab或其他软件均可），基因功能注释（WebGestalt、GO、KEGG等数据库）。

实验实例分析

=====================================================================

2. 实验方法：

2.1 表达谱数据的下载

2.1.1 从GEO数据库上下载表达谱数据

1) 网址及数据库概述

GEO主页：http://www.ncbi.nlm.nih.gov/geo/

GEO数据库中包含四种类型的条目，分别以GPLXXXX（检测平台），GSMXXXX（生物样本），GSEXXXX（基因表达系列），GDSXXXX（基因表达数据集）表示。其中GPLXXXX有SAGE、MPSS、单色芯片（Affymetrix）、双色芯片（spotcDNA/DNA）几种；GSEXXXX与GDSXXXX的区别在于：GSE是实验者一次一起提交的数据集，包含原始的数据文件，而GDS是GEO数据库的维护者根据样本和实验平台的特性进行整理的，与原有的GSE数据可能有样本量上的差异；一般GDS都有对应的GSE数据；GDS不包含单独的原始数据，

如果想获得其原始数据，需要链接到他的GSE网页上下载；GDS样本间的可比性更强，如果有GDS就先分析GDS。 2）数据下载

GEO可提供两种数据的下载，一种是整理好的soft格式数据，是一个数据矩阵，包含多个基因在多个条件下的表达值，如GDS2220.soft；另一种是单独的数据文件，每张芯片一个数据表格，如GSE3519_family.xml文件夹里的文件，就是对应GDS2220这次实验的原始数据。另外还有一个GDS2220.annot数据是提供基因描述的。具体的下载方式如下：

在GEO主页上（图1），可以通过浏览（browse）或query中输入疾病名字，如风湿性关节

炎（rheumatoid arthritis）在Datasets后，点击go进行搜索，结果如图2。

图1. GEO的主页

图2.GEO的搜索结果

之后点击你感兴趣的GDS集合，如GDS2220，就进入每套数据的页面了（图3）。

图3.GDS2220数据的浏览界面

在图3中，点击下拉菜单中的DataSet SOFT file，就能下载GDS2220.soft文件；点击Annotation SOFT file就可以下载GDS2220.annot文件；点击seriers family miniml file就可以下载GSE3519_family.xml文件夹,但这个速度较慢，这里有个小窍门，大家可以在迅雷中新建一个下载任务，粘贴地址： ftp://ftp.ncbi.nih.gov/pub/geo/DATA/MINiML/by_series/GSE139/GSE139_family.xml.tgz ，这里GSE139是可以替换的，比如要下载GDS2220配套的数据，就直接把两个GSE139都替换成GSE3519就可以直接下载了；点击series family soft file下载的文件与GDS2220.soft类似，只是样本是GSE3519的数据，可能和GDS2220的样本不同，这里是相同的。 也可以通过以下方式寻找特殊平台的数据。

3) 文件描述

（a）GDS22.soft

该文件从上到下分为三个部分：第一部分，数据集合基本描述，文字形式，以！或#开头；第二部分，表格的表头，如“ID_REF IDENTIFIER GSM80309 GSM80310 GSM80311 GSM80312

GSM80313 GSM80314 GSM80315 GSM80316 GSM80317”，以tab键分割，表示下面的数据部分每一列的含义；第三部分，数据，如GDS2220.soft中第一列为每一个基因的编号，第二列是基因名字，第三列是GSM80309 样本中个基因的表达值。从图3中我们可以看出，这个表达值是“log2Ratio”； （b）GSE3519_family.xml

从文件夹中可以看到，每个文件的名字很有规律，包含三类文件：第一个是GPL2715-tbl-1.txt，表示的是该实验检测平台的信息，第一列为基因的编号（与GDS2220.soft文件的第一列编号相同），第二列是基因的名字，最后三列分别对应每个基因在芯片上的block，row和column号。详情可参见：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL2715

2.1.2 从SMD数据库上下载基因芯片数据

主页：http://genome-www5.stanford.edu/， 如图4

图4. SMD数据库主页

SMD数据库也提供两种类型的数据：一种是与GDSsoft格式类似的整合后的表格，如7742.pcl，另外一种是单独的数据原始文件，与GSE3519_family.xml类似，文件如exptsetno_3977.tar.gz。 1）数据下载流程:

点击SMD主页（图4）上的Publication，进入数据浏览页面（图5），可进行适量的物种或发表年限筛选，选中你要下载的数据后，点击对应行的最后一列的SMD图标，进入图6的下载页面，点击view可以看见每个样本的描述，点击Rawdata可以下载原始的数据文件如exptsetno_3977.tar.gz，由于数据量太大，这里不建议下载。点击data retrievel and analysis，可以进行重复数据合并、缺失数据处理等过程，不用更改直接默认选项就可以了，最后网站会给你一个数据下载的链接，你点击就可以下载7742.pcl文件了(图7).

图5. SMD的数据下载页面 图6. 具体数据的下载页面

图7. 合并后的数据下载，点击download preclustering file

2）例子文件

共包含39729个唯一的基因的表达数据，50个样本，其中35个为RA患者，15个为正常对照组。 下载地址：

http://smd.stanford.edu/cgi-bin/publication/viewPublication.pl?pub_no=664

对应的文件分别为7742.pcl和exptsetno_3977.tar.gz。 3）文件描述

7742.pcl，第一行是表头，第二行是每一列的权重（这一列不是基因表达值，分析时要去除），第三行开始是数据。每行数据的第一列是芯片上的编号，第二列式基因的名字，这里包含多个数据库中的名字和ID，用||分隔开，第三列是基因的权重，也不是表达值，第四列开始是表达值了对应于SHFM212样本，这里每个值是log2Normorlizedmeanvalue，这个在数据整合时看见了。

2.2 将表达谱数据导入matlab软件

Matlab中数据有三种储存格式，包括结构型数据（如maStruct.mat），cell型数据(如yeastdata.mat 中的genes文件，储存的就是基因名字)和数值型数据(如yeastdata.mat 中的yeastvalues文件，储存的就是基因的表达值)。所有matlab文件都可以用 “load 文件名”的形式导入matlab文件，但首先要把我们下载的文本或excel表格文件转换成matlab可以识别的文件。

2.2.1 GDS2220.soft文件的导入

1）更新geosoftread.m程序

由于matlab对soft格式文件导入的程序有一点小bug，我们首先要将matlab文件更新。 关闭matlab，打开c:\\Program Files\\MATLAB71\oolbox\\bioinfo\\microarray，将geosoftread.m文件粘贴至此文件夹，将原有文件替换。打开matlab程序，在命令行上输入： rehash toolbox

后，新的geosoftread程序就导入了。

2）将matlab的当前目录转换为GDS2220.soft文件储存的目录 在命令行中输入：

gdsdata = geosoftread('GDS2220.soft')

就将GDS2220.soft文件导入为matlab可识别的gdsdata结构文件了，该文件是一个树形文件，还包含8个子文件，想获取任何一个子文件，只需在gdsdata后加.加子文件名字即可，如你想对基因表达矩阵分析，就用gdsdata.Data，就是一个16515*20的数值矩阵，代表GDS2220.soft文件中的表达值部分。

七个子文件描述如下：（%后为我添加的描述信息） gdsdata =

Scope: 'DATASET' %数据类型，是数据集（GDS）、样本（GSM）、系列（GSE） Accession: 'GDS2220' % GEO中的编号 Header: [1x1 struct]

ColumnDescriptions: {20x1 cell} %每一列的描述，其实就是每个样本的描述 ColumnNames: {20x1 cell} %每一列的名字，其实就是每个样本的名字

IDRef: {16512x1 cell} %每一行的编号，其实就是每个基因的编号，对应.soft文件的第一列 Identifier: {16512x1 cell}%每一行的描述，基因的名字 Data: [16512x20 double] % 基因表达矩阵，数值型数据

最后将你导入的文件保存成matlab文件，以后就可以直接打开matlab文件分析了。 存储文件的命令： save gdsdata gdsdata

2.2.2 对GPR格式的数据输入

maStruct = gprread('mouse_a1wt.gpr');

mouse_a1wt.gpr是matlab自带的一套数据，单张芯片包含原始的红光和绿光值，与我们从GEO上下载的原始数据类似，但是我们的GEO数据不可以直接用gprread命令导入。我们可以参照maStructure的结构，生成GEO导入的数据。 maStructure的结构如下： maStruct =

Header: [1x1 struct]

Data: [9504x38 double] %每一行代表一个基因在一次芯片检测中的各种值，cy3，cy5等。 Blocks: [9504x1 double] %每个基因所在的区域 Columns: [9504x1 double] %每个基因所在的列 Rows: [9504x1 double] %每个基因所在的行 Names: {9504x1 cell} % 基因标识 IDs: {9504x1 cell} % 基因编号

ColumnNames: {38x1 cell} %每一列的含义，对应Data的38列 Indices: [132x72 double] % 每个基因在芯片上的坐标 Shape: [1x1 struct] %每个block在芯片上的起始坐标

提示：若想将GEO的family数据导入，则需将单独的数据样本和单独的GPL平台文件同时导入，才可生成上述的maStruct结构文件。此部分本实验不做要求，同学可自行尝试编写程序。SMD的rawdata导入时也需要编写程序 2.2.3 SMD的整合好的数据导入（7742.pcl）

右键单击7742.pcl文件，选择用excel程序打开；

删除第二列（基因名字那列），因为这列太长了，导入matlab时容易产生错误； 存储为7742.xls

打开matlab，点击file下拉菜单，选择import，选取7742.xls就可以导入了。 生成三个文件

a) Data，对应原始的数据表，但要注意第一行和第一列不是表达值 可以通过命令去除： data(:,1)=[];data(1,:)=[];

b) textdata, cell型文件，对应基因的ID和表头 c）colheaders， 无用

2.3 补缺失值

使用k均值法补缺失值

原始命令：Imputedata=knnimpute(Data, k);

这里Data的是我们的表达谱，k为k近邻的k值，默认的距离测度是欧式距离。 如对gdsdata进行15近邻补缺失值： Imputedata=knnimpute(gdsdata.Data, 15); 可将补后的数据并入原来的gdsdata中 gdsdata.imputedata=Imputedata;

2.4 数据标准化

2.4.1 对soft格式文件的标准化

1）导入gdsdata load gdsdata

2）显示各张芯片间的变异(如图8)

maboxplot(gdsdata. imputedata,gdsdata.ColumnNames); xlabel('Samples');

图8.各芯片间的变异

3）标准化

a）使用manorm命令可以对每张芯片的几种趋势进行标准化 Xnorm = manorm(data) 默认是对均值进行标准化

b)使用quantilenorm直接对集中趋势和离散趋势进行标准化，并用图示显示： 默认为均值标准化，如改为中值标准化：

NORMDATA = QUANTILENORM(gdsdata.imputedata,'MEDIAN',1) ; 显示标准化后的结果（图9）

Normdata = quantilenorm(gdsdata.imputedata,'display',1);

图9. 标注化后的数据

将normdata并入gdsdata中 gdsdata.normdata=normdata; 2.4.2 其他类型数据的标准化

详见Bioinformatics Toolbox 中的Microarray Data Analysis，在文档最后。

2.5 差异表达基因筛选

T检验筛选

[h,significance,ci] = ttest2(x,y,alpha);

这里x是一个基因在正常样本中的表达值，y是该基因在另一类样本中的表达值，alpha是显著性水平，通常取0.05. h是0或者1，表示h0假设（无差别）还是h1假设（有差别），这是依据计算出的p-value和你的alpha比较来判断的。Significance是p-value。Ci是置信区间。

这里不要求两个向量长度相等，每次只对一个基因计算。

如果我们相对gdsdata中的16512个基因计算p-value，就要编写循环，如 for i=1:size(gdsdata.normdata,1)

[h,significance,ci] = ttest2(gdsdata.normdata(i,1:10), gdsdata.normdata(i,11:20),0.05); Sigtotal(i)=significance; end

Sigtotal文件对应的就是每个基因的显著性水平了。

2.6 选择差异表达的基因

根据你的判断，如p小于0.05的基因取出来，作为差异表达基因。 Siggene= gdsdata.Identifier(Sigtotal<=0.05);

2.7对差异表达基因进行功能注释

可以先使用我的帐号，但分析完请把你输入的数据集删除

Your username is your email address: ******************** Your new password is: VaRi7esi

We recommend you save this email for future reference.

The WebGestalt web site is http://bioinfo.vanderbilt.edu/webgestalt/

附 -- Matlab的Microarray Data Analysis

Bioinformatics Toolbox Microarray Data Analysis

MATLAB is widely used for microarray data analysis. However, the standard normalization and visualization tools that scientists use can be difficult to implement. The Bioinformatics Toolbox includes these standard functions.

Microarray data — Read Affymetrix GeneChip files (affyread) and plot data (probesetplot), ImaGene results files (imageneread), and SPOT files (sptread). Read GenePix GPR files (gprread) and GAL files (galread). Get Gene Expression Omnibus (GEO) data from the web (getgeodata) and read GEO data from files (geosoftread).

Microarray normalization and filtering — The toolbox provides a number of methods for normalizing microarray data, such as lowess normalization (malowess) and mean normalization (manorm). You can use filtering functions to clean raw data before analysis (geneentropyfilter, genelowvalfilter, generangefilter, genevarfilter), and calculate the range and variance of values (exprprofrange, exprprofvar).

Microarray visualization — The toolbox contains routines for visualizing microarray data. These routines include spatial plots of microarray data (maimage, redgreencmap), box plots (maboxplot), loglog plots (maloglog), and intensity-ratio plots (mairplot). You can also view clustered expression profiles (clustergram, redgreencmap). You can create 2–D scatter plots of principal components from the microarray data (mapcaplot).

Microarray utility functions — Use the following functions to work with Affymetrix and GeneChip data sets. Get library information for a probe (probelibraryinfo), gene information from a probe set (probesetlookup), and probe set values from CEL and CDF information (probesetvalues). Show probe set information from NetAffx (probesetlink) and plot probe set values (probesetplot).

The toolbox accesses statistical routines to perform cluster analysis and to visualize the results, and you can view your data through statistical visualizations such as dendrograms, classification, and regression trees.

蛋白质分析实验指导

实验一：

进入Pfam数据库的主页：http://pfam.sanger.ac.uk/。关于Pfam基本信息，请参考我提供的word文档。选择KEYWORD SEARCH，输入lipocalin，得到结果如下：

为了研究的方便，我们没有选择Lipocalin family，而是采用PF08212为例进行讲解，请同学们按照同样的流程以PF00061为例重新做一次。打开PF08212，得到如下内容：

从左侧可以看到该项目有summary、Domain organization等几项组成。请根据网页说明，弄清分别代表什么含义，并弄清上方方框内各个数字的含义，比如“4 architecture”代表什么。然后打开Alignments 并选择View：

得到种子序列的多重比对形式：

这里的结果显示GXW模式非常保守。在applet窗口中查看序列比对并进行分析和保存（calculat）这个applet还能显示比对家族的主成分分析结果：

这里，根据一种距离标度将多序列比对中的每一个蛋白质表示为空间中的一个点，因此很容易鉴别不属于同一类的蛋白质。可以用鼠标控制坐标轴的旋转。很明显有两个载脂蛋白离群很远，选定任一蛋白，在JalView中的多重比对中的相应序列会高亮显示。类似的信息也可以用树的形式子在applet中表示出来：

可以看出，最上面有两个序列与其他序列是分离的。

三个applet窗口是相联系的，同样的颜色可以分别在三个applet窗口中表示同一条序列，比如下图：

两条绿色的蛋白质离其他序列最远，在平均距离树中也是绿色，同时多重比对中的绿色显示了这两个蛋白质的序列信息。请多动手操作，熟悉Pfam数据库的这些功能和其他信息。

实验二：

打开第二张图左侧的Structures，得到下图：

可以看到，大部分物种的蛋白质都对应一个PDB ID，这意味着该结构域或蛋白质在PDB数据库具有结构信息。为了进一步了解结构域、蛋白质结构信息以及与蛋白质序列之间的联系，大家需要了解一下PDB数据库。关于蛋白质结构和蛋白质组学的相关内容请参看我提供的word文档。

打开http://www.rcsb.org/pdb/home/home.do 进入PDB。PDB获取码由一个数字和3个字母组成，同时PDB也支持关键词搜索。

搜索牛增味剂结合蛋白（bovine odorant binding protein）或获取码1PBO，请观察返回的页面主要包括哪些信息：

注意网页下部有，Pfam数据库和Gene Ontology的信息：

PDB数据库在结构生物学中占有中心地位，在NCBI中也可以获得蛋白质结构数据，即可以在NCBI中访问PDB记录。打开http://www.ncbi.nlm.nih.gov/sites/entrez?db=Structure：

通过关键词或PDB获取码获得增味剂结合蛋白1pbo，

点击domains中的3Ddomains可以得到同一个蛋白质不同链的记录，而点击1PBO则显示MMDB中关于这个蛋白质的记录：

下载Cn3D，并点击Structure View in Cn3D下载视图信息，然后用Cn3D打开。会打开两个窗口，一个是Cn3D Viewer，一个是OneD-Viewer：

Cn3D Viewer可以用其中格式来显示氨基酸序列,同时可以对结构进行旋转。相应的OneD-Viewer则显示此蛋白质的氨基酸序列。在Cn3D Viewer和OneD-Viewer中高亮显示某个或某一组氨基酸，在另一个查看器中的相应部分会自动高亮显示。

下面是各种形式的增味剂结合蛋白质：（worms-螺纹状模型, tubes-管状模型, wire-线状模型, ball and stick-球棍模型, space fill-填充模型）

在窗口中点击保守模体GXW的序列，蛋白质结构中的相应部分将会高亮显示，除此之外还有很多可视化选项，如显示结合位点的配体分子图像。请同学们多多发挥“主观能动性”，去探索这个软件，并通过它认识蛋白质组成。

回到MMDB的“Structure Summary”页面中找到“Related Structure”,点击后面的VAST，得到：

再点击Entire chain得到：

这个列表实际上是蛋白质矢量比对搜索工具（VAST）返回结果的一部分。除了当前的增味剂结合蛋白外，再在这个列表中选择其他蛋白，比如视黄醇结合蛋白和乳球蛋白，然后点击”View 3D Alignment”获得视图文件，利用Cn3D打开，就可以看到3个蛋白质的结构了，同时相应的多序列比对会显示在二维查看器中。如下图：

除了这个，还有其他很多基于网络的结构可视化工具，如： Swiss-PDB Viewer Chime RasMol MICE VRML

这些资源，大家如果将来做蛋白质结构研究，可以多了解一下。

计算表观遗传学实验指导

实验内容：

人类组蛋白修饰数据库（HHMD），网址：http://bioinfo.hrbmu.edu.cn/hhmd。组蛋白修饰在染色质重塑、转录调控和人类疾病中起着重要的作用。当前基于ChIP技术的实验测定了人类基因组高通量的修饰位置信息，人类组蛋白修饰数据库是迄今为止收录各种实验测定的人类基因组组蛋白修饰最为全面的数据库，共涵盖了43种人类组蛋白修饰的大通量实验数据，并提供了通过文献挖掘得到的9种癌症相关的组蛋白修饰的信息。目前，该数据库在不断更新中。 实验目的：

• 学习HHMD的搜索方式 • 学习HHMD的可视化方式

• 按照实际需要挖掘组蛋白修饰模式 实验方法：

• 熟悉HHMD的网站结构。

• 参考网站的帮助页面http://202.97.205.61/hhmd/help.jsp，利用该数据库提供的五种

搜索组蛋白修饰的方式，即组蛋白修饰类型、基因ID、功能注释、染色体定位、癌症类型，搜索相应的基因组区域中的组蛋白修饰。

• 参考网站的帮助页面http://202.97.205.61/hhmd/help.jsp，利用该数据库提供的三种

可视化组蛋白修饰的方式，即染色体图谱、基因ID、染色体定位，对相应的基因组区域中的组蛋白修饰进行可视化。

作业（5选4）：

1、 下载1p36.33的H3K9me3和H3K27me3的数据，且组织为淋巴（Lymph）。提交的内

容为提取数据的步骤（可辅以图说明）、结果截图（如跨越一屏，可用snagit软件等辅助）。及最终的数据行数（可用UltraEdit等编辑器数行数）。 2、 下载NM_153254的TSS上游1k及下游1k范围内的所有ChIP-seq测定的各组织组蛋白

修饰数据。提交的内容为提取数据的步骤（可辅以图说明）及最终的数据行数。

3、 对1p36.32的所有淋巴组织（Lymph）的数据进行可视化。提交的内容为操作步骤（可

辅以图说明）及结果截图。

4、 对NM_153254的TSS上游1k及下游1k范围内的ChIP-seq测定的H3K27me3进行可

视化。提交的内容为提取数据的步骤（可辅以图说明）及结果截图。

5、 对基因组区域：chr1：6300000-6400000的淋巴组织所有组蛋白修饰进行可视化，并在

可视化界面中下载数据。提交的内容为操作步骤（可辅以图说明）、结果截图及下载的数据的行数。

作业示例：

1、首先通过Search下的“Search by histone modification”配置搜索选项，如下图所示。

点击“Submit”，可得下图。

点击保存按钮后，可下载相应数据文件：共计3081行。

2、首先通过Search下的“Search by gene ID”配置搜索选项，如下图所示。

点击“Submit”，可得下图。

点击保存按钮后，可下载相应数据文件：共计212行。

3、首先通过HisModView下的“Start by Cytogenetic map”，点击1号染色体上的36.32区域（注意任务栏地址），配置搜索选项，如下图所示。

点击“Submit”，可得下图。

4、首先通过HisModView下的“Start by RefSeq Gene ID”，配置搜索选项，如下图所示。

点击“Submit”，可得下图。

5、

3、首先通过HisModView下的“Start by Chromosome Location”，配置搜索选项，如下图所示。

点击“Submit”，可得下图。

点击保存按钮，可下载到组蛋白修饰。经查，有224行数据。

生物学网络分析实验指导

内容一：gene ontology（GO） http://www.geneontology.org/

关于数据库的介绍，请参考提供的幻灯，其中部分为讲课中的内容。鼓励大家先用十分钟浏览一下数据库的大致资源，然后进行下面的内容。 1.

Amigo是进入这个数据库的基本工具。打开Amigo做简单的了解。

打开Amigo后再点击上方的Search，会出现高级搜索选项，了解一下这些选项的内容代表什么。

2.

上面的内容是针对单个蛋白质或少数几个蛋白质查询的，如果批量处理比如整个基因组的注释信息，需要了解下面的内容。

打开GO主页右侧中的Downloads中的Annotations（下图），从中可以看出GO是个外部数据库，依靠很多数据库到GO的注释。

下载SGD数据库对应的酵母物种的GO注释：

注释信息在annotations文件中，借助README文本了解annotations注释文件的结构，所有的GO注释都采用了类似的结构。这是功能注释的核心文件。

注意上图中出现了“IEA”，打开http://www.geneontology.org/GO.evidence.shtml了解各个证据的具体含义，重点是了解IEA。

在Downloads中的Annotations的下方，有Mapping to GO，打开看看是什么含义，可以下载里面对应的Pfam domains的信息，了解具体情况：

3.

为了能够本地处理GO数据库的文件，比如本地查找属于某个功能的所有酵母基因，除了需要下载注释文件，还要下载GO的Ontology文件，即了解GO的构成。GO的Ontology文件有很多类型，其中我们最了解的莫过于MYSQL格式。打开Downloads中的Ontology，寻找下图：

打开GOdatabase archives，寻找并下载其中的文件名类似“go_YYYYMM-termdb.tables.tar.gz” 的文件。

这个文件可以帮助你建立本地数据库（下图）。

本地化数据库实验课上不要求实现，我们的目的是了解其中的两个文件，这对于理解GO至关重要，一个是term，一个是graph_path，请下载后打开观看详情。表graph_path中的中间两列使用数字表示的，每一个数字对应term表中的第一列，因此对应的是一个GO term，或Ontology树状图中的一个节点。graph_path描述了任意两个存在一条路径的节点之间的关系，第四列表示第三列节点到第二列节点之间的距离。试想一下，怎么利用SQL语句和这两个文件实现广义的注释？比如，一个蛋白质注释到了叶绿体膜，它同时也具有叶绿体的功能，因此，怎么利用这两个文件确定一个蛋白质注释的所有功能呢？这是一个数学中的传递闭包的概念，请参看GO的帮助http://www.geneontology.org/GO.database.shtml中的“Traversing Graphs”了解具体情况。方便的话可以再了解一下“term2term”这个表格的含义。

因为所有的注释都只提供目前已知的最详细的功能，而通过网页来进行的研究又受到网速和网页服务的限制，所以为了快速和批量的处理功能注释，了解以上的几个Ontology文件是必要的。

内容二、互作网络 1. 下载数据

研究互作网络当然需要了解数据源，我们这次看一下DIP（Database of Interacting Proteins）：http://dip.doe-mbi.ucla.edu

这个数据库里包含了目前已知最全面的酵母的蛋白质互作网络。注册后下载酵母的互作网络，然后在Cytoscape中观察网络的性质。

上图为下载数据的列表，这个网页出现的过程相信大家很快就能找到。注意上面的数据分成两类，一类是full、一类是core。请浏览网络找出这两种数据的区别，并思考为什么构建core数据集。

2. 在DIP上了解数据中每一列代表什么，然后处理数据。

把互作数据简化处理成两列，位于同一行的两列表示两个蛋白质互作，然后把自环去掉。即A-A这种类型的互作去除。最好采用uniprotkb，比如用类似“P32561”的格式表示蛋白质。这是为什么呢？ 打开GO主页右侧中的Downloads中的Annotations，找到位于网页最下方的下图：

打开后下载酵母的gp2protein文件，并对比SGD到GO的注释文件。然后你会发现，使用uniprotkb可以更方便寻找蛋白质的功能注释。为了处理的数据量考虑，请选择core文件进行简化处理。

3. 进一步的分析之前，我们先对Cytoscape进行处理。

首先，安装插件，把BiNGO.jar和NetworkAnalyzer.zip的文件解压缩到Cytoscape的安装目录plugins下。利用help文件熟悉BiNGO.jar和NetworkAnalyzer两个插件：

BiNGO：http://www.psb.ugent.be/cbd/papers/BiNGO/

NetworkAnalyzer：NetworkAnalyzer-2.6.1-help.zip（见附件）

这两个文件在我提供的附件里都有，如果版本跟你的不兼容，请自行到Cytoscape的网站下载对应文件并安装。

把处理后的数据输入Cytoscape中并显示。如果不会处理和简化数据，请打开附件中提供的“yeast20071202.txt”，该文件为DIP以前提供的文件，格式非常标准，挑选第三列和第六列就达到了简化的目的。这个文件以orf表示蛋白质，打开SGD到GO的注释文件，就可以明白orf也可以直接进行功能注释。

上图是利用orf表示的蛋白质互作在Cytoscape中显示的局部。利用插件中的network analysis可以对蛋白质网络进行性质分析。

在help文件的帮助下，通过这个插件了解蛋白质网络的各个基本性质。 同样，还可以进行网络的功能富集分析，比如打开：

然后选择你感兴趣的节点进行分析，比如下面的高亮的节点：

下表反应的是高亮的这些节点显著富集的GO中的生物学过程（也可以选择其他分支或都选择）：

选择一个功能后，属于这个功能的网络中的节点就会显示出来。比如选择上表中的第二个功能后，网络显示：

这些功能在GO中的组织关系也会显示出来：

节点的颜色表示功能富集的显著性。

好了，关于Cytoscape我们就说到这里，我没有对细节做太多介绍。因为听说你们学过Cytoscape，那么理论上对细节可能比我掌握的更好，有欠缺的可以到课下在做补充。有时间的话还可以看看下面两个内容。 4.遗传互作。

biogrid： http://www.thebiogrid.org/ 从这个数据库下载遗传互作的数据

通过同样的方式，也可以对遗传互作数据做分析。 5. 下载致死性数据：

http://www-sequence.stanford.edu/group/yeast_deletion_project/deletions3.html

利用Cytoscape这个可视化软件，我们可以做很多事情。比如，把遗传互做中的hub节点以其他醒目颜色在物理互作网络中显示出来，观看它们的联系；把致死基因找出来，在物理互作中显示出来，观看它们的联系。

序列比对与BLAST实验指导

进入ncbi网站：http://www.ncbi.nlm.nih.gov/

实验一：腔棘鱼类Latimeria chalumnae被称为“活化石”。人们很早就认为它们灭绝于至少9000万年以前，现在在海洋中发现了几个标本。令人惊讶的是，对线粒体序列的分析表明腔棘鱼类比鲱鱼更接近与人类。在人类、腔棘鱼类和鲱鱼的线粒体中找一些DNA序列，并得到它们的accession number，然后执行双序列比对，看看结果是否如此。

提示：1.利用Entrez搜索物种的线粒体，比如鲱鱼可选择Genome然后搜索 clapea harengus Mitochondrion：

2. 进入下面的界面后找出一段DNA序列以做比对分析：

3.由NCBI主页进入BLAST（点击BLAST）：

然后选择BLASTN：

点击Align two or more sequences：

最后界面为：

4.可以输入整个线粒体的基因组编号，比如Clupea harengus mitochondrion 为NC_009577（如第二步图中所示），然后利用输入框右侧的范围限制你要比对的DNA序列。也可以直接输入从第二步中找到的某个基因的编号。当然，这里有两个输入框，分别输入人类线粒体编号和Clupea harengus的线粒体编号才能进行比较。试着看看，腔棘鱼类的线粒体序列是不是并不比鲱鱼关系更远。选择线粒体中的一个蛋白质进行比较，看看有什么分别。是不是蛋白质序列差异更小，效果不太明显。

做完上述工作之后，你已经对资源有了基本的了解，下面可以做个更有意思的事情（看时间选做，它们的序列可能不太容易找到）：

非洲的白氏斑马已经灭绝，这是一种既像马又像斑马的动物。在1972年，人们拍到了一张白氏斑马的照片（不是下面这张）。从博物馆的白氏斑马的标本里人们得到了它的线粒体的DNA并进行了测序。将白氏斑马（Equus quagga boehmi）、马（Equuscaballus）和斑马（Equus burchelli）的线粒体DNA进行比对，看看白氏斑马与谁更为接近？

实验二、

利用人类RBP蛋白质为例子了解PSI-BLAST和BLASTP的流程，该蛋白质编号为NP_006735。看看第一次的PSI-BLAST的循环所返回的结果和BLASTP是不是相同。注意每次BLAST的结果都可以保存，其中PSI-BLAST循环所产生的PSSM矩阵也可以保存。

了解基本流程后，可以做下面实验：寄生虫Plasmodium vivax有一个多基因家族vir为其所特有。这些基因有600～1000个拷贝，他们可能会通过抗原变异导致慢性感染。选择vir1在非冗余数据库中blastp搜索，然后用同样的输入进行PSI-BLAST搜索。 A、 在最初的搜索中，大约有多少个蛋白质的E值小于0.002，多少大于0.002？ B、 PSI-BLAST搜索第一轮和第二轮循环之间增加的最好的新序列的分值是多少？

实验三、

以人类视黄醇结合蛋白为例学习PHI-BLAST搜索。请选择PHI-BLAST，输入NP_006735进行搜索，第一次先不要输入PHI模式。然后再做第二次PHI-BLAST搜索，这次在输入框中加入下列模式：GXW[YF][EA][IVLM]。仔细观擦两次返回的结果有什么不同。PHI-BLAST的返回结果能不能作为PSI-BLAST的输入？这样做有什么好处？

尝试去prosite数据库中寻找一个感兴趣的家族，找出这个家族的信号序列，然后进行PHI-BLAST和PSI-BLAST搜索，看看能不能找到这个家族的新的成员。

实验四、

1. 创造一个人工蛋白质序列，包含两个不同的结构域。将这个联合序列进行BLAST搜索和PSI-BLAST搜索。自然界存在同时含有你选择的两个结构域的蛋白质吗？ PSI-BLAST能够在搜索多结构域蛋白质时，表现如何？

2. 挑选一个多结构域蛋白质，如HIV-1 pol（NP_057849），利用它来做blast和利用它的某一个结构域做blast有什么不同？PSI-BLAST搜索中又有什么不同？

序列信息资源实验指导

第一部分：课堂知识复习

形式：题目问答，请同学单独回答以下问题，并陈述答题理由。 1． 判断对错：两条序列的同源程度为60%。

2． 判断对错：两条序列的相似性很高，所以它们一定是同源序列。 3． 计算下面两条序列的海明距离：

AGCAAACACACTA ACATAAGCACACA

4． 通过字符编辑操作将序列s转换成t

s：AG-CA t：ACAC-

5． 分析两条序列的关系时，（）方法可以通过观察矩阵对角线迅速发现可能的序列比对。

A Dot-plot B Pairwise-Alignment C BLAST D FASTA E Score Matrix 6． 下列哪些是核酸序列数据库（）。

A GenBank B PDB C Entrez D EMBL-Bank E DDBJ

7． 下列哪些是蛋白质序列数据库（）。

A PIR B SWISS-PROT C TrEMBL D EPD E PDB

8． PDB文件的显示序列信息中，关键字（）作为显式序列标记，以该关键字打头的每一行都是关于序列的信息。

A HEADER B REMARK C SEQRES D EXPDTA 9． 下列哪些是可以用来显示分子结构的软件（）。

A GCG-DS Visualizer B RasMol C ChemView D DSSP

10． Entrez数据库集成系统中集成了NCBI中哪些数据库中的信息（）。

A 核酸序列 B 蛋白质序列 C 生物大分子结构 D 基因组数据 E 生物分类数据库 F 孟德尔人类遗传学数据（OMIM） G Pubmed

11． 排序：一个基本的DNA序列分析方案的分析顺序是：

1 发现重复元素 2分析功能位点 3数据库搜索 4综合分析 5 序列组成统计分析

12． 序列ATTCGATCGCAA的三种ORF顺序是：

第二部分：例题演示及讲解

形式：通过对例题的演示以及讲解，使同学能够更好的理解序列分析中数据的获得及处理方法，开发同学的分析问题解决问题的能力。

例题1、登陆到NCBI，查看Entrez包括哪些数据库，了解Entrez系统的使用方法。

附：查询OMIM向导实例：

1、What human genes are related to hypertension? Which of those genes are on chromosome 17?

2、List the OMIM entries that describe genes on chr10.

3、List the OMIM entries that contain information about allelic variants.

4、 Retrieve the OMIM record for the cystic fibrosis transmembrane conductance regulator (CFTR), and link to related protein sequence records in Entrez.

CFTR: is an ABC transporter-class protein and ion channel that transports chloride ions across epithelial cell membranes.

5、 Find the OMIM record for the p53 tumor protein, and linkout to related information in Entrez Gene and the p53 Mutation Database.

例题2、熟悉EBI的CLUSTALW程序， http://www.ebi.ac.uk/clustalw/

将下列三种蝎毒蛋白用CLUSTALW进行比对，将返回结果写入Word文档，并说明三种蛋白中哪两种的同源性更强。 1. 1SN1:A 2. 1SN4:A 3. 1SNB:_

从PDB中下载序列的FASTA格式，然后提交到CLUSTALW 例题3、序列分析

在真核生物启动子数据库中下载启动子数据：

TSS上游-250bp到TSS下游49bp总长度为300bp的人类启动子数据 手工剔出含N的序列

自己编写程序，计算出ATCG4碱基频率表以及频率矩阵。 用相同的操作对已知的非启动子数据集进行分析。

在Excel或其它软件中，对启动子和非启动子的频率表做统计图，发现ATCG4个碱基在不同位置出现频率的差异。

第三部分：习题作业

形式：学生综合掌握例题中的知识点及方法后，仿照例题完成以下习题。

习题1、任选一个数据库，熟悉其查询和比对的操作，将操作方法和查询比对的例子写入word文档。

提示：阅读数据库首页的说明以及帮助。

习题2、试从Genbank中查找10条SARS冠状病毒（coronovirus）的基因组序列，利用软件clustalw进行多重比对。 Clustalw for windows

ftp://ftp.ebi.ac.uk/pub/software/clustalw2/

习题2、试通过生物大分子结构数据库PDB服务器取出一个蛋白质结构数据，并用相应的分子图形学软件显示该结构。 GCG-setupdsv201.exe PDB Rasmol Chemview

习题3、编写程序，对例3 中的序列进行处理。统计出相邻两个碱基的频率表，并对序列进行替换，将原来的字符序列矩阵转换成频率矩阵。

同样在Excel中绘图，观察16种碱基的两两组和在300个位置上出现的频率特征。