转录组有参考生物信息分析结题报告模版-V2.0

北京诺禾致源生物信息科技有限公司

转录组有参考基因组生物信息分析结题报告

获得原始测序序列（Sequenced Reads）后，并且其相应的基因组参考序列( Reference Genome )

可以获得的情况下，可以用有参考基因组信息分析流程对数据进行详细的分析，分析流程图如下：

1

北京诺禾致源生物信息科技有限公司

1. 原始序列数据

高通量测序（如Illunima HiSeq

TM

2000/ Miseq等测序平台）测序得到的原始图像数据文件经碱基

识别（Base Calling）分析转化为原始测序序列（Sequenced Reads），我们称之为Raw Data或Raw Reads，

结果以FASTQ（简称为fq）文件格式存储，其中包含测序序列（reads）的序列信息以及其对应的测

序质量信息。测序样品中真实数据随机截取结果如下：

@HWI-ST1106:227:D14F6ACXX:1:1101:1202:2188 1:N:0:GCCAAT

CGGATGATCTTCTTAATCTCTCCTTGCATAGTTATGAAACAGTCCGTGGACTTGCTGGAAAATCTCTCTTGAAGATGATGAAGAGATGGCCCTCTACAAT

+

CCCFFFDFFHHHHJJJJJIJIGGGIGICIGIIJEIIJIIJJI@DHEDHECFGGAHGGJGHIICGEEIEHGGGIECEEHH@HE>C@EBBE@CCDDCCCDDC

@HWI-ST1106:227:D14F6ACXX:1:1101:1237:2217 1:N:0:GCCAAT

GAAGGTGAGTCTGAGGAGGCCAAGGAGGGAATGTTTGTGAAAGGATATGTCTACTAAGATATTAGAAAGTATGTACTACTACTACTACTACATGTTTTCA

+

@@@FDADDFDHFHIIIDHIIJJJGICGGGCGHGFIGHBHEHHGI;BDHHCFGCHIIIIEHGIGHHIJJE7??ACHCDFFFFFEEECCEE>C>ACCCDC>@

@HWI-ST1106:227:D14F6ACXX:1:1101:1382:2195 1:N:0:GCCAAT

TTTTGCAACAATGGCTTCCACCATGATGACTACTCTACCACAGTTCAATGGACTCAAACCCCAACCTTTCTCAGCTTCTCCAATTCAAGGCTTGGTGGCA

+

@@@DD3DDFFFF:CDGI@GIEEDH@AAD;;;@@#######

@HWI-ST1106:227:D14F6ACXX:1:1101:1255:2239 1:N:0:GCCAAT

CGGATTTTCAAGGGCCGCCGGGAGCGCACCGGACACCACGCGACGTGCGGTGCTCTTCCAGCCGCTGGACCCTACCTCCGGCTGAGCCGATTCCAGGGTG

+

