无参考基因的转录组分析

基因注释到 GO 条目结果文件示例
GO 条目与 All-Unigene 对应结果文件示例
5、Unigene 代谢通路分析
KEGG 是系统分析基因产物在细胞中的代谢途径以及这些基因产物的功能的数据库，利 用 KEGG 可以进一步研究基因在生物学上的复杂行为。 根据 KEGG 注释信息我们能进一步得到 Unigene 的 Pathway 注释。
表 1 测序错误率与测序质量值简明对应关系
测序错误率
测序质量值
对应字符
5%
13
M
1%
20
T
0.1%
30
^
去除杂质数据
某些原始序列带有 adaptor 序列， 或含有少量低质量序列。 我们首先经过一系列数据处理以 去除杂质数据，得到 Clean reads。 数据处理的步骤： 1. 2. 3. 4. 去除含 adaptor 的 reads 去除 N 的比例大于 10%的 reads 去除低质量 reads（质量值 Q≤5 的碱基数占整个 read 的 50％以上） 获得 Clean reads
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;
注释到代谢通路结果文件示例
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'）和氨基酸序列。'
Clean Reads 数据
原始序列数据经过去除杂质后得到的数据。产量统计和后续信息分析分析都基于 Clean Reads。
测序产量统计表格示例
Samples
Total Reads
Total Nucleotides (nt)
Q20 percentage
N percentage
GC percentage *
无参考基因的转录组分析 一、实验流程
提取样品总 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 进行测序。
每个序列共有 4 行，第 1 行和第 3 行是序列名称（有的 fq 文件为了节省存储空间会省 略第三行“＋”后面的序列名称），由测序仪产生；第 2 行是序列；第 4 行是序列的测序质 量，每个字符对应第 2 行每个碱基，第四行每个字符对应的 ASCII 值减去 64，即为该碱基 的测序质量值，比如 c 对应的 ASCII 值为 99，那么其对应的碱基质量值是 35。从 Illumina GA Pipeline v1.3 开始（目前为 v1.6），碱基质量值范围为 2 到 35。表 1 为测序错误率与 测序质量值简明对应关系。具体地，如果测序错误率用 E 表示，碱基质量值用 sQ 表示，则 有下列关系： sQ = -10lgE
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.件 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 我们给出组装软件得到的序列。
[5]
组装出来的序列长度是组装质量的一个评估标准。我们会对组装出来的 Contig、Scaffold、Unigene 做一个长度分布 统计。如下图所示，给出的 bar 图统计 Contig 的长度分布。横坐标是组装出来的 Contig 的长度， 纵坐标是
对
应
长
度
的
Contig
的
数
目
。
组装成功的 Contig 结果在文件夹 1.Contig， Scaffold 相关结果在文件夹 2.Scaffold， Unigene 相关结果在文 件夹 3.Unigene。文件的详细意义可见各个文件夹下面对应的 readme。注：文件夹下面 svg 图可能需要安装 svg 插件才能打开
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!
二、信息分析流程
1、产量统计
原始序列数据
测序得到的原始图像数据经 base calling 转化为序列数据，我们称之为 raw data 或 raw reads，结果以 fastq 文件格式存储，fastq 文件为用户得到的最原始文件，里面存储 reads 的序列以及 reads 的测序质量。在 fastq 格式文件中每个 read 由四行描述： \@FC61FL8AAXX:1:17:1012:19200#GCCAAT/1 CCACTGTCATGTGAACATCACAGAGACATTTCTTGA + bbbbbbbbbbbbbbbbbbbbbbbbbaaaaaaaaa_\\
然后， 我们对差异检验的 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 分析。