诺禾致源有参转录组分析流程

一、建库测序流程

从RNA 样品到最终数据获得，样品检测、建库、测序每一个环节都会对数据质量和数量产生影响，而数据质量又会直接影响后续信息分析的结果。为了从源头上

保证测序数据的准确性、可靠性，诺禾致源对样品检测、建库、测序每一个生产步骤都严格把控，从根本上确保了高质量数据的产出。流程图如下：

1 Total RNA样品检测

诺禾致源对RNA 样品的检测主要包括4种方法：(1) 琼脂糖凝胶电泳分析RNA 降解程度以及是否有污染(2) Nanodrop检测RNA 的纯度（OD260/280比值）(3) Qubit对RNA 浓度进行精确定量(4) Agilent 2100精确检测RNA 的完整性

2 文库构建

样品检测合格后，用带有Oligo （dT ）的磁珠富集真核生物mRNA （若为原核生物，则通过试剂盒去除rRNA 来富集mRNA ）。随后加入fragmentation buffer将mRNA 打断成短片段，以mRNA 为模板，用六碱基随机引物（random hexamers）合成一链cDNA ，然后加入缓冲液、dNTPs 和DNA polymerase I合成二链cDNA ，随后利用AMPure XP beads纯化双链cDNA 。纯化的双链cDNA 再进行末端修复、加A 尾并连接测序接头，然后用AMPure XP beads进行片段大小选择，最后进行PCR 富集得到最终的cDNA

文库。构建原理图如下：

3 库检

文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent 2100对文库的insert size进行检测，insert size符合预期后，使用Q-PCR 方法对文库的有效浓度进行准确定量（文库有效浓度＞2nM ），以保证文库质量。

4 上机测序

库检合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling 后进行HiSeq/MiSeq测序。

二、生物信息分析流程

获得原始测序序列(Sequenced Reads)

后，在有相关物种参考序列或参考基因组的情况下，通过如下流程进行生物信息分析：

其中，DEU 分析仅针对有生物学重复样品，若样品无生物学重复，则不进行此项分析。对于蛋白互作网络分析，若其存在于合同信息分析内容中，而且选择了相应的分析物种或者近缘物种，则进行此项分析；若不存在，则不进行。

三、结果展示及说明

1　原始序列数据

高通量测序(如Illumina HiSeqTM 2500/MiseqTM ) 得到的原始图像数据文件经CASAVA 碱基识别(Base Calling)分析转化为原始测序序列（Sequenced Reads），我们称之为 Raw Data或Raw Reads，结果以 FASTQ (简称为fq) 文件格式存储，其中包含测序序列（reads ）的序列信息以及其对应的测序质量信息。

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

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT +

#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

其中第一行以“@”开头，随后为Illumina 测序标识符(Sequence Identifiers)和描述文字(选择性部分) ；第二行是碱基序列；

第三行以“+”开头，随后为Illumina 测序标识符(选择性部分) ；

第四行是对应碱基的测序质量，该行中每个字符对应的 ASCII 值减去 33，即为对应第二行碱基的测序质量值。Illumina 测序标识符详细信息如下：HWI-ST127671

C1162ACXX [1**********]581N 0CGATGT

Instrument – unique identifier of the sequencerrun number – Run number on instrumentFlowCell ID – ID of flowcellLaneNumber – positive integerTileNumber – positive integer

X – x coordinate of the spot. Integer which can be negativeY – y coordinate of the spot. Integer which can be negativeReadNumber - 1 for single reads; 1 or 2 for paired ends

whether it is filtered - NB：Y if the read is filtered out, not in the delivered fastq file, N otherwisecontrol number - 0 when none of the control bits are on, otherwise it is an even numberIllumina index sequences

　测序数据质量评估

2.1　测序错误率分布检查

每个碱基测序错误率是通过测序Phred 数值(Phred score, Qphred ) 通过公式（公式1 ：Qphred = -10log10(e)）转化得到，而Phred 数值是在碱基识别(Base Calling)过程中通过一种概率模型计算得到，这种模型可以准确地预测碱基判别的错误率。Phred 分值，不正确的碱基识别率，碱基正确识别率以及Q-score 的对应关系如下表所显示：

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

Phred 分值

10203040

不正确的碱基识别

1/101/1001/10001/10000

碱基正确识别率

90%99%99.9%99.99%

Q-sorce

Q10Q20Q30Q40

测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于RNA-seq 技术，测序错误率分布具有两个特点：

(1)测序错误率会随着测序序列(Sequenced Reads)长度的增加而升高，这是由于测序过程中化学试剂的消耗导致的，并且为illumina 高通量测序平台都具有的特征。(2)前6个碱基的位置也会发生较高的测序错误率，而这个长度也正好等于在RNA-seq 建库过程中反转录所需要的随机引物的长度。所以前6个碱基测序错误率较高的原因为随机引物和RNA 模版的不完全结合(Jiang et al.)。

图1　测序错误率分布图

横坐标为reads 的碱基位置，纵坐标为单碱基错误率

2.2

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

GC 含量分布检查用于检测有无AT 、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。

在illumina 测序平台的转录组测序中，反转录成cDNA 时所用的6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。而这种偏好性与测序的物种和实验室环境无关，但会影响转录组测序的均一性程度(Hansen et al.)。除此之外，理论上普通文库的G 和C 碱基及A 和T 碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线，而对于链特异性建库会出现GC 分离的现象。对于DGE 测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read 前6-7个碱基有较大的波动，这种波动属于正常情况。

图2　GC 含量分布图

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

2.4　测序数据质量情况汇总

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

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

Sample name

sampleA1_1sampleA1_2sampleA2_1sampleA2_2sampleB1_1sampleB1_2sampleB2_1sampleB2_2sampleC1_1sampleC1_2sampleC2_1sampleC2_2

数据质量情况详细内容如下：

Raw reads

[***********][***********][***********][***********][***********]521158

Clean reads

[***********][***********][***********][***********][***********]035877

Clean bases

4.19G 4.19G 4.14G 4.14G 3.92G 3.92G 3.65G 3.65G 3.55G 3.55G 3.5G 3.5G

Error rate(%)

0.030.040.030.040.030.040.030.040.030.040.030.04

Q20(%)

95.6194.1395.8893.9895.9094.2395.9294.0595.9594.3995.9494.28

Q30(%)

91.3289.0791.8788.8491.8889.2491.9288.9491.9889.5092.0089.38

GC content(%)

48.4348.4748.4248.4449.3049.3048.8748.8848.5748.5747.7247.72

(1) Sample name：样品名，1为左端reads ，2为右端reads 。样品的 clean reads 总数为左端+右端。(2) Raw reads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。

(3) Clean reads：计算方法同 Raw Reads，只是统计的文件为过滤后的测序数据。后续的生物信息分析都是基于Clean reads。(4) Clean bases：Clean reads的个数乘以长度，并转化为以G 为单位。(5) Error rate：通过公式1计算得到。

(6) Q20、Q30：分别计算 Phred 数值大于20、30的碱基占总体碱基的百分比。(7) GC content：计算碱基G 和C 的数量总和占总的碱基数量的百分比。

2.5 质量评估Q&A

问：测序错误率会随着测序序列长度的增加而升高，错误率在多少是可以接受的范围？

答：诺禾的测序会进行严格的数据质量把控。一般情况下，单个碱基位置的测序错误率应该低于1%，最高在6%左右可以接受。问：诺禾质控的标准是什么？是否严格？

答：为保证后续分析的质量，诺禾会严格把控cleandata 的筛选标准，具体标准如下： (1) 去除带接头(adapter)的reads ；

(2) 去除N(N表示无法确定碱基信息) 的比例大于10%的reads ；

