无参考基因的转录组分析

注释到代谢通路结果文件示例
6、预测编码蛋白框（CDS）
首先，我们按 nr、Swiss-Prot、KEGG 和 COG 的优先级顺序将 Unigene 序列与以上蛋白 库做 blastx 比对（evalue<0.00001），如果某个 Unigene 序列比对上高优先级数据库中的 蛋白，则不进入下一轮比对，否则自动跟下一个库做比对，如此循环直到跟所有蛋白库比对 完。我们取 blast 比对结果中 rank 最高的蛋白确定该 Unigene 的编码区序列，然后根据标 准密码子表将编码区序列翻译成氨基酸序列， 从而得到该 Unigene 编码区的核酸序列 （序列 方向 5'->3'）和氨基酸序列。 最后， 跟以上蛋白库皆比对不上的 Unigene 我们用软件 ESTScan （Iseli, Jongeneel et al. 1999）预测其编码区，得到其编码区的核酸序列（序列方向 5'->3'）和氨基酸序列。'
[5]
组装出来的序列长度是组装质量的一个评估标准。我们会对组装出来的 Contig、Scaffold、Unigene 做一个长度分布 统计。如下图所示，给出的 bar 图统计 Contig 的长度分布。横坐标是组装出来的 Contig 的长度， 纵坐标是
对
应
长
度
的
Contig
的
数
目
。
组装成功的 Contig 结果在文件夹 1.Contig， Scaffold 相关结果在文件夹 2.Scaffold， Unigene 相关结果在文 件夹 3.Unigene。文件的详细意义可见各个文件夹下面对应的 readme。注：文件夹下面 svg 图可能需要安装 svg 插件才能打开
3、Unigene 功能注释
功能注释信息给出 Unigene 的蛋白功能注释、COG 功能注释。 首先，通过 blastx 将 Unigene 序列比对到蛋白数据库 nr、Swiss-Prot t、KEGG 和 COG （evalue<0.00001） ， 得到跟给定 Unigene 具有最高序列相似性的蛋白， 从而得到该 Unigene 的蛋白功能注释信息。 COG 是对基因产物进行直系同源分类的数据库，每个 COG 蛋白都被假定来自祖先蛋白， COG 数据库是基于细菌、藻类、真核生物具有完整基因组的编码蛋白、系统进化关系进行构 建的， 我们将 Unigene 和 COG 数据库进行比对， 预测 Unigene 可能的功能并对其做功能分类 统计，从宏观上认识该物种的基因功能分布特征。
已知， 样本一能比对到所有 Unigene 的总 reads 数为 N1 ， 样本二能比对到所有 Unigene 的总 reads 数为 N 2 ，基因 A 在样本一中对应的 reads 数为 x ，在样本二中对应的 reads 数 为 y ，则基因 A 在两样本中表达量相等的概率可由以下公式计算 ：Clean Rea Nhomakorabeas 数据
原始序列数据经过去除杂质后得到的数据。产量统计和后续信息分析分析都基于 Clean Reads。
测序产量统计表格示例
Samples
Total Reads
Total Nucleotides (nt)
Q20 percentage
N percentage
GC percentage *
Sample_A
1,634,670
122,600,250
89.47%
0.00%
48.50%
* Total Nucleotides = Total Reads1 x Read1 size + Total Reads2 x Read2 size; Nucleotides are actually clean reads and clean nucleotides;
y y
[1]
2 p(i | x) （当 p(i | x) 0.5 时）
i 0 i 0
或者
y y
2(1 p(i | x)) （当 p(i | x) 0.5 时）
i 0 i 0
其中
p(i | x) (
N2 i ) N1
( x i )! N x! y!(1 2 ) ( x i 1) N1
[1]
H 0 : 某一个基因在两个样本中表达量相同 H a : 某一个基因在两个样本中表达量不同
假设观测到基因 A 对应基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布 ：
p( x)
e x （ 为基因 A 的真实转录数） x!
然后， 我们对差异检验的 p-value 作多重假设检验校正， 通过控制 FDR （False Discovery Rate）来决定 p value 的域值。假设挑选了 R 个差异表达基因，其中 S 个是真正有差异表达 的基因，另外 V 个是其实没有差异表达的基因，为假阳性结果。希望错误比例 Q＝V/R 平均 而言不能超过某个可以容忍的值（比如 1％），则在统计时预先设定 FDR 不能超过 0.01 （Benjamini, Yekutieli. 2001）。 在得到差异检验的 FDR 值同时，我们也根据基因的表达量（RPKM 值）计算该基因在不 同样本间的差异表达倍数。FDR 值越小，差异倍数越大，则表明表达差异越显著。在我们的 分析中，差异表达基因定义为 FDR≤0.001 且倍数差异在 2 倍以上的基因。 得到差异表达基因之后，我们对差异表达基因做 GO 功能分析和 KEGG Pathway 分析。
基因注释到 GO 条目结果文件示例
GO 条目与 All-Unigene 对应结果文件示例
5、Unigene 代谢通路分析
KEGG 是系统分析基因产物在细胞中的代谢途径以及这些基因产物的功能的数据库，利 用 KEGG 可以进一步研究基因在生物学上的复杂行为。 根据 KEGG 注释信息我们能进一步得到 Unigene 的 Pathway 注释。
Unigene 的 COG 功能注释结果样式示例
4、Unigene 的 GO 分类
根据 nr 注释信息我们能得到 GO 功能注释。Gene Ontology（简称 GO）是一个国际标准 化的基因功能分类体系，提供了一套动态更新的标准词汇表（controlled vocabulary）来 全面描述生物体中基因和基因产物的属性。GO 总共有三个 ontology，分别描述基因的分子 功能（molecular function）、所处的细胞位置（cellular component）、参与的生物过程 （biological process） 。 我们根据 nr 注释信息， 使用 Blast2GO 软件 （Conesa, Gotz et al. 2005）得到 Unigene 的 GO 注释信息。Blast2GO 已被其它文献引用超过 150 次，是同行广泛 认可的 GO 注释软件。得到每个 Unigene 的 GO 注释后，我们用 WEGO 软件（Ye, Fang et al. 2006）对所有 Unigene 做 GO 功能分类统计，从宏观上认识该物种的基因功能分布特征。
无参考基因的转录组分析 一、实验流程
提取样品总 RNA 后，用带有 Oligo（dT）的磁珠富集真核生物 mRNA（若为原核生物， 则用试剂盒去除 rRNA 后进入下一步） 。 加入 fragmentation buffer 将 mRNA 打断成短片段， 以 mRNA 为模板，用六碱基随机引物（random hexamers）合成第一条 cDNA 链，然后加入缓 冲液、dNTPs、RNase H 和 DNA polymerase I 合成第二条 cDNA 链，在经过 QiaQuick PCR 试 剂盒纯化并加 EB 缓冲液洗脱之后做末端修复、加 poly（A）并连接测序接头，然后用琼脂 糖凝胶电泳进行片段大小选择00 进行测序。
7、Unigene 表达差异分析
差异表达分析找出在不同样本间存在差异表达的基因， 并对差异表达基因做 GO 功能分 析和 KEGG Pathway 分析。
差异表达基因的筛选
参照 Audic S.等人发表在 Genome Research 上的基于测序的差异基因检测方法 ，我 们开发了严格的算法筛选两样本间的差异表达基因。 筛选差异表达基因所用的到的假设检验的零假设和备择假设如下：
Total Reads and Total
Q20 percentage is proportion of
nucleotides with quality value larger than 20; N percentage is proportion of unknown nucleotides in clean reads; GC percentage is proportion of guanidine and cytosine nucleotides among total nucleotides.
2、组装结果
我们使用短 reads 组装软件 SOAPdenovo 做转录组从头组装。SOAPdenovo 首先将具有 一定长度 overlap 的 reads 连成更长的片段，这些通过 reads overlap 关系得到的不含 N 的组装片段我们称之称为 Contig。 然后， 我们将 reads 比对回 Contig， 通过 paired-end reads 能确定来自同一转录本的不同 Contig 以及这些 Contig 之间的距离，SOAPdenovo 将这些 Contig 连在一起，中间未知序列用 N 表示，这样就得到 Scaffold。进一步利用 paired-end reads 对 Scaffold 做补洞处理，最后得到含 N 最少，两端不能再延长的序列，我们称之为 Unigene。如果同一物种做了多个样品测序，则不同样品组装得到的 Unigene 可通过序列聚 类软件做进一步序列拼接和去冗余处理，得到尽可能长的非冗余 Unigene。 最后，将 Unigene 序列与蛋白数据库 nr、Swiss-Prot、KEGG 和 COG 做 blastx 比对 （evalue<0.00001），取比对结果最好的蛋白确定 Unigene 的序列方向。如果不同库之间的 比对结果有矛盾，则按 nr、Swiss-Prot、KEGG 和 COG 的优先级确定 Unigene 的序列方向， [3] 跟以上四个库皆比不上的 Unigene 我们用软件 ESTScan 预测其编码区并确定序列的方向。 对于能确定序列方向的 Unigene 我们给出其从 5'到 3'方向的序列，对于无法确定序列方向 的 Unigene 我们给出组装软件得到的序列。