CCCDFFFFHHH?FHIIIJJJJJIGBEHHJJBHBDDCDAC??@@BDBBBBD8BDDCDDACC@A?@BBB@<(8>?395?4:(:<@##

@HWI-ST1106:227:D14F6ACXX:1:1101:1423:2239 1:N:0:GCCAAT

CTTGTATTGCTCTCCCACAACCCCGTTTTCACGGTTTAGGCTGCTCCCATTTCGCTCGCCGCTACTACGGGAATCGCTTTTGCTTTCTTTTCCTCTGGCT

+

CCCFDFFFHHHHHJJIJJJJJIJJGGIHIIGIIJGIGGIJJGGGJGIJ>FGIIGHGGBEHBCCBBDDD@BB@@@C::@C

2

北京诺禾致源生物信息科技有限公司

2.

测序数据质量评估

2.1 测序错误率分布检查

测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。通常测序

序列（Sequenced Reads）5’端前几个碱基的错误率相对较高，随着序列的延伸，3’端碱基错误率会

不断升高，这是由高通量测序的技术特点决定的。项目结果见图1。

图1 测序错误率分布图

横坐标为reads的碱基位置，纵坐标为单碱基错误率。其中前100个碱基位置为双端

测序序列的第一端测序Reads的分布情况，随后100bp是另一端测序reads的分布情况。

3

北京诺禾致源生物信息科技有限公司

2.2 A/T/G/C含量分布检查

对于RNA-seq来说，因随机性打断及G/C和A/T含量分别相等的原则，理论上GC及AT含量每

个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。项目结果见图2。

图2 GC含量分布图

横坐标为reads的碱基位置，纵坐标为单碱基所占的比例。不同颜色代表不同的碱基类型

4

北京诺禾致源生物信息科技有限公司

2.3 测序数据过滤

测序得到的原始测序序列（Sequenced Reads）或者raw reads，里面含有带接头的、低质量的

reads，为了保证信息分析质量，必须对raw reads过滤，得到clean reads，后续分析都基于clean reads。

项目结果见图3。

图３ 原始数据过滤结果

5

北京诺禾致源生物信息科技有限公司

2.4 测序数据质量情况汇总

样品测序产出数据质量评估情况详见表1。

表1 数据产出质量情况一览表

Sample

TS-1_1

TS-1_2

TS-2_1

TS-2_2

TR-3_1

TR-3_2

TR-5_1

TR-5_2

Raw reads

48900437

48900437

50753113

50753113

37877095

37877095

55854530

55854530

Clean reads

48857403

48857403

50709069

50709069

37819080

37819080

55791168

55791168

Bases

4.89G

4.89G

5.07G

5.07G

3.78G

3.78G

5.58G

5.58G

Error(%)

0.03

0.03

0.03

0.03

0.03

0.04

0.03

0.03

Q20(%)

98.24

96.96

98.26

97.11

97.41

95.91

97.9

96.66

Q30(%)

94.22

91.59

94.3

91.92

92.21

89.18

93.38

90.83

GC(%)

43.69

43.68

43.41

43.44

47.79

47.81

45.57

45.54

Dup(%)

79.61

78.21

79.26

77.83

82.78

81.45

81.97

80.75

Raw Reads：由测序得到的原始图像数据经base calling 转化而来的原始序列reads。

Clean reads：将Raw Reads过滤得到的reads。

Bases (Clean bases)： 过滤得到的数据的总碱基数。

Error (Error rate)：指测序错误率，与碱基质量值之间有一定的对应关系。

Q20：测序错误率≤1%的碱基数目比例。

Q30：测序错误率≤0.1%的碱基数目比例。

GC content：G+C的数量占总的碱基数量的百分比。

Dup (Duplication level)：重复的reads数占总reads数的比例。

6

北京诺禾致源生物信息科技有限公司

3.

参考序列比对分析

3.1 RNA-Seq reads参考基因组比对统计

如果参考基因组选择合适并且相关实验不存在污染的情况下，实验所产生的测序序列的定位的

百分比正常情况下会高于70% (Total Mapped Reads or Fragments)，其中具有多个定位的测序序列

(Multiple Mapped Reads or Fragments)占总体的百分比通常不会超过10%。项目结果见表2。

表2 参考基因组比对的统计情况一览表

Sample name

Total reads

Total mapped

Multiple mapped

Uniquely mapped

Read-1

Read-2

Reads map to '+'

Reads map to '-'

Non-splice reads

Splice reads

Reads mapped in

proper pairs

TS1

97714806

88921431(91%)

961182 (0.98%)

87960249 (90.02%)

44157413 (45.19%)

43802836 (44.83%)

43944185 (44.97%)

44016064 (45.05%)

64211264 (65.71%)

23748985 (24.3%)

77892308 (79.71%)

TS2

101418138

92225043 (90.94%)

1053580 (1.04%)

91171463 (89.9%)

45745753 (45.11%)

45425710 (44.79%)

45558208 (44.92%)

45613255 (44.98%)

67345511 (66.4%)

23825952 (23.49%)

82296934 (81.15%)

TR3

75638160

54554812 (72.13%)

1422941 (1.88%)

53131871 (70.24%)

26677577 (35.27%)

26454294 (34.97%)

26479698 (35.01%)

26652173 (35.24%)

42758444 (56.53%)

10373427 (13.71%)

45364262 (59.98%)

TR5

111582336

91773760 (82.25%)

1481392 (1.33%)

90292368 (80.92%)

45296200 (40.59%)

44996168 (40.33%)

45055089 (40.38%)

45237279 (40.54%)

68211989 (61.13%)

22080379 (19.79%)

80126002 (71.81%)

7

北京诺禾致源生物信息科技有限公司

3.2 RNA-Seq reads参考基因组比对分布图

定位到基因组上的测序序列分布统计，用于检测测序序列基因组上的来源。项目结果见图4。

图４ RNA-Seq 测序得到的reads比对到参考基因组不同区域上的分布情况

8

北京诺禾致源生物信息科技有限公司

3.3 RNA-Seq reads参考序列密度分布图

对定位到基因组上的测序序列完成染色体密度分布统计，用于检测染色体上测序序列分布的异

常情况。项目结果见图5。

图５ RNA-Seq 测序得到的reads比对到参考基因组不同染色体上的分布情况

横坐标为染色体的长度信息（以百万碱基为单位），纵坐标为log2（reads的密度的中位数）

9

北京诺禾致源生物信息科技有限公司

4. 可变剪切分析

对该物种及其相应的测序样品进行可变剪切事件的统计。项目结果见图6。

图６ 可变剪切类型分析

横坐标为可变剪切事件的五种分类缩写，纵坐标为该种事件下可变剪切的数量，

不同颜色代表不同的样品组合或者已知的基因模型

（1）

（2）

（3）

（4）

（5）

Skipped exon （SE）;外显子跳跃

Retained intron（RI）; 内含子滞留

Alternative 5’ splicing stie（A5SS）; 可变5’端剪切

Alternative 3’ splicing site（A3SS）; 可变3’端剪切

Mutually exclusive exon （MEX）; 互相排斥的外显子

10

北京诺禾致源生物信息科技有限公司

5．SNP分析

图7

11

北京诺禾致源生物信息科技有限公司

6. 新转录本预测

对所分析的物种在已知的基因模型的基础上，用所有测序的数据对新转录区域进行预测，并对

新转录区域的表达水平进行统计分析，项目结果见图8,9,10。

图8 新转录本的RPKM累积分布图

12

北京诺禾致源生物信息科技有限公司

图9 新转录本的RPKM盒形图

13

北京诺禾致源生物信息科技有限公司

图10 新转录本的RPKM密度分布图

14

北京诺禾致源生物信息科技有限公司

7. 基因表达水平分析

在RNA-技术中，RPKM（Reads Per Kilo bases per Million mapped Reads）是一种表示基因表达水

平的通用方法，代表每百万reads中来自于某基因每千碱基长度的reads数。项目结果见表3，4，

表4详细内容见 ./Results/DEseq/。

表3 不同表达水平区间的基因数量统计表

RPKM Interval

0 - 0.01

0.1 - 0.3

0.3 - 3.57

3.57 - 15

15 - 60

> 60

Gene Counts

TS1

13724 (36.73%)

2032 (5.44%)

9783(26.18%)

7641 (20.45%)

3041 (8.14%)

1143 (3.06%)

Gene Counts

TS2

13459 (36.02%)

1909 (5.11%)

9786 (26.19%)

7981 (21.36%)

3116 (8.34%)

1113 (2.98%)

Gene Counts

TR3

14515 (38.85%)

2101 (5.62%)

9984 (26.72%)

7075 (18.94%)

2668 (7.14%)

1021 (2.73%)

Gene Counts

TR5

13728 (36.74%)

1982 (5.30%)

10379 (27.78%)

7424 (19.87%)

2800 (7.49%)

1051 (2.81%)

表4 基因表达水平统计表 （部分）

chromosome

SL2.40ch00

SL2.40ch00

SL2.40ch00

SL2.40ch00

SL2.40ch00

SL2.40ch00

gene_id

Novo_00001

Novo_00002

Novo_00003

Novo_00004

Novo_00005

Novo_00006

start

876120

1146195

1230237

4277120

4641283

4640465

end

876253

1147157

1234310

4277288

4642496

4641341

RPKM

(TS1)

0.965

0.282

0.387

17.526

1.332

0.546

RPKM

(TS2)

3.255

0.440

0.493

20.280

1.365

1.108

RPKM

(TR3)

0.985

0.000

0.142

22.794

0.000

0.000

RPKM

(TR5)

0.403

0.000

0.351

27.304

0.000

0.000

15

北京诺禾致源生物信息科技有限公司

8

. RNA-seq整体质量评估

8.1 均一性分布检查

根据转录组建库实验的特点，转录本其产生的测序序列（reads）实际覆盖度的分布特点见下图：

距离转录本的5'端和3'端越近，平均测序深度越低，但总体的均一化程度比较高。项目结果见图11。

图11 不同表达水平的转录本的reads密度分布图

High：高表达量转录本；Medium：中度表达量转录本；Low：低表达量转录本。

横坐标为距离转录本5’端的相对位置（以百分比表示），纵坐标为覆盖深度的平均值。

8.2 饱和曲线检查

定量饱和曲线检查反映了基因表达水平定量对数据量的要求。表达量高的基因，就越容易被准

确定量；反之，表达量低的基因，需要较大的测序数据量才能被准确定量。项目结果见图12。

图12 定量饱和曲线检查分布图

横坐标代表定位到基因组上的reads数占总reads数的百分比，纵坐标代表定量误差

在10%以内的基因的比例。n代表在此RPKM范围内的基因数。

16

北京诺禾致源生物信息科技有限公司

8.3 重复相关性检查

样品间基因表达水平相关性是检验实验可靠性和样本选择是合理性的重要指标。相关系数越接

近1，表明样品之间表达模式的相似度越高。通常生物学重复要求 R

2

> 0.92。

图13 RNA-Seq 生物学重复相关性分析

对取自同一时间点均进行RNA-Seq测序的样本进行生物学重复相关性分析

17

北京诺禾致源生物信息科技有限公司

9. 差异表达分析

9.1 不同试验条件下，基因表达水平对比图

图14 在不同试验条件，所有基因的表达水平的箱线图

18

北京诺禾致源生物信息科技有限公司

9.2 样本间差异基因筛选

图15 样本间基因差异表达分析火山图

有显著性差异表达的基因用红色点表示；横坐标代表基因在不同样本中表达倍数变化；

纵坐标代表基因表达量变化差异的统计学显著性，红色圆点表示有显著性差异的基因。

19

北京诺禾致源生物信息科技有限公司

9.3 差异表达基因列表

项目结果见 ./Results/DEseq/。

表5 差异基因列表（部分）

chromosome

SL2.40ch06

SL2.40ch10

SL2.40ch06

SL2.40ch11

SL2.40ch06

SL2.40ch04

gene_id

Solyc06g009960

Novo_02114

Novo_01304

Solyc11g028010

Novo_01303

Novo_00935

RPKM1

0.107986

0.04235

0.015863

19.92548

0.006346

0.017206

RPKM2

158.3342

63.33414

22.18282

0.021828

31.31789

41.81995

log2FoldChange

10.7241

10.77632

10.67336

-9.55392

12.53154

11.47674

padj

2.46E-221

3.06E-200

4.24E-126

2.38E-125

1.50E-121

3.76E-119

20

北京诺禾致源生物信息科技有限公司

9.4 差异基因维恩图

维恩图展示了样品两两间差异基因的个数，以及共同差异基因的个数

图16 多样品差异基因维恩图

不同样品相互比较，找出差异基因（仅针对多个样品）

21

北京诺禾致源生物信息科技有限公司

9.5 差异基因聚类分析

差异基因聚类分析用于判断不同实验条件下调控模式的聚类模式，不同样品比较组合之间求

RPKM的差异倍数，用于聚类分析。项目结果见图17。

图17 差异基因聚类图

（红色表示上调，绿色表示下调）

22

北京诺禾致源生物信息科技有限公司

10. 差异基因Gene Ontology富集分析

10.1 差异表达基因GO富集分析

项目结果见 ./Results/GOseq/。

表6 样品中差异基因的Gene Ontology 分类（部分）

GO ID

GO:0005813

GO:0005634

GO:0005351

GO:0042645

GO:0043533

GO:0005737

GO:0032792

GO:0005896

GO:0016853

GO:0005515

GO Term

centrosome

nucleus

sugar:hydrogen symporter activity

mitochondrial nucleoid

inositol 1,3,4,5 tetra kisphosphate

binding

cytoplasm

negative regulation of CREB

transcription factor activity

interleukin-6 receptor complex

isomerase activity

protein binding

Term

type

C

C

M

C

M

C

P

C

M

M

P-Value

6.50E-06

5.47E-05

0.000130

0.000156

Correct DEG

P-Value item

0.07377

0.31026

0.44389

0.44389

17

108

3

5

Bg/ref

item

252

4569

6

31

2

4488

3

3

111

8310

0.000312 0.67799 2

0.000358 0.67799 105

0.000542 0.87928 2

0.000626 0.88749 2

0.000714 0.90057 7

0.000803 0.91167 172

上述表格为差异基因的Gene Ontology富集分析结果表格。

GO ID: Gene Ontology数据库中唯一的标号信息

GO Term: Gene Ontology功能的描述信息

Term type为该GO的类别（C:细胞组分；P：生物路径；M：分子功能）

P-value: 富集分析统计学显著水平，一般情况下，P-value < 0.05 该功能为富集项

Correct P-Value 矫正后的P-Value

DEG item：差异基因中与该Term相关的基因数

DEG item：所有（bg）基因中与该Term相关的基因数

23

北京诺禾致源生物信息科技有限公司

10.2 GO有向无环图分析

GO有向无环图能直观展示差异基因富集的GO term及其层级关系。项目结果见图15。

图18 差异基因在生物学功能方面的GO富集有向无环图

每个节点代表一个GO术语，方框代表的是富集程度为TOP10的GO，颜色的深浅代表富集程度，

颜色越深就表示富集程度越高，每个节点上展示了该TERM的名称及富集分析的P-value

24

北京诺禾致源生物信息科技有限公司

11. 差异基因KEGG富集分析

11.1 差异基因KEGG显著性富集列表

项目结果见 ./Results/GOseq/。

表7 样品中差异基因的KEGG富集（部分）

Term Database Id Sample

number

Photosynthesis

Inositol phosphate metabolism

Tyrosine metabolism

Indole alkaloid biosynthesis

RNA polymerase

Linoleic acid metabolism

KEGG PATHWAY

KEGG PATHWAY

KEGG PATHWAY

KEGG PATHWAY

KEGG PATHWAY

KEGG PATHWAY

ath00195

ath00562

ath00350

ath00901

ath03020

ath00591

5

3

2

1

2

1

Background

number

77

58

30

7

35

9

0.011839

0.084135

0.099162

0.118239

0.128205

0.149395

P-Value Corrected

P-Value

1

1

1

1

1

1

25

北京诺禾致源生物信息科技有限公司

11.2 差异基因KEGG富集散点图

差异基因KEGG富集分析的散点图，其中四个象限内的KEGG通路的富集程度不同。

图19 差异基因代谢通路分析的四象限图

Rich factor 是差异表达的基因中位于该pathway 条目的基因数目与所有有注释基因中位于该pathway 条目的基因数的比值。Rich factor越大，表

示富集的程度越大。-log10（Qvalue）中Qvalue是做过多重假设检验校正之后的pvalue，-log10（Qvalue）越大，表示富集越显著。

26

北京诺禾致源生物信息科技有限公司

12．差异基因蛋白互作网络分析

差异基因的蛋白互相作用的信息，结果如图20所示：

图20 差异基因相互作用网络分析

27

北京诺禾致源生物信息科技有限公司

1. 原始序列数据

FASTQ格式文件中每个read由四行描述，如下：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG

GCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCTGT

+

@@CFFFDEHHHHFIJJJ@FHGIIIEHIIJBHHHIJJEGIIJJIGHIGHCCF

其中第一行以“@”开头，随后为Illumina 测序标识别符（Sequence Identifiers）和描述文字（选

择性部分）。第二行是碱基序列；第三行以“+”开头，随后为Illumina 测序标识别符（选择性部分）；

第四行是对应序列的测序质量；(Cock et al.)

Illumina 测序标识别符（Sequence Identifiers）详细信息如下：

EAS139

136

FC706VJ

2

2104

15343

197393

1

Y

18

ATCACG

Unique instrument name

Run ID

Flowcell ID

Flowcell lane

Tile number within the flowcell lane

'x'-coordinate of the cluster within the tile

'y'-coordinate of the cluster within the tile

Member of a pair, 1 or 2 (paired-end or mate-pair reads only)

Y if the read fails filter (read is bad), N otherwise

0 when none of the control bits are on, otherwise it is an even number

Index sequence

第四行中每个字符对应的ASCII值减去33，即为对应第二行碱基的测序质量值。如果测序错误

率用e表示，Illunima HiSeq

TM

2000/ Miseq的碱基质量值用Q

phred

表示，则有下列关系：

公式1:

Illunima Casava 1.8版本测序错误率与测序质量值简明对应关系

测序错误率

5%

1%

0.1%

0.01%

测序质量值

13

20

30

40

对应字符

.

5

?

I

碱基识别（Base Calling）分析软件：Illunima Casava 1.8版本

测序参数：双端测序（Paired end）；测序序列读长：100bp （或者单位为循环数(cycle)）

28

北京诺禾致源生物信息科技有限公司

2. 测序数据质量评估

2.1 测序错误率分布检查

每个碱基测序错误率是通过测序Phred数值（Phred score ,Q

phred

）通过公式1转化得到，而Phred

数值是在碱基识别(Base Calling)过程通过一种预测碱基判别发生错误概率模型计算得到的，对应关

系如下表所显示：

Illunima Casava 1.8版本碱基识别与Phred分值之间的简明对应关系

Phred分值

10

20

30

40

不正确的碱基识别

1 / 10

1 / 100

1 / 1000

1 / 10000

碱基正确识别率

90%

99%

99.9%

99.99%

Q-sorce

Q10

Q20

Q30

Q40

对于RNA-seq技术，测序错误率分布具有两个特点：

（1）测序错误率会随着测序序列（Sequenced Reads）的长度的增加而升高，这是由于测序过

程中化学试剂的消耗而导致的，并且为illumina高通量测序平台都具有的特征（Erlich and Mitra, 2008;

Jiang et al.）。

（2）前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq建库

过程中反转录所需要的随机引物的长度。所以推测这部分较高的测序错误率原因为随机引物和RNA

测序错误率分布检查用于检测在测序长度范围内，有无异常的碱基位置存在高错误率，比如中间位

置的碱基的测序错误率显著的高于其他位置。一般情况下，每个碱基位置的测序错误率应该低于

0.5 %。

2.2 A/T/G/C含量分布检查

在illumina测序平台的转录组测序中，反转录成cDNA时所用的6bp 的随机引物会引起前几个

位置的核苷酸的组成存在一定的偏好性。而这种偏好型独立于所测序的物种和实验室，并会影响转

29

北京诺禾致源生物信息科技有限公司

录组测序的均一化程度(Hansen et al.)。

除此之外，理论上G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定

不变，呈水平线。

A/T/G/C含量分布检查用于检测有无AT、GC 分离现象, 而这种现象可能是测序或者建库所带来

的，并且会影响后续的定量分析。

2.3 测序数据过滤

数据处理的步骤如下：

1） 去除带接头（adopter）的reads。

2） 去除N（N表示无法确定碱基信息）的比例大于10%的reads。

3） 去除低质量reads（质量值sQ <= 5的碱基数占整个read的50％以上的reads）。

RNA-seq 的接头（Adopter, Oligonucleotide sequences for TruSeq

TM

Small RNA Sample Prep Kit）信息：

RNA 5’ Adapter (RA5), part # 15013205: 5’ -GUUCAGAGUUCUACAGUCCGACGAUC-3’

RNA 3’ Adapter (RA3), part # 15013207: 5’- TGGAATTCTCGGGTGCCAAGG-3’

原始测序序列（Sequenced Reads）的数据处理（过滤）结果：

不同颜色代表对不同过滤条件的统计

30

北京诺禾致源生物信息科技有限公司

2.4 测序数据质量情况汇总

Raw reads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列（Sequenced

Reads）的个数。

Clean reads: 计算方法同 Raw Reads、Raw bases，只是统计的文件为过滤后的测序数据。后续

的生物信息分析都是基于Clean reads。

Clean bases： 测序序列（Sequenced Reads）的个数乘以测序序列的长度，并转化为以G为单

位。Clean bases 要大于合同上规定的测序数据量。

Error rate：通过公式1计算得到。

Q20、Q30：分别计算 Phred 数值大于 20、30的碱基占总体碱基的百分比。

GC content：计算碱基G和C的数量总和占总的碱基数量的百分比。

Sequence duplication level：重复的reads数占总reads数的比例，重复的reads的定义为两个reads

在序列水平上完全相同。相对于基因组测序，RNA-seq的较高的测序序列重复水平是由于PCR 扩增

带来的(Mamanova et al.)。除了PCR-free实验（例如FRT-seq）产生的RNA-seq数据，reads重复率一

般都会在 70% 左右。而简单的去除这部分测序重复的RNA-seq数据并不会提高定量的准确性，所

以现有基因表达水平的分析方法是保留这部分数据的基础上进行的，而在做 SNP 分析的时候需要

除去这部分数据(Chepelev et al., 2009)。

3.

参考序列比对分析

3.1 RNA-Seq reads参考基因组比对统计

测序序列定位算法：根据不同的基因组的特征，我们选取相对合适的软件（动植物用TopHat

(Trapnell et al., 2009)、真菌或者基因密度较高的物种用BLAT (Kent, 2002)），合适参数设置（如最大

的内含子长度，会根据已知的该物种的基因模型来进行统计分析），将过滤后的测序序列进行基因

组定位分析，下图为TOPHAT 的算法示意图：

31

北京诺禾致源生物信息科技有限公司

Tophat的算法主要分为两个部分：

（1） 将测序序列整段比对到外显子上。

（2） 将测序序列分段比对到两个外显子上。

比对结果统计表：

Sample name

Total reads

Total mapped

Multiple mapped

Uniquely mapped

Read-1

Read-2

Reads map to '+'

Reads map to '-'

Non-splice reads

Splice reads

Reads mapped in

proper pairs

TS1

97714806

88921431(91%)

961182 (0.98%)

87960249 (90.02%)

44157413 (45.19%)

43802836 (44.83%)

43944185 (44.97%)

44016064 (45.05%)

64211264 (65.71%)

23748985 (24.3%)

77892308 (79.71%)

TS2

101418138

92225043 (90.94%)

1053580 (1.04%)

91171463 (89.9%)

45745753 (45.11%)

45425710 (44.79%)

45558208 (44.92%)

45613255 (44.98%)

67345511 (66.4%)

23825952 (23.49%)

82296934 (81.15%)

TR3

75638160

54554812 (72.13%)

1422941 (1.88%)

53131871 (70.24%)

26677577 (35.27%)

26454294 (34.97%)

26479698 (35.01%)

26652173 (35.24%)

42758444 (56.53%)

10373427 (13.71%)

45364262 (59.98%)

TR5

111582336

91773760 (82.25%)

1481392 (1.33%)

90292368 (80.92%)

45296200 (40.59%)

44996168 (40.33%)

45055089 (40.38%)

45237279 (40.54%)

68211989 (61.13%)

22080379 (19.79%)

80126002 (71.81%)

（1）Total reads: 测序序列经过测序数据过滤后的数量统计（Clean data）。

（2）Total mapped: 能定位到基因组上的测序序列的数量的统计; 一般情况下，如果不存在污

染并且参考基因组选择合适的情况下，这部分数据的百分比大于 70%。

（3）Multiple mapped: 在参考序列上有多个比对位置的测序序列的数量的统计; 这部分数据的

百分比一般会小于10%。

（4）Uniquely mapped: 在参考序列上有唯一比对位置的测序序列的数量的统计。

（5）Read-1, Read-2: 双端测序序列，其中两端分别能定位到基因组上的数量的统计; 两部分的

统计的比例应该大体相同。

（6）Reads map to '+', Reads map to '-': 测序序列比对到基因组上正链和负链的统计。

（7）Splice reads:（2）中，分段比对到两个外显子上的测序序列（也称为Junction reads）的统

计，Non-splice reads为整段比对到外显子的将测序序列的统计，Splice reads的百分比取决于测序片

段的长度。

（8）Reads mapped in proper pairs，双端测序序列定位到基因组上相对的距离符合测序片段的

长度分布。

3.2 RNA-Seq reads参考基因组比对分布图

对Total mapped的比对到基因组上的各个部分的情况进行统计，如下图所示，定位区域分为Exon

（外显子）, Intron（内含子）和 Intergenic（基因间隔区域）。

32

北京诺禾致源生物信息科技有限公司

正常情况下，Exon (外显子)区域的测序序列定位的百分比含量应该最高，定位到Intron (内含子)

区域的测序序列可能是由于非成熟的mRNA的污染或者基因组注释不完全导致的，而定位到

Intergenic（基因间隔区域）的测序序列可能为基因组注释不完全以及背景噪音。

3.3 RNA-Seq reads参考序列密度分布图

对Total mapped的比对到基因组上的各个染色体的密度进行统计，如下图所示，图中只保留前

30个最长的染色体，如果基因组中的染色体的数量不够30的情况则全部显示，具体作图的方法为用

滑动窗口（window size）为1K，计算窗口内部比对到碱基位置上的reads的中位数，并转化成 log

2

。

正常情况下，整个染色体长度越长，该染色体内部定位的reads总数会越多(Marquez et al.)。

33

北京诺禾致源生物信息科技有限公司

4. 可变剪切分析

用 Astalavista 软件 (Foissac and Sammeth, 2007; Sammeth, 2009; Sammeth et al., 2008) 对基因

组的已知的基因模型、Cufflinks (Trapnell et al.) 预测出的基因模型（包括单个实验组中所有重复样品

的数据, 所有实验样品的数据）, 具体算法流程以及结果图片如下所示：

可变剪切事件图例 (Wang et al., 2008)

Skipped exon （SE）;外显子跳跃

（1）

（2）

Retained intron（RI）; 内含子滞留

（3）

Alternative 5’ splicing stie（A5SS）; 可变5’端剪切

（4） Alternative 3’ splicing site（A3SS）; 可变3’端剪切

（5） Mutually exclusive exon （MEX）; 互相排斥的外显子

34

北京诺禾致源生物信息科技有限公司

5．SNP分析

具体缺相应的图表信息。

6. 新转录本预测

将所有数据的基因组定位结果放到一起，用 Cufflinks 进行组装，然后用 Cuffcompare 和 已知

的基因模型进行比较，发现新的转录区域。具体的算法流程和可变剪切分析里到 Astalavista 软件分

析之前是一样的。

35

北京诺禾致源生物信息科技有限公司

上图为不同区域的RPKM的累积分布（包括外显子，内含子，新转录本区域以及基因间隔区域），

从图上看，累积分布线上升的越快则代表相应区域的整体RPKM数值比其他上升较慢的区域低。（整

体的RPKM数值 外显子 > 新基因区域 > 内含子 > 基因间隔区域 ）。

上图为不同区域的RPKM的盒形图（包括外显子，内含子，新转录本区域以及基因间隔区域）。

每个区域的盒形图对五个统计量（至上而下分别为最大值，上四分位数，中值，下四分位数和最小

值）。

36

北京诺禾致源生物信息科技有限公司

上图为不同区域的RPKM的密度分布图（包括外显子，内含子，新转录本区域以及基因间隔区域）。

7. 基因表达水平分析

一个基因表达水平的直接体现就是其转录本的丰度情况，转录本丰度程度越高，则基因表达水

平越高。在RNA-seq分析中，我们可以通过定位到基因组区域或基因外显子区的测序序列（reads）

的计数来估计基因的表达水平。Reads计数除了与基因的真实表达水平成正比外，还与基因的长度、

测序深度成正相关。为了使不同基因、不同实验间估计的基因表达水平具有可比性，人们引入了

RPKM的概念，RPKM（Reads Per Kilo bases per Million reads）是每百万reads中来自某一基因每千碱基

长度的reads数目。RPKM同时考虑了测序深度和基因长度对reads计数的影响，是目前最为常用的基

因表达水平估算方法 (Mortazavi et al., 2008)。

不同表达水平区间的基因数量统计表

RPKM Interval

0 – 0.01

0.1 - 0.3

0.3 - 3.57

3.57 - 15

15 - 60

> 60

Gene Counts

Sample1

13724 (36.73%)

2032 (5.44%)

9783(26.18%)

7641 (20.45%)

3041 (8.14%)

1143 (3.06%)

Gene Counts

Sample2

13459 (36.02%)

1909 (5.11%)

9786 (26.19%)

7981 (21.36%)

3116 (8.34%)

1113 (2.98%)

Gene Counts

Sample3

14515 (38.85%)

2101 (5.62%)

9984 (26.72%)

7075 (18.94%)

2668 (7.14%)

1021 (2.73%)

Gene Counts

Sample4

13728 (36.74%)

1982 (5.30%)

10379 (27.78%)

7424 (19.87%)

2800 (7.49%)

1051 (2.81%)

37

北京诺禾致源生物信息科技有限公司

基因表达水平统计表

chromosome

SL2.40ch00

SL2.40ch00

SL2.40ch00

SL2.40ch00

SL2.40ch00

SL2.40ch00

gene_id

Novo_00001

Novo_00002

Novo_00003

Novo_00004

Novo_00005

Novo_00006

start

876120

1146195

1230237

4277120

4641283

4640465

end

876253

1147157

1234310

4277288

4642496

4641341

RPKM

(TS1)

0.965

0.282

0.387

17.526

1.332

0.546

RPKM

(TS2)

3.255

0.440

0.493

20.280

1.365

1.108

RPKM

(TR3)

0.985

0.000

0.142

22.794

0.000

0.000

RPKM

(TR5)

0.403

0.000

0.351

27.304

0.000

0.000

上表分别统计了不同表达水平下基因的数量以及单个基因的表达水平。一般情况下，RPKM 数

值0.1 或者 0.01 作为判断基因是否表达的阈值, 不同的文献所采用的阈值不同。RPKM的数据值为

1代表一个细胞中含有一个RNA 分子(Mortazavi et al., 2008)；RPKM 在 0.1 到 3.75 之间可以认为

是低丰度表达水平的基因；RPKM 在 3.75 到 15 之间为中等丰富度表达水平基因；RPKM大于15

之间为高丰富度表达水平基因。

8. RNA-seq整体质量评估

8.1 均一性分布检查

理想条件下，对于RNA-seq技术来说，测序序列（reads）之间为独立抽样并且reads在所有的表

达的转录本的上的分布应该呈现均一化分布。然而很多研究表明，很多偏好型的因素都会影响这种

均一化的分布(Dohm et al., 2008)。例如，在RNA-seq建库过程中，片段破碎和RNA反转录的顺序不一

样会导致RNA-seq最终的数据呈现严重的3’偏好性。其他因素还包括转录区域的GC含量不同、随机引

物等等，并且生物体内从5’或者3’的降级过程同样会导致不均一性分布。

如上图显示，左边为正常的RNA-seq的均一性分布；右边存在严重的3’端的偏好性，而右图的数

据会对的表达水平的定量分析产生影响。

38

北京诺禾致源生物信息科技有限公司

8.2 表达水平的饱和曲线检查

表达水平的饱和曲线的具体算法描述如下：

分别对10%, 20%,30% … 90%的总体测序数据分别的单独进行基因定量分析，并把所有数据条件

下得到的基因的表达水平作为最终的数值。用每个百分比条件下求出的单个基因的RPKM数值和最终

对应基因的表达水平数值进行比较，如果差异小于15%, 则认为这个基因在这个条件下定量是准确

的。

如上图所示，RPKM数值大于3.75的条件下基因的在该数据量的情况下对于基因定量分析而言

是饱和的。

8.3 重复相关性检查

生物学重复是任何生物学实验所必须的，高通量测序技术也不例外(Hansen et al.)。生物学重复

主要有两个用途：一个是证明所涉及的生物学实验操作是可以重复的且变异不大，另一个为后续的

差异基因分析所需要的。

如上图所示，两个生物学重复样品间的RPKM数值大部分都落到对角线上，并且皮尔逊相关系

数的平方大于0.92（Encode计划建议），如果小于0.92，需要对样品做出合适的解释，否则需要重新

进行实验。

39

北京诺禾致源生物信息科技有限公司

9. 差异表达分析

差异表达分析找出在不同样本间存在差异表达的基因。

参照Simon Anders等人于2010年发表在Genome Biology上的基于测序序列数据的差异基因检

测方法（DEseq）进行筛选两样本间的差异表达基因，判断条件为 FDR < 0.05。(Anders and Huber)

9.1 不同实验条件下，基因表达水平对比图

通过所有基因的RPKM的分布图以及盒形图对不同实验条件下的基因表达水平两两进行比较。对

于同一实验条件下的重复样品，最终的RPKM为所有重复数据的平均值。

RPKM分布图（左图）的横坐标为log

10

RPKM, 纵坐标为基因的密度。

RPKM盒形图（左图）的横坐标为样品名称，纵坐标为log

10

RPKM，每个区域的盒形图对五个

统计量（至上而下分别为最大值，上四分位数，中值，下四分位数和最小值）。

实例图上两个不同实验的RPKM分布差异不显著。

9.2 样本间差异基因筛选

用火山图可以推断差异基因的整体分布情况，从差异倍数(Fold change)和显著水平(P-value)两

个水平进行评估，如下图所示：

40

北京诺禾致源生物信息科技有限公司

9.3 差异表达基因列表

差异基因列表

chromosome

SL2.40ch06

SL2.40ch10

SL2.40ch06

SL2.40ch11

SL2.40ch06

SL2.40ch04

gene_id

Solyc06g009960

Novo_02114

Novo_01304

Solyc11g028010

Novo_01303

Novo_00935

RPKM1

0.107986

0.04235

0.015863

19.92548

0.006346

0.017206

RPKM2

158.3342

63.33414

22.18282

0.021828

31.31789

41.81995

log2FoldChange

10.7241

10.77632

10.67336

-9.55392

12.53154

11.47674

padj

2.46E-221

3.06E-200

4.24E-126

2.38E-125

1.50E-121

3.76E-119

差异基因分析表主要包括的内容：

(1) 基因信息：

包含基因名称（如无参考基因组，不提供）、基因标号、两组实验条件下均一化的测序片段的

度数、RPKM 或者 FPKM 数值 （表达水平）、log

2

Fold Change (相对差异倍数，计算方法用两个

实验条件下的表达水平数值相除) 和 特定分布下的统计量（如 z-score 等，取决于假设分布）。

(2) 统计学检验：

p-value：不同的分析方法所采用的统计学模型不一样，如DESeq（在有重复的条件下，检测差

异基因的效果最好）在检验差异基因时采用负二项分布，p-value 是指在原假设认为在两组实

验条件下的所有基因都不存在差异的情况下出现极端情况的概率。

Corrected p-value : 基因差异表达分析是对于成千上万的基因进行独立的假设统计学检验，这种

多重检验会存在总体假阳性偏高的问题，因此为了降低假阳性，需要对p-value进行校正。q-value

(Benjamini et al. 1995), q-value(Storey et al. 2003)是采用不同的方法对p-value进行校正后的结

果。q-value越低，基因表达差异越显著。

Significant (q-value(Storey et al. 2003) < 0.05) : 当q-value小于默认阈值0.05时，为TRUE，表示

基因表达差异显著，反正为FALSE，无统计学差异。

41

北京诺禾致源生物信息科技有限公司

9.4 差异基因维恩图

差异基因的维恩图代表不同实验条件之间比对得到的差异基因共有的和特异的部分。

如上图所示，不同部分分别代表不同差异基因集共有以及特有的部分。

9.5 差异基因聚类分析

不同的实验条件之间的差异倍数作为输入，用K均值（K-means）或者层次聚类（hierarchical

clustering）的方法对不同表达调控模式进行分类，如下图所示：

不同的颜色的区域代表不同的聚类分组信息，同组内的基因表达调控模式类似，可能参与到相

同的生物学过程中。

42

北京诺禾致源生物信息科技有限公司

10. 差异基因Gene Ontology富集分析

得到差异表达基因之后，我们对差异表达基因做Gene Ontology 富集分析。

Gene Ontology（简称GO）是基因功能国际标准分类体系。根据实验目的筛选差异基因后，研

究差异基因在Gene Ontology中的分布状况将阐明实验中样本差异在基因功能上的体现。GO功能显

著性富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO功能条目，从而给出差异

表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向Gene Ontology数据库

（/）的各个term映射，计算每个term的基因数目，然后应用超几何

检验或者其他统计学方法，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，

下图为超几何分布分计算公式：

公式2：

其中，N为所有基因中具有GO注释的基因数目；n为N中差异表达基因的数目；M为所有基

因中注释为某特定GO term的基因数目；m为注释为某特定GO term的差异表达基因数目。计算得

到的p-value通过Bonferroni校正之后，以corrected p-value≤0.05为阈值，满足此条件的GO term定

义为在差异表达基因中显著富集的GO term。通过GO功能显著性富集分析能确定差异表达基因行

使的主要生物学功能。

分析Gene Ontology的富集分析的软件为GOseq (Young et al.), 除了完成上述分析功能以外，

GOseq还消除了转录本长度对于功能富集结果的影响。

10.1 差异表达基因GO富集分析

样品中差异基因的Gene Ontology 分类

GO ID

GO:0005813

GO:0005634

GO:0005351

GO:0042645

GO:0043533

GO:0005737

GO:0032792

GO:0005896

GO:0016853

GO:0005515

GO Term

centrosome

nucleus

sugar:hydrogen symporter activity

mitochondrial nucleoid

inositol 1,3,4,5 tetra kisphosphate

binding

cytoplasm

negative regulation of CREB

transcription factor activity

interleukin-6 receptor complex

isomerase activity

protein binding

Term

type

C

C

M

C

M

C

P

C

M

M

P-Value

6.50E-06

5.47E-05

0.000130

0.000156

Correct DEG

P-Value item

0.07377

0.31026

0.44389

0.44389

17

108

3

5

Bg/ref

item

252

4569

6

31

2

4488

3

3

111

8310

0.000312 0.67799 2

0.000358 0.67799 105

0.000542 0.87928 2

0.000626 0.88749 2

0.000714 0.90057 7

0.000803 0.91167 172

43

北京诺禾致源生物信息科技有限公司

上述表格为差异基因的Gene Ontology富集分析结果表格。

GO ID: Gene Ontology数据库中唯一的标号信息

GO Term: Gene Ontology功能的描述信息

Term type为该GO的类别（C:细胞组分；P：生物路径；M：分子功能）

P-value: 富集分析统计学显著水平，一般情况下，P-value < 0.05 该功能为富集项

Correct P-Value 矫正后的P-Value

DEG item：差异基因中与该Term相关的基因数

DEG item：所有（bg）基因中与该Term相关的基因数

10.2 GO有向无环图分析

有向无环图为差异基因GO富集分析的结果图形化展示方式，分支代表包含关系，从上至下所

定义的功能范围越来越小，一般选取GO富集分析的结果前10 位作为有向无环图的主节点，并通

过包含关系, 将相关联的GO Term一起展示，颜色的深浅代表富集程度。

11. 差异基因KEGG富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达

基因参与的最主要生化代谢途径和信号转导途径。KEGG（Kyoto Encyclopedia of Genes and Genomes）

是有关Pathway的主要公共数据库（Kanehisa,2008）。Pathway显著性富集分析以KEGG Pathway为单

44

北京诺禾致源生物信息科技有限公司

位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。该

分析的计算公式同公式2中GO功能显著性富集分析，在这里N为所有基因中具有Pathway注释的

基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定Pathway的基因数目；m

为注释为某特定Pathway的差异表达基因数目。FDR≤0.05的Pathway定义为在差异表达基因中显著

富集的Pathway。

11.1 差异基因KEGG显著性富集列表

差异基因KEGG显著性富集列表

Term ID Sample

number

VEGF signaling pathway

Leishmaniasis

Endocytosis

Glycine, serine and threonine metabolism

SNARE interactions in vesicular transport

mmu04370

mmu05140

mmu04144

mmu00260

mmu04130

3 / 58

4 / 58

7 / 58

3 / 58

3 / 58

Background

number

76 / 6865

65 / 6865

224 / 6865

35 / 6865

35 / 6865

0.0260412

0.002117

0.0026149

0.0030921

0.0030921

P-Value Corrected

P-Value

0.0069057

0.0027138

0.0027138

0.0027138

0.0027138

Term：KEGG通路的描述信息。

ID：KEGG数据库中通路唯一的编号信息。

Sample Number：在该通路下的差异基因的个数 / 涉及到KEGG通路中的差异基因的个数。

Background Number：所有通路下的差异基因的个数 / 所有通路下的基因的个数。

P-value：富集分析统计学显著水平。

Corrected P-value：矫正后的统计学显著水平，一般情况下，P-value < 0.05 该功能为富集项。

11.2 差异基因KEGG富集散点图

差异基因KEGG富集散点图中，KEGG富集程度的分布到不同的象限中，横坐标Rich factor 是差

异表达的基因中位于该pathway 条目的基因数目与所有有注释基因中位于该pathway 条目的基因

数的比值。Rich factor越大，表示富集的程度越大。纵坐标-log10（Qvalue）中Qvalue是做过多重假

设检验校正之后的pvalue，-log10（Qvalue）越大，表示富集越显著。我们挑选了富集最显著的20

条pathway条目在该图中进行展示。我们给出的参考重要程度排序如下： 1 Phase > 2 > Phase > 3

Phase > 4 Phase。

45

北京诺禾致源生物信息科技有限公司

12．差异基因蛋白互作网络分析

对差异分析中所产生的差异基因列表，在STRING数据库中找出这些差异基因间的互作关系，并将得

到的互作数据导入 Cytoscape 软件实现互作网络的可视化。蛋白质互作数据来源于STRING数据库

（/）。

互作网络如下图所示，其中节点（node）的大小与此节点的度（degree）成正比，即与此节点相连

的边越多，它的度越大，节点也就越大。节点的颜色与此节点的聚集系数（clustering coefficient）相

关，颜色梯度由绿到红对应聚集系数的值由低到高；聚集系数表示此节点的邻接点之间的连通性好

坏，聚集系数值越高表示此节点的邻接点之间的连通性越好。边（edge）的宽度表示此边连接的两

个节点间的互相作用的关系强弱，互相作用的关系越强，边越宽。

46

北京诺禾致源生物信息科技有限公司

图11 差异基因互作网络图

Anders, S., and Huber, W. Differential expression analysis for sequence count data. Genome Biol 11, R106.

Chepelev, I., Wei, G., Tang, Q., and Zhao, K. (2009). Detection of single nucleotide variations in expressed exons

of the human genome using RNA-Seq. Nucleic acids research 37, e106-e106.

Cock, P.J.A., Fields, C.J., Goto, N., Heuer, M.L., and Rice, P.M. The Sanger FASTQ file format for sequences with

quality scores, and the Solexa/Illumina FASTQ variants. Nucleic acids research 38, 1767-1771.

Dohm, J.C., Lottaz, C., Borodina, T., and Himmelbauer, H. (2008). Substantial biases in ultra-short read data sets

from high-throughput DNA sequencing. Nucleic acids research 36, e105-e105.

Erlich, Y., and Mitra, P.P. (2008). Alta-Cyclic: a self-optimizing base caller for next-generation sequencing. Nature

methods 5, 679-682.

47

北京诺禾致源生物信息科技有限公司

Foissac, S., and Sammeth, M. (2007). ASTALAVISTA: dynamic and flexible analysis of alternative splicing events

in custom gene datasets. Nucleic acids research 35, W297-W299.

Hansen, K.D., Brenner, S.E., and Dudoit, S. Biases in Illumina transcriptome sequencing caused by random

hexamer priming. Nucleic acids research 38, e131-e131.

Hansen, K.D., Wu, Z., Irizarry, R.A., and Leek, J.T. Sequencing technology does not eliminate biological variability.

Nature biotechnology 29, 572-573.

Jensen L J, Kuhn M, Stark M, Chaffron S, Creevey C, Muller J, Doerks T, Julien P, Roth A, Simonovic M, Bork P and

von Mering C. 2009. STRING 8-a global view on proteins and their functional interactions in 630 organisms.

Nucleic Acids Research 37: D412-D416.

Jiang, L., Schlesinger, F., Davis, C.A., Zhang, Y., Li, R., Salit, M., Gingeras, T.R., and Oliver, B. Synthetic spike-in

standards for RNA-seq experiments. Genome research 21, 1543-1551.

Kent, W.J. (2002). BLAT

鈥攖

he BLAST-like alignment tool. Genome research 12, 656-664.

Mamanova, L., Andrews, R.M., James, K.D., Sheridan, E.M., Ellis, P.D., Langford, C.F., Ost, T.W.B., Collins, J.E.,

and Turner, D.J. FRT-seq: amplification-free, strand-specific transcriptome sequencing. Nature methods 7,

130-132.

Marquez, Y., Brown, J.W.S., Simpson, C.G., Barta, A., and Kalyna, M. Transcriptome survey reveals increased

complexity of the alternative splicing landscape in Arabidopsis. Genome research.

Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying

mammalian transcriptomes by RNA-Seq. Nature methods 5, 621-628.

Sammeth, M. (2009). Complete alternative splicing events are bubbles in splicing graphs. Journal of

Computational Biology 16, 1117-1140.

Sammeth, M., Foissac, S., and Guig

贸

, R. (2008). A general definition and nomenclature for alternative splicing

events. PLoS computational biology 4, e1000147.

Trapnell, C., Pachter, L., and Salzberg, S.L. (2009). TopHat: discovering splice junctions with RNA-Seq.

Bioinformatics 25, 1105-1111.

Trapnell, C., Williams, B.A., Pertea, G., Mortazavi, A., Kwan, G., Van Baren, M.J., Salzberg, S.L., Wold, B.J., and

Pachter, L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform

switching during cell differentiation. Nature biotechnology 28, 511-515.

Wang, E.T., Sandberg, R., Luo, S., Khrebtukova, I., Zhang, L., Mayr, C., Kingsmore, S.F., Schroth, G.P., and Burge,

C.B. (2008). Alternative isoform regulation in human tissue transcriptomes. Nature 456, 470-476.

Young, M.D., Wakefield, M.J., Smyth, G.K., and Oshlack, A. goseq: Gene Ontology testing for RNA-seq datasets.

48