(3) 去除低质量reads(质量值sQ

adapter ：接头，用于上机测序。建库时引入的接头序列与测序芯片（flow cell）上固定的接头相互识别。 index ：测序的标签，用于测定混合样本，通过每个样本添加的不同标签进行数据区分，鉴别测序样品。 Q20,Q30：Phred 数值大于20、30的碱基占总体碱基的百分比，其中Phred=-10log10(e). raw data/raw reads：测序下机的原始数据。

clean data/clean reads：对原始数据进行过滤后，剔除了低质量数据的剩余数据。后续分析均基于clean data。

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

3　参考序列比对分析

测序序列定位算法：根据不同的基因组的特征，我们选取相对合适的软件(动植物用TopHat2、细菌或者基因密度较高的物种用Bowtie2, 二者mismatch 均设为2，其余选用默认参数) ，合适的参数设置(如最大的内含子长度，会根据已知的该物种的基因模型来进行统计分析) ，将过滤后的测序序列进行基因组定位分析。下图为TopHat2

的算法示意图：

TopHat2的算法主要分为三个部分：(1) 将测序序列和转录组进行比对(可选) (2) 将测序序列整段比对到基因组外显子上(3) 将测序序列分段比对到基因组的两个外显子上

如果参考基因组选择合适，而且相关实验不存在污染，实验所产生的测序序列的定位的百分比正常情况下会高于70% (Total Mapped Reads or Fragments)，其中具有多个定位的测序序列(Multiple Mapped Reads or Fragments)占总体的百分比通常不会超过10%。

3.2

　Reads 在参考基因组不同区域的分布情况

将比对到基因组上的reads 分布情况进行统计，定位区域分为Exon(外显子) 、Intron(内含子) 和Intergenic(基因间区) 。

在基因组注释较为完全的物种中，比对到Exon （外显子）的reads 含量最高，比对到Intron （内含子）区域的reads 来源于pre-mRNA 的残留及可变剪切过程中发生的内含子滞留事件导致的，而比对到Intergenic （基因间区）的reads 是因为基因组注释不完全。

图3.2　Reads 在参考基因组不同区域的分布情况

3.3

　Reads 在染色体上的密度分布情况

对Total mapped reads的比对到基因组上的各个染色体（分正负链）的密度进行统计，如下图所示，具体作图的方法为用滑动窗口(window size)为1K ，计算窗口内部比对到碱基位置上的reads 的中位数，并转化成 log2 。正常情况下，整个染色体长度越长，该染色体内部定位的reads 总数会越多(Marquez et al. 2012)。从定位到染色体上的reads 数与染色体长度的关系图中，可以更加直观看出染色体长度和reads 总数的关系。

图3.3　Reads 在染色体上的密度分布图

每个样品两张图；图一：横坐标为染色体的长度信息(以百万碱基为单位) ，纵坐标为log 2(reads的密度的中位数) ，绿色为正链，红色为负链图二: 横坐标为染色体的长度信息(单位为

Mb) ，纵坐标为mapped 到染色体上的reads 数(单位为M) ，图中灰色区域表示95%的置信区间

3.4　Reads 比对结果可视化

我们提供RNA-seq Reads在基因组上比对结果的bam 格式文件，部分物种还提供相应的参考基因组和注释文件，并推荐使用IGV (Integrative Genomics Viewer) 浏览器对bam 文件进行可视化浏览。IGV 浏览器具有以下特点：(1)能在不同尺度下显示单个或多个读段在基因组上的位置，包括读段在各个染色体上的分布情况和在注释的外显子、内含子、剪接接合区、基因间区的分布情况等；(2)能在不同尺度下显示不同区域的读段丰度，以反映不同区域的转录水平；(3)能显示基因及其剪接异构体的注释信息；(4)能显示其他注释信息；(5)既可以从远程服务器端下载各种注释信息，又可以从本地加载注释信息。IGV 浏览器使用方法可参考我们提供的使用说明文档(IGVQuickStart.pdf )

。

图3.4　IGV 浏览器界面

3.5 参考序列比对分析 Q&A问：有参分析都需要什么文件？

答：相应的参考基因组及基因结构注释文件（gtf/gff/gff3/bed等格式，推荐gtf ，gff ）、基因的GO 注释文件的直接下载链接以及基因功能描述文件。问：造成mapping rate较低的原因可能有哪些？

答：TopHat 比对时，默认为2个mismatch ，即：reads 和reference 在2mismatch 之内，就算mapping 到了。当mappingrate 较低时主要可能有2个原因：（1）由于reference 组装不好，或者所测物种与reference 的亲缘关系较远；（2）由于样品的特殊前处理或者相对于参考基因组此样品本身的变异太大，导致mapping rate相对较低。

问：mapping 时用的是read 全长，还是头尾有处理？

答：实验方面，我们使用标准的RNA-seq 试剂盒，其index 处于Adapter 中间，在测序中由Index read完成，由此测序得到的Read 1和Read 2的各个碱基全都是样本的序列，因此mapping 时，头尾不处理。信息分析方面，我们会将过滤得到的clean reads的全长进行mapping 。问：Read-1,Read-2的具体含义是什么？

答：双端测序会在cDNA 片段的两端先后读取一定长度的碱基（如pe125，就是分别测125bp ）得到一对reads ，其中一条称为Read-1，另一条称为Read-2。大量文献和实际项目都表明，转录组分析中使用双端reads 对于序列拼接和定量都大有裨益。

4　可变剪切分析

用ASprofile 软件对Cufflinks (Trapnell et al.)预测出的基因模型对每个样品的可变剪切事件分别进行分类和表达量统计。可变剪切分析流程及ASprofile 中的可变剪切

事件分类如下图所示：

12类可变剪切事件定义如下:

(1) TSS: Alternative 5' first exon (transcription start site) 第一个外显子可变剪切 (2) TTS: Alternative 3' last exon (transcription terminal site)　最后一个外显子可变剪切 (3) SKIP: Skipped exon

(SKIP_ON,SKIP_OFF pair)　单外显子跳跃 (4) XSKIP: Approximate SKIP (XSKIP_ON,XSKIP_OFF pair)　单外显子跳跃（模糊边界） (5) MSKIP: Multi-exon SKIP (MSKIP_ON,MSKIP_OFF pair)　多外显子跳跃 (6) XMSKIP: Approximate MSKIP (XMSKIP_ON,XMSKIP_OFF pair)　多外显子跳跃（模糊边界） (7) IR: Intron retention (IR_ON, IR_OFF pair)　单内含子滞留 (8) XIR: Approximate IR (XIR_ON,XIR_OFF pair)　单内含子滞留（模糊边界） (9) MIR: Multi-IR (MIR_ON, MIR_OFF pair)　多内含子滞留 (10) XMIR: Approximate MIR (XMIR_ON, XMIR_OFF pair)　多内含子滞留（模糊边界） (11)AE: Alternative exon ends (5', 3', or both)　可变 5'或3' 端剪切 (12) XAE: Approximate AE　可变 5'或3' 端剪切（模糊边界）

图4.1　AS 分类和数量统计

纵轴为可变剪切事件的分类缩写，横轴为该种事件下可变剪切的数量，不同样品用不同子图和颜色区分

表4.2　AS 结构和表达量统计

event_idevent_typegene_idchromevent_startevent_end

[***********]0031000004

TSS TTS SKIP_ONSKIP_OFF

CUFF.1000CUFF.1000CUFF.1000CUFF.1000

1111

[***********][***********]

event_pattern

[***********]130657314,130657463-130657479,[***********],130661448

strand

----

fpkm ref_id

0.0000000000ENSMUST[1**********]0.0000000000ENSMUST[1**********]0.0000000000ENSMUST[1**********]0.0000000000ENSMUST[1**********]

(1) event_id: AS事件编号

(2) event_type: AS事件类型 (TSS, TTS, SKIP_{ON,OFF}, XSKIP_{ON,OFF}, MSKIP_{ON,OFF}, XMSKIP_{ON,OFF}, IR_{ON ,OFF}, XIR_{ON,OFF}, AE, XAE)(3) gene_id: cufflink组装结果中的基因编号(4) chrom: 染色体编号

(5) event_start: AS事件起始位置(6) event_end: AS事件结束位置

(7) event_pattern: AS事件特征 (for TSS, TTS - inside boundary of alternative marginal exon; for *SKIP_ON,the coordinates of the skipped exon(s); for *SKIP_OFF, the coordinates of the enclosing introns; for*IR_ON, the end coordinates of the long, intron-containing exon; for *IR_OFF, the listing of coordinates of all the exons along the path containing the retained intron; for *AE, the coordinates of the exonvariant)

(8) strand: 基因正负链信息

(9) fpkm: 此AS 类型该基因表达量

(10) ref_id: 此基因在参考注释文件中的编号

4.3 可变剪切分析 Q&A

问：什么是可变剪切(alternative splicing)？

答：大多数真核基因转录产生的mRNA 前体是按一种方式剪接产生出一种mRNA ，因而只产生一种蛋白质。但有些基因产生的mRNA 前体可按不同的方式剪接，产生出两种或更多种mRNA ，即可变剪接(alternative splicing)。

5　新转录本预测

5.1　新转录本预测

将所有测序reads 数据的基因组定位结果放到一起，用 Cufflinks 进行组装，然后用Cuffcompare 和已知的基因模型进行比较，可以:(1)发现新的未知基因（相对于原有基因注释文件）；(2)发现已知基因新的外显子区域；(3)对已知基因的起始和终止位置进行优化。新基因和新外显子区域预测结果为GTF 格式的注释文件。GTF 格式的详细说明可参考(http://mblab.wustl.edu/GTF22.html)

表5.1　新转录本结构注释结果

seqname

1111

source

novelGene novelGene novelGene novelGene

feature

exon exon exon exon

start

[***********][1**********]178

end

[***********][1**********]237

score strand frame

. . . .

++++

. . . .

attributes

gene_id "Novel00001"; transcript_id "Novel00001.1"; exon_number "2";gene_id "Novel00001"; transcript_id "Novel00001.1"; exon_number "3";gene_id "Novel00001"; transcript_id "Novel00001.1"; exon_number "4";gene_id "Novel00001"; transcript_id "Novel00001.2"; exon_number "1";

(1) seqname：染色体编号

(2) source：来源标签，这里的novelGene 指新基因(3) feature：区域类型，目前我们预测外显子区域(4) start：起始坐标(5) end：终止坐标(6) score：不必关注(7) strand：正负链信息(8) frame：不必关注

(9) attributes：属性，包括基因编号、转录本编号等信息

5.2　基因结构优化

表5.2　已知基因结构优化

Gene_id

Chromosome

ENSMUSG[1**********]3ENSMUSG[1**********]X ENSMUSG[1**********]11ENSMUSG[1**********]

(1) Gene_id：原注释文件中基因命名编号(2) Chromosome：染色体编号(3) Strand：正负链信息

(4) Original_span：原注释文件中基因起始位置~终止位置

(5) Assembled_span：转录组拼接结果中基因起始位置~终止位置

Strand

Original_span

Assembled_span

-108107280~[***********]~108154921+161117193~[***********]~161258395+121237253~[***********]~121258656+

17281185~17289115

17280908~17289115

5.3 新转录本预测 Q&A问：新转录本预测的意义？

答：我们使用cufflinks 拼接得到的基因注释文件，与原有基因注释文件进行比较，找出原有注释中没有包含的基因并对基因的位置进行优化，补充并修改了原有的注释文件。

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

6.2 SNP和InDel Q&A

问：SNP 分析的参考序列是什么，即REF 是指什么？

答：参考序列REF 是选取的参考基因组序列，SNP 是通过将reads 比对到参考基因组上从而进行SNP calling。问：SNP 中一列中两个数字分别代表支持REF 和ALT 两种碱基的reads 数目，为什么有些会是0呢？

答：以0,12为例，表示有0个reads 支持REF 的碱基，即没有支持该碱基的reads ，12个reads 支持ALT 的碱基。问：SNP 具体指什么？其与InDel 有何区别？

答：一般而言，SNP 是指变异频率大于1%的单核苷酸变异。InDel(insertion-deletion)则是插入或者缺失，insert 或者deletion 。问：SNP 和InDel 相关名词的解释

SNP ：SNP （Single Nucleotide Polymorphisms）单核苷酸多态性。

SNP calling：查找NGS 数据与参考序列区别的过程，称为SNP calling。其中包含统计矩阵的计算，以筛选出最可能的SNP 。 InDel ：插入缺失。是指相对于参考基因组，样本中发生的小片段的插入缺失，该插入缺失可能含一个或多个碱基。

7.2 基因表达水平分析Q&A问：基因表达水平如何计算？

答：在RNA-seq 技术中，FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) 是每百万fragments 中来自某一基因每千碱基长度的fragments 数目,FPKM 时考虑了测序深度和基因长度对reads 计数的影响，是目前最为常用的基因表达水平估算方法。问：认为基因表达的阈值是多少？为什么设置为这个阈值？

答：有参转录组当中，认为FPKM>1是基因表达的。这个阈值是主流杂志推荐的，也能够很好的反应基因的表达水平。基因表达水平分析相关名词的解释:

FPKM ：expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced，是每百万fragments 中来自某一基因每千碱基长度的fragments 数目。

RPKM ：expected number of reads Per Kilobase of transcript sequence per Millions base pairs sequenced，是每百万reads 中来自某一基因每千碱基长度的reads 数目。

　RNA-seq 整体质量评估

8.1 RNA-Seq相关性检查

生物学重复是任何生物学实验所必须的，高通量测序技术也不例外(Hansen et al.)。生物学重复主要有两个用途：一个是证明所涉及的生物学实验操作是可以重复的且变异不大，另一个是为了确保后续的差异基因分析得到更可靠的结果。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1，表明样品之间表达模式的相似度越高。Encode 计划建议皮尔逊相关系数的平方(R2) 大于0.92(理想的取样和实验条件下) 。具体的项目操作中，我们要求R 2至少要大于0.8，否则需要对样品做出合适的解释，或者重新进行实验。

图8.1　RNA-Seq 相关性检查

热图：样品间相关系数热图；散点图(若样品多于4组，则仅展示生物学重复之间的散点图) ：样品间的相关系数散点图，R 2:pearson相关系数的平方。

8.2 RNA-seq整体质量评估Q&A

问：样品间的相关性有何意义？如何计算?

答：样品间的相关性反应了样品间的相似程度，即不同处理或组织的样品在表达水平方面的相似度。相关系数越接近1，样品间的相似度越高，样品间的差异基因也越少。生物学重复间的样品的相关系数应大于生物学重复外的样品的相关系数。相关系数的计算方法有三种：A. Pearson correlation; B. Spearman rankcorrelation; C. Kendall’s τ。诺禾使用R 语言进行Pearson 相关系数的计算。

差异表达分析

9.1　基因表达水平对比

通过所有基因的FPKM 分布图以及盒形图对不同实验条件下的基因表达水平进行比较。对于同一实验条件下的重复样品，最终的FPKM 为所有重复数据的平均值。

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

图一：FPKM 盒形图，横坐标为样品名称，纵坐标为log 10(FPKM+1)，每个区域的盒形图对五个统计量(至上而下分别为最大值，上四分位数，中值，下四分位数和最小值) 图二：FPKM 分

布图，横坐标为log 10(FPKM+1), 纵坐标为基因的密度

9.2　差异表达基因列表

基因差异表达的输入数据为基因表达水平分析中得到的readcount 数据。对于有生物学重复的样品，我们采用DESeq （Anders et al, 2010）进行分析：该分析方法基于的模型是负二项分布，第 i 个基因在第 j 个样本中的 read count 值为K ij ，则有

K ij ～ NB(µij , σij 2)

对于无生物学重复的样品，先采用TMM 对read count数据进行标准化处理，之后用DEGseq 进行差异分析。差异表达基因列表如下：

表9.2　差异基因列表

Gene Id

ENSMUSG[1**********]ENSMUSG[1**********]ENSMUSG[1**********]ENSMUSG[1**********]

差异基因列表主要包括的内容：(1) Gene id: 基因编号

(2) Sample1：校正后样品1的readcount 值(3) Sample2：校正后样品2的readcount 值(4) log2FoldChange: log2(Sample1/Sample2)

sampleA

27459.[**************].[1**********]97.[**************]22.5175799308

sampleB

3588.[**************]00.[1**********]33.[**************].[1**********]

log2FoldChange

2.936-2.789-2.86772.3387

pval

1.5913e-1709.2315e-1086.3681e-1089.3935e-106

p-adjusted

3.2181e-1666.223e-1046.223e-1044.7491e-102

(5) pvalue(pval): 统计学差异显著性检验指标

(6) qvalue(p-adjusted): 校正后的pvalue 。qvalue 越小，表示基因表达差异越显著

9.3

　差异表达基因筛选

用火山图可以推断差异基因的整体分布情况，对于无生物学重复的实验，为消除生物学变异，从差异倍数和显著水平两个方面进行评估，对差异基因进行筛选，阈值设定一般为: |log2(FoldChange)| > 1 且 qvalue

图9.3　差异基因火山图

有显著性差异表达的基因用红色点(上调) 和绿色点（下调）表示，无显著性差异表达的基因用蓝色点表示；横坐标代表基因在不同样本中表达倍数变化；纵坐标代表基因表达量变化差异

的统计学显著性

9.4

　差异基因聚类分析

聚类分析用于判断差异基因在不同实验条件下的表达模式；通过将表达模式相同或相近的基因聚集成类，从而识别未知基因的功能或已知基因的未知功能；因为这些同类的基因可能具有相似的功能，或是共同参与同一代谢过程或细胞通路。以不同实验条件下的差异基因的FPKM 值为表达水平，做层次聚类(hierarchicalclustering) 分析，不同颜色的区域代表不同的聚类分组信息，同组内的基因表达模式相近，可能具有相似的功能或参与相同的生物学过程。

除了差异基因表达量FPKM 层次聚类分析，我们还分别用H-cluster 、K-means 和SOM 等三种方法对差异基因的相对表达水平值log 2(ratios)进行聚类。不同的聚类算法分别将差异基因分为若干cluster ，同一cluster 中的基因在不同的处理条件下具有相似的表达水平变化趋势。

图9.4　差异基因聚类图

图一：整体FPKM 层次聚类图，以log 10(FPKM+1)值进行聚类，红色表示高表达基因，蓝色表示低表达基因。颜色从红到蓝，表示log 10(FPKM+1)从大到小图二：log 2(ratios)折线图每个子图中的灰色线条表示一个cluster 中的基因在不同实验条件下相对表达量，蓝色线条表示这个cluster 中的所有基因在不同实验条件下相对表达量的平均值，x 轴表示实验条件，y 轴表示相对

表达量

9.5

　差异基因维恩图

差异基因维恩图展示了各比较组间差异基因的个数，以及比较组间的重叠关系。(仅提供两组、三组和四组比较的venn 图)

图9.5　差异基因维恩图

每个大圆圈中的数字之和代表该比较组合的差异基因总个数，圆圈交叠的部分表示组合之间共有的差异基因。

9.6 差异表达分析Q&A

问：如何判断一个基因是否是差异基因？如果是差异基因，如何判断该基因的表达量是上调还是下调？

答：无生物学重复时，先采用TMM 对read count数据进行标准化处理，再用DEGseq 进行差异分析, 筛选阈值为qvalue

答：在做差异分析时，诺禾是采用readcount 数据，通过DESeq 或者TMM 标准化后，进行差异分析。FPKM 实际上也是对readcount 进行标准化处理的一种方法，各种标准化方法优劣势比较见下图(Dillies, M. A. et al, 2013)，可以看出，DESeq 和TMM 的标准化效果最好，FPKM 的标准化效果较差，所以，不推荐使用FPKM 进行

差异分析。

问：如何判断差异基因在两个样品间的差异大小？

答：差异的显著情况可通过矫正后的pvalue 来看，矫正后的pvalue 越小，差异越显著。也可通过|log2Foldchange|来判断差异的大小情况，|log2Foldchange|越大，差异倍数越大。

问：某基因在两个样本中表达量差别很大，却不存在与显著差异的基因列表中，这是为何？

答：差异基因的筛选是基于统计学意义的，不能直观的通过两个数值的大小判断差异基因的是否：首先：受测序深度的影响，有些样品的测序深度较深，可能导致该样品的首先：readcount 数值较高，做差异分析的第一步就是要消除测序深度的影响，对原始数据进行标准化处理（我们在有重复项目中，使用DESeq 自带的标准化方法；无重复项目中，使用TMM 标准化方法）其次：在差异分析过程中，需要对其次：readcount 的分布进行估计，经验表明，readcount 服从负二项分布。在有重复的项目中，重复的好坏也会对差异基因与否产生影响。如果重复较差，组内差异情况会屏蔽掉部分组间的差异。在估计完参数后，需要用特定检验方法来判断差异基因与否再次：在计算完再次：pvalue 以后，需要对pvalue 进行多重假设检验校正，来减少假阳性。这个过程会使得padj 会大于原来的pvalue ，使得部分通过pvalue 阀值的基因，无法通过padj 的阀值。

问：差异基因筛选条件最大能设的阈值是多少？很多客户希望通过调整差异基因筛选阈值来找相关基因是否有必要？

答：一般来说，等级较高的文章阈值的设置会比较严格；而在某些文章中，差异基因筛选阈值会适当放宽：如在一些无生物学重复的文章中，只将qvalue 作为差异基因筛选标准，不考虑log2foldchange ；有的文章则将pvalue 作为差异基因的筛选标准。问：某基因readcount 值为0，但是也有foldchange 以及pvalue 、qvalue 值？

答：在DESeq 中，如果某基因的在一个样品中的校正后的readcount 为0，而在另一个样品中不为0，foldchange 会为INF 或者-INF ；如果两个数值均为

0，log2foldchange 以及pvalue 、qvalue 值均为NA ；在DEGseq 中，如果某基因的在一个样品中的校正后的readcount 为0，软件会默认的把0进行轻微的校正，校正成一个接近于0，但不为0的值，故会产生foldchange 以及pvalue 、qvalue 值。

问：差异基因列表中，readcount 一个为0，另一个不为0，能否说明一个表达，一个不表达？答：在有参项目中，一般默认rpkm>1时，基因表达。一般不推荐看readcount 的值看判断表达与否。问：能否提取部分基因来做差异分析？

答：不能。差异分析是基于整体来做的。差异分析软件的作者推荐用全部readcount 进行差异分析，若使用部分基因做分析，会毁坏掉数据整体的特点，如测序深度、reads 分布特征。所以不推荐老师抽取部分来做差异分析。问：按照指定的顺序画聚类图可以吗？

答：h_cluster,k-means和som 的聚类图可按照指定顺序绘制，但聚类热图的顺序不能调整。问：聚类分析是怎么做的？

答：聚类使用的为R 中的聚类软件包pheatmap ，所针对的数据为union_for_cluster（差异基因的并集），以基因的相对表达水平值log2(ratios) 进行聚类。其采用相应的距离算法，算出每个基因之间的距离，然后通过反复迭代，计算基因之间的相对距离，最后根据基因的相对距离远近来分成不同的subcluster ，从而实现聚类。该软件包是免费的，只需通过R 来运行。H-cluster 、K-means 和SOM 都是聚类的方法，均采用的是R 语言相关代码和函数实现的，也有一些免费的软件可以做这些聚类分析，例如gene_cluster等。问：为什么进行聚类分析？

答：聚类分析用于判断差异基因在不同实验条件下的表达模式，将表达模式相同或相近的基因聚集成类，进而识别未知基因的功能或已知基因的未知功能，这些同类基因可能具有相似的功能，共同参与同一代谢过程或存在于同一细胞通路中。问：为什么要用校正后的p 值，能直接用p_value吗？

答：校正后的p 值（padj/qvalue），是对p 值进行了多重假设检验，校正后的p 值比原有的p 值要大，能更好地控制假阳性率。

10　差异基因GO 富集分析

Gene Ontology（简称 GO, http://www.geneontology.org/）是基因功能国际标准分类体系。作为基因本体联合会(Gene Onotology Consortium)所建立的数据库，它旨在建立一个适用于各种物种的，对基因和蛋白质功能进行限定和描述的，并能随着研究不断深入而更新的语言词汇标准。GO 分为分子功能（Molecular Function）、生物过程（biological process）、和细胞组成（cellular component）三个部分。基因或蛋白质可以通过ID 对应或者序列注释的方法找到与之对应的GO 编号，而GO 编号可用于对应到Term ，即功能类别或者细胞定位。根据实验目的筛选差异基因后，富集分析研究差异基因在 Gene Ontology 中的分布状况以期阐明实验中样本差异在基因功能上的体现。普通GO 富集分析的原理为超几何分布：根据挑选出的差异基因计算这些差异基因同GO 分类中某几个特定的分支的超几何分布关系，通过假设验证得到一个特定p-value 。小的p 值表示差异基因在该GO 中出现了富集。

我们在分析中GO 富集分析采用的软件方法为GOseq （Young et al, 2010）, 此方法基于 Wallenius non-central hyper-geometric distribution。相对于普通的超几何分布(Hyper-geometric distribution)，此分布的特点是从某个类别中抽取个体的概率与从某个类别之外抽取一个个体的概率是不同的，而这种概率的不同是通过对基因长度的偏好性进行估计得到的，从而能更为准确地计算出 GOterm 被差异基因富集的概率。10.1　差异基因GO 富集列表

表10.1　差异基因GO 富集列表

GO accession

GO:0005488GO:0036094GO:0000166GO:1901265

Description

binding

small molecule bindingnucleotide bindingnucleoside phosphate binding

Term type

molecular_functionmolecular_functionmolecular_functionmolecular_function

Over represented p-Value

4.847e-071.7728e-052.9533e-052.9533e-05

Corrected p-ValueDEG itemDEG list

0.00197760.0210210.0210210.021021

[1**********]63

[**************]9

结果表格详细内容如下：

(1) GO accession：Gene Ontology数据库中唯一的标号信息(2) Description：Gene Ontology功能的描述信息

(3) Term type：该GO 的类别(cellular_component：细胞组分；biological_process：生物学过程；molecular_function：分子功能) (4) Over represented p-Value：富集分析统计学显著水平

(5) Corrected p-Value：矫正后的P-Value ，一般情况下，Corrected_pValue

10.2

　差异基因GO 富集柱状图

差异基因GO 富集柱状图，直观的反映出在生物过程(biological process)、细胞组分(cellular component)和分子功能(molecular function)富集的GO term上差异基因的个数分布情况。我们挑选了富集最显著的30个GO term在图中展示，如果不足30条，则全部展示。

图10.2　GO 富集柱状图

每组两张图；图一：纵坐标为富集的GO term，横坐标为该term 中差异基因个数。不同颜色用来区分生物过程、细胞组分和分子功能，带“*”为显著富集的GOterm 图二：对图一中的GO ，

按生物过程、细胞组分和分子功能三大类别及差异基因上下调分类画的三个子图

10.3

　差异基因GO 富集DAG 图

有向无环图(Directed Acyclic Graph，DAG) 为差异基因GO 富集分析结果的图形化展示方式。图中，分支代表包含关系，从上至下所定义的功能范围越来越小，一般选取GO 富集分析的结果前10位作为有向无环图的主节点，并通过包含关系，将相关联的GO Term一起展示，颜色的深浅代表富集程度。我们的项目中分别绘制生物过程(biological process)、分子功能(molecular function)和细胞组分(cellular component)的DAG 图。

图10.3　GO 富集有向无环图

每个节点代表一个GO 术语，方框代表的是富集程度为TOP10的GO ，颜色的深浅代表富集程度，颜色越深就表示富集程度越高，每个节点上展示了该TERM 的名称及富集分析的Corrected

p-value

10.4 差异基因GO 富集分析Q&A问：GO 富集分析所使用的软件是什么？

答：GO 富集分析均采用R 包，富集采用的为GOseq, 有向无环图采用的为topGO 。问：GO 富集分析一般分析到二级，还能继续分析吗？

答：GO 富集分析是针对所有注释的GO 进行统计检验，任何等级的都有。问：GO 富集和GO 分类有何区别？

答：GO 分类是将每个基因与其对应的GO 功能联系起来，以获取基因的GO 注释信息，而GO 富集分析则是将GO 功能相似的基因集通过统计学检验算法富集到一起，从而方便研究具有某一类GO 功能的基因。

11　差异基因KEGG 富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway 显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。KEGG （Kyoto Encyclopedia of Genes and Genomes）是系统分析基因功能、基因组信息数据库，它有助于研究者把基因及表达信息作为一个整体网络进行研究。作为Pathway 相关的主要公共数据库(Kanehisa,2008）) ，KEGG 提供的整合代谢途径 (pathway)查询十分出色，包括碳水化合物、核苷、氨基酸等的代谢及有机物的生物降解，不仅提供了所有可能的代谢途径，而且对催化各步反应的酶进行了全面的注解，包含有氨基酸序列、PDB 库的链接等等，是进行生物体内代谢分析、代谢网络研究的强有力工具。Pathway 显著性富集分析以KEGG 数据库中Pathway 为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway 。计算公式如下

在这里N 为所有基因中具有Pathway 注释的基因数目； n为N 中差异表达基因的数目；M 为所有基因中注释为某特定Pathway 的基因数目； m 为注释为某特定Pathway 的差异表达基因数目。FDR≤0.05 时，表示差异基因在该Pathway 中显著富集，我们使用KOBAS （2.0）进行Pathway 富集分析。11.1　差异基因KEGG 富集列表

表11.1　差异基因KEGG 富集列表

#Term

Axon guidanceGlutamatergic synapseRetrograde endocannabinoid

signaling

Circadian entrainment结果表格详细内容如下：

Database

KEGG PATHWAY KEGG PATHWAY KEGG PATHWAY KEGG PATHWAY

mmu04360mmu04724mmu04723mmu04713

Sample number

53423735

Background number

[1**********]

P-Value

2.[1**********]e-060.[**************]0.[**************]0.[1**********]844

Corrected P-Value

0.[**************]0.[1**********]160.[1**********]660.[1**********]66

(1) #Term：KEGG 通路的描述信息。(2) Database：KEGG 数据库。

(3) ID：KEGG 数据库中通路唯一的编号信息。(4) Sample number：该通路下差异基因的个数。

(5) Background number：该通路下注释基因的个数。(6) P-value：富集分析统计学显著水平。

(7) Corrected P-value：矫正后的统计学显著水平，Corrected P-value

11.2

　差异基因KEGG 富集散点图

散点图是KEGG 富集分析结果的图形化展示方式。在此图中，KEGG 富集程度通过Rich factor、Qvalue 和富集到此通路上的基因个数来衡量。其中Rich factor指该pathway 中富集到的差异基因个数与注释基因个数的比值。Rich factor越大，表示富集的程度越大。Qvalue 是做过多重假设检验校正之后的Pvalue ，Qvalue 的取值范围为[0,1]，越接近于零，表示富集越显著。我们挑选了富集最显著的20条pathway 条目在该图中进行展示，若富集的pathway 条目不足20条，则全部展示。

图11.2　差异基因KEGG 富集散点图

纵轴表示pathway 名称，横轴表示Rich factor，点的大小表示此pathway 中差异表达基因个数多少，而点的颜色对应于不同的Qvalue 范围

11.3

富集KEGG 通路图

为便于查看差异基因在通路图中的分布情况，我们将差异基因标注到通路图中，查看方法如下：拿到全部分析结果后，打开results 结果目录中

DEG_KEGGEnrichment/DEG_KEGGPath文件夹下相应比较组合的html 文件，点击不同的通路名称会弹出相应的通路图。其中，包含上调基因的KO 节点标红色，包含下调基因的KO 节点标绿色，包含上下调的标黄色。鼠标悬停于标记的KO 节点，弹出差异基因细节框，标色同上，括号中数字为log2(Fold change)。以上步骤可脱机实现，如连接互联网，点击各个节点，可以连接到KEGG 官方数据库中各个KO 的具体信息页。

图11.3 显著富集的KEGG pathway代谢通路图

11.4 差异基因KEGG 富集分析Q&A

问：为什么编码同一个酶的基因，会有的上调有的下调？

答：这些编号的基因存在着多个条目，也可能包含了一个家族的多个基因，它们间的调控机制可能尚不清楚，反映在图上会有部分上调，部分下调的现象，这是比较常见的现象。

问：KEGG Pathway图片中节点及底色有何含义？

答：KEGG Pathway共有5种类型, 分别为： map: Reference pathway ko: Reference pathway (KO) ec: Reference pathway (EC)

rn: Reference pathway (Reaction) org: Organism-specific pathway map 图中节点及底色的含义如下：

1、"map" pathway: 节点代表某一基因、该基因编码的酶及这个酶参与的反应，将鼠标悬停在某个节点上，可以看到上述信息，底色以无色表示。

2、"ko/ec/rn" pathway： ko 通路中的节点只代表基因；ec 通路中的节点只代表相关的酶；rn 通路中的节点只表示该点参与的某个反应、反应物及反应类型。底色以蓝色表示。

3、"org" pathway：物种特异性通路图。节点含义同map pathway。物种相关的节点底色以绿色表示。

每个数据路的pathway 都有相应的唯一编号，如map00010，据此可在kegg 数据库官网查询。有参考基因组项目中得到的通路图都是Organism-specific pathway 。

12　蛋白互作网络分析

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

我们主要应用STRING 蛋白质互作数据库(http://string-db.org/)中的互作关系进行差异基因蛋白互作网络的分析。，针对于数据库中包含的物种，直接从数据库中提取出目标基因集(比如差异基因list) 的互作关系来构建网络；针对于数据库中不包含的物种，我们首先将目标基因集中的序列应用blastx 比对到string 数据库中包含的参考物种的蛋白质序列上，并利用比对上的该参考物种的蛋白质互作关系构建互作网络。

我们提供差异基因蛋白互作网络数据文件，此文件可以直接导入Cytoscape 软件进行可视化编辑。Cytoscape 软件使用方法可参考我们提供的使用说明文档

(CytoscapeQuickStart.pdf)。客户可以针对一些网络的拓扑属性进行统计和标示作图，比如：互作网络图中节点(node)的大小与此节点的度(degree)成正比，即与此节点相连的边越多，它的度越大，节点也就越大，这些节点在网络中可能处于较为核心的位置。节点的颜色与此节点的聚集系数(clustering coefficient)相关，颜色梯度由绿到红对应聚集系数的值由低到高；聚集系数表示此节点的邻接点之间的连通性好坏，聚集系数值越高表示此节点的邻接点之间的连通性越好等等。根据不同的研究目的和需求，客户还可以在网络图中进行调整节点位置和颜色、标注表达量水平等操作。需要注意的是，通过blast 比对得到的结果不能保证较好的准确性，这部分的工作仅供参考，辅助客户发现一些重要的候选基因。按我们提供的使用说明将文件导入Cytoscape

软件后的效果图如下：

图12　Cytoscape 软件界面

12.2 差异基因蛋白互作网络分析Q&A

问：PPI 分析是用的什么软件？用的什么数据库？

答：通过blastx 比对，在STRING 数据库中找出这些差异基因间的互作关系，再将得到的互作数据导入Cytoscape 软件实现互作网络的可视化。PPI 分析相关名词的解释:

聚集系数：聚集系数聚集系数:(clustering coefficient)，聚集系数是图中的点倾向于集聚在一起的程度的一种度量。

四、备注

1　文件目录列表

为方便您查看项目结果文件，诺禾致源转录调控团队为您提供了结果文件的目录列表（配有相应的中英文说明），项目结果文件可通过此目录列表查看，点击目录列表即可弹出相应结果文件夹（注注：请保证report 文件夹和results 文件夹在同一文件夹下）。

文件目录列表：html

2　软件列表

有参转录组软件介绍

Analysis

Mapping 基因表达水平分析可变剪切预测新转录本预测SNP detection

Software

Tophat HTSeq cufflinks ASprofile cufflinks GATK2DEGSeq

Version

v2.0.12v0.6.12.1.11.02.1.1v3.21.12.01.10.13.0.8Release2.12

v2.0

Parameters

mismatch = 2-m union默认参数默认参数默认参数

QUAL 1 &&

qvalue

padj

Corrected P-Value

Remarks

与参考基因组进行比对

差异表达分析

DESeq edgeR

对于有重复的样品使用DESeq ，无重复的样品使用DEGSeq ，特殊情况下使用

edgeR 通过 hmmscan 得到新基因的

GO 注释文件

GO 富集KEGG 富集

GOSeq ，topGO,hmmscan

KOBAS

蛋白互作分析BLAST v2.2.28

e-value = 1e-10 && string

score >700

若物种存在于数据库string 中，则直接提取相应的互作信息；若无，则提取近缘物种的互作

信息

3　Methods 英文版

Methods 英文版：PDF

4　Novofinder 软件使用说明PDF 版

为了便于大家查询、整合RNA-seq 分析中产生的多种多样的数据表格，我们新研发推出了试行版NovoFinder 软件，欢迎大家使用！NovoFinder 使用说明PDF 版：PDF

5　结果文件

结果文件解析说明：

*.tar.gz形式的压缩文件：

Unix/Linux/Mac用户Windows 用户Unix/Linux/Mac用户Windows 用户Unix/Linux/Mac用户Windows 用户

使用tar -zxvf *.tar.gz命令

使用解压缩软件如WinRAR 、7-Zip 等使用命令gzip –d *.gz命令

使用解压缩软件如WinRAR 、7-Zip 等使用命令unzip *.zip命令

使用解压缩软件如WinRAR 、7-Zip 等

*.gz形式的压缩文件：

*.zip形式的压缩文件：

更多解压命令可参见网络资料：http://www.php100.com/html/webkaifa/Linux/2009/1213/3652.html。

2. 结果文件查看说明:

*.fasta

序列文件，fasta 格式，一般为基因序列或者基因组序列。因文件一般较大，打开较为困难，为您提供了sampled*.fasta文件（部分*.fasta中的序列），方便您查看文件格式。

unix/Linux/Mac用户使用 less 或 more 命令查看；windows 使用高级文本编辑器Editplus/Notepad++等查看。

unix/Linux/Mac用户使用 less 或 more 命令查看；windows 用户使用高级文本编辑器Editplus/Notepad++等查看。

unix/Linux/Mac用户使用 less 或 more 命令查看；

*.xls,*.txt结果数据表格文件，文件以制表符（Tab ）分隔。

windows 用户使用高级文本编辑器

Editplus/Notepad++ 等，也可以用Microsoft Excel打开。

unix/Linux/Mac用户使用display 命令打开。

*.png

结果图像文件，位图，无损压缩。

windows 用户可以使用图片浏览器打开，如photoshop 等。

unix/Linux用户使用evince 命令打开。

windows/Mac用户可以使用Adobe Reader/福昕阅读器/网页浏览器等打开。

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

序列文件，fastq 格式，一般为reads 序列；因文件一般较大，打开较为困难。为您提供了

*.fq/fastq

sampled*.fastq文件（部分*.fastq中的序列），方便您查看文件格式。

*.pdf

结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用Adobe Illustrator 进行图片编辑，用于文章发表等。

五、参考文献

Marioni, J. C., C. E. Mason, et al. (2008). RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome research.Mortazavi, A., B. A. Williams, et al. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature methods.Wang, Z., M. Gerstein, et al. (2009). RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews Genetics.

Langmead, B., Trapnell, C., Pop, M. & Salzberg, S.L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol.(Bowtie)Langmead, B. and S. L. Salzberg (2012). Fast gapped-read alignment with Bowtie 2. Nature methods.(Bowtie 2)

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

Kim, D., G. Pertea, et al. (2012).TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.(TopHat2)

Marquez Y, Brown JW, Simpson C, Barta A, Kalyna M. (2012). Transcriptome survey reveals increased complexity of the alternative splicing landscape in Arabidopsis.Anders, S.(2010). HTSeq: Analysing high-throughput sequencing data with Python.(HTSeq)

Trapnell, C., A. Roberts, et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. nature protocols.(Tophat & Cufflinks)

Trapnell, C. et al. (2010).Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat.Biotechnol.(Cufflinks)

McKenna, A, Hanna, M, Banks, E, Sivachenko, A, Cibulskis, K, Kernytsky, A, Garimella, K, Altshuler, D, Gabriel, S, Daly, M, DePristo, MA. 2010. The Genome Analysis Toolkit: aMapReduce framework for analyzing next-generation DNA sequencing data. Genome Research.(GATK)

Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol.(DESeq)Anders, S. and Huber, W. (2012). Differential expression of RNA-Seq data at the gene level-the DESeq package.(DESeq)

Wang, L.Feng, Z.Wang, X.Zhang, X. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics.(DEGseq)

Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics.(edgeR)Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010).Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology.(GOseq)Kanehisa, M., M. Araki, et al. (2008). KEGG for linking genomes to life and the environment. Nucleic acids research.(KEGG)

Mao, X., Cai, T., Olyarchuk, J.G., Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary.Bioinformatics.(KOBAS)

一、建库测序流程

保证测序数据的准确性、可靠性，诺禾致源对样品检测、建库、测序每一个生产步骤都严格把控，从根本上确保了高质量数据的产出。流程图如下：

1 Total RNA样品检测

2 文库构建

文库。构建原理图如下：

3 库检

4 上机测序

库检合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling 后进行HiSeq/MiSeq测序。

二、生物信息分析流程

获得原始测序序列(Sequenced Reads)

后，在有相关物种参考序列或参考基因组的情况下，通过如下流程进行生物信息分析：

三、结果展示及说明

1　原始序列数据

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

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT +

#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

其中第一行以“@”开头，随后为Illumina 测序标识符(Sequence Identifiers)和描述文字(选择性部分) ；第二行是碱基序列；

第三行以“+”开头，随后为Illumina 测序标识符(选择性部分) ；

第四行是对应碱基的测序质量，该行中每个字符对应的 ASCII 值减去 33，即为对应第二行碱基的测序质量值。Illumina 测序标识符详细信息如下：HWI-ST127671

C1162ACXX [1**********]581N 0CGATGT

Instrument – unique identifier of the sequencerrun number – Run number on instrumentFlowCell ID – ID of flowcellLaneNumber – positive integerTileNumber – positive integer

X – x coordinate of the spot. Integer which can be negativeY – y coordinate of the spot. Integer which can be negativeReadNumber - 1 for single reads; 1 or 2 for paired ends

　测序数据质量评估

2.1　测序错误率分布检查

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

Phred 分值

10203040

不正确的碱基识别

1/101/1001/10001/10000

碱基正确识别率

90%99%99.9%99.99%

Q-sorce

Q10Q20Q30Q40

测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于RNA-seq 技术，测序错误率分布具有两个特点：

图1　测序错误率分布图

横坐标为reads 的碱基位置，纵坐标为单碱基错误率

2.2

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

GC 含量分布检查用于检测有无AT 、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。

图2　GC 含量分布图

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

2.4　测序数据质量情况汇总

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

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

Sample name

sampleA1_1sampleA1_2sampleA2_1sampleA2_2sampleB1_1sampleB1_2sampleB2_1sampleB2_2sampleC1_1sampleC1_2sampleC2_1sampleC2_2

数据质量情况详细内容如下：

Raw reads

[***********][***********][***********][***********][***********]521158

Clean reads

[***********][***********][***********][***********][***********]035877

Clean bases

4.19G 4.19G 4.14G 4.14G 3.92G 3.92G 3.65G 3.65G 3.55G 3.55G 3.5G 3.5G

Error rate(%)

0.030.040.030.040.030.040.030.040.030.040.030.04

Q20(%)

95.6194.1395.8893.9895.9094.2395.9294.0595.9594.3995.9494.28

Q30(%)

91.3289.0791.8788.8491.8889.2491.9288.9491.9889.5092.0089.38

GC content(%)

48.4348.4748.4248.4449.3049.3048.8748.8848.5748.5747.7247.72

(6) Q20、Q30：分别计算 Phred 数值大于20、30的碱基占总体碱基的百分比。(7) GC content：计算碱基G 和C 的数量总和占总的碱基数量的百分比。

2.5 质量评估Q&A

问：测序错误率会随着测序序列长度的增加而升高，错误率在多少是可以接受的范围？

答：为保证后续分析的质量，诺禾会严格把控cleandata 的筛选标准，具体标准如下： (1) 去除带接头(adapter)的reads ；

(2) 去除N(N表示无法确定碱基信息) 的比例大于10%的reads ；

(3) 去除低质量reads(质量值sQ

clean data/clean reads：对原始数据进行过滤后，剔除了低质量数据的剩余数据。后续分析均基于clean data。

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

3　参考序列比对分析

的算法示意图：

3.2

　Reads 在参考基因组不同区域的分布情况

将比对到基因组上的reads 分布情况进行统计，定位区域分为Exon(外显子) 、Intron(内含子) 和Intergenic(基因间区) 。

图3.2　Reads 在参考基因组不同区域的分布情况

3.3

　Reads 在染色体上的密度分布情况

图3.3　Reads 在染色体上的密度分布图

Mb) ，纵坐标为mapped 到染色体上的reads 数(单位为M) ，图中灰色区域表示95%的置信区间

3.4　Reads 比对结果可视化

。

图3.4　IGV 浏览器界面

3.5 参考序列比对分析 Q&A问：有参分析都需要什么文件？

问：mapping 时用的是read 全长，还是头尾有处理？

4　可变剪切分析

事件分类如下图所示：

12类可变剪切事件定义如下:

图4.1　AS 分类和数量统计

纵轴为可变剪切事件的分类缩写，横轴为该种事件下可变剪切的数量，不同样品用不同子图和颜色区分

表4.2　AS 结构和表达量统计

event_idevent_typegene_idchromevent_startevent_end

[***********]0031000004

TSS TTS SKIP_ONSKIP_OFF

CUFF.1000CUFF.1000CUFF.1000CUFF.1000

1111

[***********][***********]

event_pattern

[***********]130657314,130657463-130657479,[***********],130661448

strand

----

fpkm ref_id

0.0000000000ENSMUST[1**********]0.0000000000ENSMUST[1**********]0.0000000000ENSMUST[1**********]0.0000000000ENSMUST[1**********]

(1) event_id: AS事件编号

(5) event_start: AS事件起始位置(6) event_end: AS事件结束位置

(8) strand: 基因正负链信息

(9) fpkm: 此AS 类型该基因表达量

(10) ref_id: 此基因在参考注释文件中的编号

4.3 可变剪切分析 Q&A

问：什么是可变剪切(alternative splicing)？

5　新转录本预测

5.1　新转录本预测

表5.1　新转录本结构注释结果

seqname

1111

source

novelGene novelGene novelGene novelGene

feature

exon exon exon exon

start

[***********][1**********]178

end

[***********][1**********]237

score strand frame

. . . .

++++

. . . .

attributes

(1) seqname：染色体编号

(9) attributes：属性，包括基因编号、转录本编号等信息

5.2　基因结构优化

表5.2　已知基因结构优化

Gene_id

Chromosome

ENSMUSG[1**********]3ENSMUSG[1**********]X ENSMUSG[1**********]11ENSMUSG[1**********]

(1) Gene_id：原注释文件中基因命名编号(2) Chromosome：染色体编号(3) Strand：正负链信息

(4) Original_span：原注释文件中基因起始位置~终止位置

(5) Assembled_span：转录组拼接结果中基因起始位置~终止位置

Strand

Original_span

Assembled_span

-108107280~[***********]~108154921+161117193~[***********]~161258395+121237253~[***********]~121258656+

17281185~17289115

17280908~17289115

5.3 新转录本预测 Q&A问：新转录本预测的意义？

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

6.2 SNP和InDel Q&A

问：SNP 分析的参考序列是什么，即REF 是指什么？

答：以0,12为例，表示有0个reads 支持REF 的碱基，即没有支持该碱基的reads ，12个reads 支持ALT 的碱基。问：SNP 具体指什么？其与InDel 有何区别？

答：一般而言，SNP 是指变异频率大于1%的单核苷酸变异。InDel(insertion-deletion)则是插入或者缺失，insert 或者deletion 。问：SNP 和InDel 相关名词的解释

SNP ：SNP （Single Nucleotide Polymorphisms）单核苷酸多态性。

7.2 基因表达水平分析Q&A问：基因表达水平如何计算？

答：有参转录组当中，认为FPKM>1是基因表达的。这个阈值是主流杂志推荐的，也能够很好的反应基因的表达水平。基因表达水平分析相关名词的解释:

FPKM ：expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced，是每百万fragments 中来自某一基因每千碱基长度的fragments 数目。

RPKM ：expected number of reads Per Kilobase of transcript sequence per Millions base pairs sequenced，是每百万reads 中来自某一基因每千碱基长度的reads 数目。

　RNA-seq 整体质量评估

8.1 RNA-Seq相关性检查

图8.1　RNA-Seq 相关性检查

热图：样品间相关系数热图；散点图(若样品多于4组，则仅展示生物学重复之间的散点图) ：样品间的相关系数散点图，R 2:pearson相关系数的平方。

8.2 RNA-seq整体质量评估Q&A

问：样品间的相关性有何意义？如何计算?

差异表达分析

9.1　基因表达水平对比

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

布图，横坐标为log 10(FPKM+1), 纵坐标为基因的密度

9.2　差异表达基因列表

K ij ～ NB(µij , σij 2)

对于无生物学重复的样品，先采用TMM 对read count数据进行标准化处理，之后用DEGseq 进行差异分析。差异表达基因列表如下：

表9.2　差异基因列表

Gene Id

ENSMUSG[1**********]ENSMUSG[1**********]ENSMUSG[1**********]ENSMUSG[1**********]

差异基因列表主要包括的内容：(1) Gene id: 基因编号

(2) Sample1：校正后样品1的readcount 值(3) Sample2：校正后样品2的readcount 值(4) log2FoldChange: log2(Sample1/Sample2)

sampleA

27459.[**************].[1**********]97.[**************]22.5175799308

sampleB

3588.[**************]00.[1**********]33.[**************].[1**********]

log2FoldChange

2.936-2.789-2.86772.3387

pval

1.5913e-1709.2315e-1086.3681e-1089.3935e-106

p-adjusted

3.2181e-1666.223e-1046.223e-1044.7491e-102

(5) pvalue(pval): 统计学差异显著性检验指标

(6) qvalue(p-adjusted): 校正后的pvalue 。qvalue 越小，表示基因表达差异越显著

9.3

　差异表达基因筛选

图9.3　差异基因火山图

的统计学显著性

9.4

　差异基因聚类分析

图9.4　差异基因聚类图

表达量

9.5

　差异基因维恩图

差异基因维恩图展示了各比较组间差异基因的个数，以及比较组间的重叠关系。(仅提供两组、三组和四组比较的venn 图)

图9.5　差异基因维恩图

每个大圆圈中的数字之和代表该比较组合的差异基因总个数，圆圈交叠的部分表示组合之间共有的差异基因。

9.6 差异表达分析Q&A

问：如何判断一个基因是否是差异基因？如果是差异基因，如何判断该基因的表达量是上调还是下调？

答：无生物学重复时，先采用TMM 对read count数据进行标准化处理，再用DEGseq 进行差异分析, 筛选阈值为qvalue

差异分析。

问：如何判断差异基因在两个样品间的差异大小？

问：某基因在两个样本中表达量差别很大，却不存在与显著差异的基因列表中，这是为何？

问：差异基因筛选条件最大能设的阈值是多少？很多客户希望通过调整差异基因筛选阈值来找相关基因是否有必要？

答：在DESeq 中，如果某基因的在一个样品中的校正后的readcount 为0，而在另一个样品中不为0，foldchange 会为INF 或者-INF ；如果两个数值均为

答：h_cluster,k-means和som 的聚类图可按照指定顺序绘制，但聚类热图的顺序不能调整。问：聚类分析是怎么做的？

答：校正后的p 值（padj/qvalue），是对p 值进行了多重假设检验，校正后的p 值比原有的p 值要大，能更好地控制假阳性率。

10　差异基因GO 富集分析

表10.1　差异基因GO 富集列表

GO accession

GO:0005488GO:0036094GO:0000166GO:1901265

Description

binding

small molecule bindingnucleotide bindingnucleoside phosphate binding

Term type

molecular_functionmolecular_functionmolecular_functionmolecular_function

Over represented p-Value

4.847e-071.7728e-052.9533e-052.9533e-05

Corrected p-ValueDEG itemDEG list

0.00197760.0210210.0210210.021021

[1**********]63

[**************]9

结果表格详细内容如下：

(1) GO accession：Gene Ontology数据库中唯一的标号信息(2) Description：Gene Ontology功能的描述信息

(5) Corrected p-Value：矫正后的P-Value ，一般情况下，Corrected_pValue

10.2

　差异基因GO 富集柱状图

图10.2　GO 富集柱状图

按生物过程、细胞组分和分子功能三大类别及差异基因上下调分类画的三个子图

10.3

　差异基因GO 富集DAG 图

图10.3　GO 富集有向无环图

p-value

10.4 差异基因GO 富集分析Q&A问：GO 富集分析所使用的软件是什么？

答：GO 富集分析均采用R 包，富集采用的为GOseq, 有向无环图采用的为topGO 。问：GO 富集分析一般分析到二级，还能继续分析吗？

答：GO 富集分析是针对所有注释的GO 进行统计检验，任何等级的都有。问：GO 富集和GO 分类有何区别？

11　差异基因KEGG 富集分析

表11.1　差异基因KEGG 富集列表

#Term

Axon guidanceGlutamatergic synapseRetrograde endocannabinoid

signaling

Circadian entrainment结果表格详细内容如下：

Database

KEGG PATHWAY KEGG PATHWAY KEGG PATHWAY KEGG PATHWAY

mmu04360mmu04724mmu04723mmu04713

Sample number

53423735

Background number

[1**********]

P-Value

2.[1**********]e-060.[**************]0.[**************]0.[1**********]844

Corrected P-Value

0.[**************]0.[1**********]160.[1**********]660.[1**********]66

(1) #Term：KEGG 通路的描述信息。(2) Database：KEGG 数据库。

(3) ID：KEGG 数据库中通路唯一的编号信息。(4) Sample number：该通路下差异基因的个数。

(5) Background number：该通路下注释基因的个数。(6) P-value：富集分析统计学显著水平。

(7) Corrected P-value：矫正后的统计学显著水平，Corrected P-value

11.2

　差异基因KEGG 富集散点图

图11.2　差异基因KEGG 富集散点图

纵轴表示pathway 名称，横轴表示Rich factor，点的大小表示此pathway 中差异表达基因个数多少，而点的颜色对应于不同的Qvalue 范围

11.3

富集KEGG 通路图

为便于查看差异基因在通路图中的分布情况，我们将差异基因标注到通路图中，查看方法如下：拿到全部分析结果后，打开results 结果目录中

图11.3 显著富集的KEGG pathway代谢通路图

11.4 差异基因KEGG 富集分析Q&A

问：为什么编码同一个酶的基因，会有的上调有的下调？

问：KEGG Pathway图片中节点及底色有何含义？

答：KEGG Pathway共有5种类型, 分别为： map: Reference pathway ko: Reference pathway (KO) ec: Reference pathway (EC)

rn: Reference pathway (Reaction) org: Organism-specific pathway map 图中节点及底色的含义如下：

1、"map" pathway: 节点代表某一基因、该基因编码的酶及这个酶参与的反应，将鼠标悬停在某个节点上，可以看到上述信息，底色以无色表示。

3、"org" pathway：物种特异性通路图。节点含义同map pathway。物种相关的节点底色以绿色表示。

每个数据路的pathway 都有相应的唯一编号，如map00010，据此可在kegg 数据库官网查询。有参考基因组项目中得到的通路图都是Organism-specific pathway 。

12　蛋白互作网络分析

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

我们提供差异基因蛋白互作网络数据文件，此文件可以直接导入Cytoscape 软件进行可视化编辑。Cytoscape 软件使用方法可参考我们提供的使用说明文档

软件后的效果图如下：

图12　Cytoscape 软件界面

12.2 差异基因蛋白互作网络分析Q&A

问：PPI 分析是用的什么软件？用的什么数据库？

聚集系数：聚集系数聚集系数:(clustering coefficient)，聚集系数是图中的点倾向于集聚在一起的程度的一种度量。

四、备注

1　文件目录列表

文件目录列表：html

2　软件列表

有参转录组软件介绍

Analysis

Mapping 基因表达水平分析可变剪切预测新转录本预测SNP detection

Software

Tophat HTSeq cufflinks ASprofile cufflinks GATK2DEGSeq

Version

v2.0.12v0.6.12.1.11.02.1.1v3.21.12.01.10.13.0.8Release2.12

v2.0

Parameters

mismatch = 2-m union默认参数默认参数默认参数

QUAL 1 &&

qvalue

padj

Corrected P-Value

Remarks

与参考基因组进行比对

差异表达分析

DESeq edgeR

对于有重复的样品使用DESeq ，无重复的样品使用DEGSeq ，特殊情况下使用

edgeR 通过 hmmscan 得到新基因的

GO 注释文件

GO 富集KEGG 富集

GOSeq ，topGO,hmmscan

KOBAS

蛋白互作分析BLAST v2.2.28

e-value = 1e-10 && string

score >700

若物种存在于数据库string 中，则直接提取相应的互作信息；若无，则提取近缘物种的互作

信息

3　Methods 英文版

Methods 英文版：PDF

4　Novofinder 软件使用说明PDF 版

为了便于大家查询、整合RNA-seq 分析中产生的多种多样的数据表格，我们新研发推出了试行版NovoFinder 软件，欢迎大家使用！NovoFinder 使用说明PDF 版：PDF

5　结果文件

结果文件解析说明：

*.tar.gz形式的压缩文件：

Unix/Linux/Mac用户Windows 用户Unix/Linux/Mac用户Windows 用户Unix/Linux/Mac用户Windows 用户

使用tar -zxvf *.tar.gz命令

使用解压缩软件如WinRAR 、7-Zip 等使用命令gzip –d *.gz命令

使用解压缩软件如WinRAR 、7-Zip 等使用命令unzip *.zip命令

使用解压缩软件如WinRAR 、7-Zip 等

*.gz形式的压缩文件：

*.zip形式的压缩文件：

更多解压命令可参见网络资料：http://www.php100.com/html/webkaifa/Linux/2009/1213/3652.html。

2. 结果文件查看说明:

*.fasta

unix/Linux/Mac用户使用 less 或 more 命令查看；windows 使用高级文本编辑器Editplus/Notepad++等查看。

unix/Linux/Mac用户使用 less 或 more 命令查看；windows 用户使用高级文本编辑器Editplus/Notepad++等查看。

unix/Linux/Mac用户使用 less 或 more 命令查看；

*.xls,*.txt结果数据表格文件，文件以制表符（Tab ）分隔。

windows 用户使用高级文本编辑器

Editplus/Notepad++ 等，也可以用Microsoft Excel打开。

unix/Linux/Mac用户使用display 命令打开。

*.png

结果图像文件，位图，无损压缩。

windows 用户可以使用图片浏览器打开，如photoshop 等。

unix/Linux用户使用evince 命令打开。

windows/Mac用户可以使用Adobe Reader/福昕阅读器/网页浏览器等打开。

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

序列文件，fastq 格式，一般为reads 序列；因文件一般较大，打开较为困难。为您提供了

*.fq/fastq

sampled*.fastq文件（部分*.fastq中的序列），方便您查看文件格式。

*.pdf

结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用Adobe Illustrator 进行图片编辑，用于文章发表等。

五、参考文献

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

Kim, D., G. Pertea, et al. (2012).TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.(TopHat2)

Trapnell, C., A. Roberts, et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. nature protocols.(Tophat & Cufflinks)

Trapnell, C. et al. (2010).Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat.Biotechnol.(Cufflinks)

Wang, L.Feng, Z.Wang, X.Zhang, X. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics.(DEGseq)

Mao, X., Cai, T., Olyarchuk, J.G., Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary.Bioinformatics.(KOBAS)

诺禾致源有参转录组分析流程

相关文章