目录 一 建库测序流程... 2 1.1 Total RNA 样品检测... 2 1.2 文库构建... 2 1.3 库检... 2 1.4 上机测序... 2 二 生物信息分析流程... 3 三 结果展示及说明... 4 3.1 原始序列数据... 4 3.2 测序数据质量评估... 4 3.2.1 测序错误率分布检查... 4 3.2.2 A/T/G/C 含量分布检查... 5 3.2.3 测序数据过滤... 5 3.2.4 测序数据质量情况汇总... 5 3.2.5 质量评估 Q&A... 6 3.3 参考序列比对分析... 6 3.3.1 Reads 与参考基因组比对情况统计... 6 3.3.2 Reads 在参考基因组不同区域的分布情况... 6 3.3.3 Reads 在染色体上的密度分布情况... 7 3.3.4 Reads 比对结果可视化... 7 3.3.5 参考序列比对分析 Q&A... 7 3.4 基因表达水平分析... 8 3.4.1 基因表达水平分析... 8 3.4.2 基因表达水平对比... 8 3.4.3 基因表达水平分析 Q&A... 8 3.5 RNA-seq 整体质量评估... 9 3.5.1 RNA-Seq 相关性检查... 9 3.5.2 RNA-seq 整体质量评估 Q&A... 9 3.6 基因差异表达分析... 9 3.6.1 差异表达基因列表... 9 3.6.2 差异表达基因筛选... 10 3.6.3 差异基因聚类分析... 10 3.6.4 基因表达维恩图... 40 3.6.5 差异表达分析 Q&A... 40 3.7 差异基因 GO 富集分析... 41 3.7.1 差异基因 GO 富集列表... 41 3.7.2 差异基因 GO 富集柱状图... 41 3.7.3 差异基因 GO 富集 DAG 图... 42 3.7.4 差异基因 GO 富集分析 Q&A... 42 3.8 差异基因 KEGG 富集分析... 42 3.8.1 差异基因 KEGG 富集列表... 42 3.8.2 差异基因 KEGG 富集散点图... 43 3.8.3 富集 KEGG 通路图... 43 3.8.4 差异基因 KEGG 富集分析 Q&A... 43 3.9 SNP 和 InDel 分析... 44 3.9.1 SNP 和 InDel 分析... 44 3.9.2 SNP 和 InDel Q&A... 44 3.10 新转录本预测... 44 3.11 基因结构分析... 45 3.11.1 操纵子预测... 45 3.11.2 TSS 和 TTS 预测... 45 3.11.3 启动子预测... 45 3.12 UTR 分析... 55 3.12.1 UTR 预测和长度分布... 55 3.12.2 5 UTR SD 序列预测... 55 3.12.3 3 UTR 不依赖 ρ 因子的终止子预测... 55 Gene_ID Term_start Term_end Strand 5'_tail 5'_stem Loop 3'_stem 3'_tail... 55 3.13 反义转录本预测... 55 3.14 srna 分析... 56 3.14.1 srna 预测和长度分布... 56 3.14.2 srna 二级结构预测... 56 3.14.3 srna 靶基因预测... 56 3.14.4 srna 表达水平分析... 56 四 参考文献... 57
一 建库测序流程 从 RNA 样品到最终数据获得, 样品检测 建库 测序每一个环节都会对数据质量和数量产生影响, 而数据质量又会直接影响后续信息分析的结果 为了从源头上保证测序数据的准确性 可靠性, 必须要对样品检测 建库 测序每一个生产步骤都严格把控, 从根本上确保了高质量数据的产出 流程图如下 : 1.1 Total RNA 样品检测 对 RNA 样品的检测主要包括 4 种方法 : (1) 琼脂糖凝胶电泳分析 RNA 降解程度以及是否有污染 (2) Nanodrop 检测 RNA 的纯度 (OD 260/280 比值 ) (3) Qubit 对 RNA 浓度进行精确定量 (4) Agilent 2100 精确检测 RNA 的完整性 1.2 文库构建 样品检测合格后, 通过 Ribo-zero 试剂盒去除 rrna 来富集 mrna 随后加入 fragmentation buffer 将 mrna 打断成短片段, 以 mrna 为模板, 用六碱基随机引物 (random hexamers) 合成一链 cdna, 然后加入缓冲液 dntps(dntp 中的 dttp 用 dutp 取代 ) 和 DNA polymerase I 和 RNase H 合成二链 cdna, 再用 AMPure XP beads 纯化双链 cdna, 之后用 USER 酶降解含有 U 的 cdna 第二链 纯化的双链 cdna 先进行末端修复 加 A 尾并连接测序接头, 再用 AMPure XP beads 进行片段大小选择 最后进行 PCR 扩增, 并用 AMPure XP beads 纯化 PCR 产物, 得到最终的文库 构建原理图如下 : 1.3 库检 文库构建完成后, 先使用 Qubit2.0 进行初步定量, 稀释文库至 1ng/ul, 随后使用 Agilent 2100 对文库的插入片段长度 (insert size) 进行检测,insert size 符合预期后, 使用 Q-PCR 方法对文库的 有效浓度进行准确定量 ( 文库有效浓度 >2nM), 以保证文库质量 1.4 上机测序 库检合格后, 把不同文库按照有效浓度及目标下机数据量的需求 pooling 后进行 HiSeq/MiSeq 测序
二 生物信息分析流程 获得原始测序序列 (Sequenced Reads) 后, 在有相关物种参考序列或参考基因组的情况下, 通过如下流程进行生物信息分析 :
三 结果展示及说明 3.1 原始序列数据 高通量测序 ( 如 Illumina HiSeq TM 2500/Miseq TM ) 得到的原始图像数据文件经 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-ST1276 Instrument unique identifier of the sequencer 71 run number Run number on instrument C1162ACXX FlowCell ID ID of flowcell 1 LaneNumber positive integer 1101 TileNumber positive integer 1208 X x coordinate of the spot. Integer which can be negative 2458 Y y coordinate of the spot. Integer which can be negative 1 ReadNumber - 1 for single reads; 1 or 2 for paired ends N whether it is filtered - NB:Y if the read is filtered out, not in the delivered fastq file, N otherwise 0 control number - 0 when none of the control bits are on, otherwise it is an even number CGATGT Illumina index sequences 3.2 测序数据质量评估 3.2.1 测序错误率分布检查 每个碱基测序错误率是通过测序 Phred 数值 (Phred score, Q phred ) 通过公式 ( 公式 1 :Q phred = -10log10(e)) 转化得到, 而 Phred 数值是在碱基识别 (Base Calling) 过程中通过一种概率模型计算得到, 这种模型可以准确地预测碱基判别的错误率 Phred 分值, 不正确的碱基识别率, 碱基正确识别率以及 Q-score 的对应关系如下表所显示 : Illumina Casava 1.8 版本碱基识别与 Phred 分值之间的简明对应关系 : Phred 分值 不正确的碱基识别 碱基正确识别率 Q-sorce 10 1/10 90% Q10 20 1/100 99% Q20 30 1/1000 99.9% Q30 40 1/10000 99.99% Q40 测序错误率与碱基质量有关, 受测序仪本身 测序试剂 样品等多个因素共同影响 对于 RNA-seq 技术, 测序错误率分布具有两个特点 : (1) 测序错误率随着测序序列 (Sequenced Reads) 长度的增加而升高 这是由测序过程中化学试剂的消耗导致的, 为 Illumina 高通量测序平台所具有的特征 (2) 前 6 个碱基具有较高的测序错误率, 此长度恰好为 RNA-seq 建库过程中反转录所需的随机引物长度 前 6 个碱基测序错误率较高是因为随机引物和 RNA 模版的不完全结合 (Jiang et al.) 图 1 测序错误率分布图横坐标为 reads 的碱基位置, 纵坐标为单碱基错误率
3.2.2 A/T/G/C 含量分布检查 GC 含量分布检查用于检测有无 AT GC 分离现象, 而这种现象可能是测序或者建库所带来的, 并且会影响后续的定量分析 在 Illumina 测序平台的转录组测序中, 反转录成 cdna 时所用的 6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性 而这种偏好偏好性与测序的物种和实验室环境无关, 但会影响转录组测序的均一性程度 (Hansen et al.) 此外, 理论上普通文库的 G 和 C 碱基及 A 和 T 碱基含量每个测序循环上应分别相等, 且整个测序过程稳定不变, 呈水平线 ; 而对于链特异性建 库则会出现 GC 分离的现象 对于 Illumina 测序来说, 由于随机引物扩增偏差等原因, 测序所得每个 read 前 6-7 个碱基出现较大波动属正常情况 3.2.3 测序数据过滤 图 2 GC 含量分布图横坐标为 reads 的碱基位置, 纵坐标为单碱基所占的比例 ; 不同颜色代表不同的碱基类型 测序得到的原始测序序列, 里面含有带接头的 低质量的 reads, 为了保证信息分析质量, 必须对 raw reads 进行过滤, 得到 clean reads, 后续分析都基于 clean reads 数据处理的步骤如下 : (1) 去除带接头 (adapter) 的 reads; (2) 去除 N(N 表示无法确定碱基信息 ) 的比例大于 10% 的 reads; (3) 去除低质量 reads( 质量值 Q pred <= 5 的碱基数占整个 read 长度的 50% 以上的 reads) RNA-seq 的接头 (Adapter, Oligonucleotide sequences for TruSeq TM RNA and DNA Sample Prep Kits) 信息 : RNA 5 Adapter (RA5), part # 15013205: 5 -AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-3 RNA 3 Adapter (RA3), part # 15013207: 5 -GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCTCGTATGCCGTCTTCTGCTTG-3 图 3 原始数据组成 3.2.4 测序数据质量情况汇总 样品测序产出数据质量评估情况详见表 1
表 1 数据产出质量情况一览表 Sample name Raw reads Clean reads Clean bases Error rate(%) Q20(%) Q30(%) GC content(%) fjoligo 10542762 10217260 1.53G 0.02 96.73 91.34 51.25 fjh2o 10232386 9908402 1.49G 0.02 96.43 90.67 51.53 3.2.5 质量评估 Q&A 问 : 测序错误率会随着测序序列长度的增加而升高, 错误率在多少是可以接受的范围? 答 : 测序会进行严格的数据质量把控 一般情况下, 单个碱基位置的测序错误率应该低于 1%, 最高在 6% 左右可以接受 问 : 质控的标准是什么? 是否严格? 答 : 为保证后续分析的质量, 基因帮会严格把控 cleandata 的筛选标准, 具体标准如下 : 去除带接头 (adapter) 的 reads; 去除 N(N 表示无法确定碱基信息 ) 的比例大于 10% 的 reads; 去除低质量 reads( 质量值 Q pred <= 20 的碱基数占整个 read 的 50% 以上的 reads) 相关名词解释 : adapter: 接头, 用于上机测序 建库时引入的接头序列与测序芯片 (flow cell) 上固定的接头相互识别 index: 测序的标签, 用于测定混合样本, 通过每个样本添加的不同标签进行数据区分, 鉴别测序样品 Q20,Q30:Phred 数值大于 20 30 的碱基占总体碱基的百分比, 其中 Phred=-10log10(e). raw data/raw reads: 测序下机的原始数据 clean data/clean reads: 对原始数据进行过滤后, 剔除了低质量数据的剩余数据 后续分析均基于 clean data 3.3 参考序列比对分析 针对细菌等基因密度较高的生物, 我们用 Bowtie2 将过滤后的测序序列进行基因组定位分析 如果参考基因组选择合适, 而且相关实验不存在污染, 实验所产生的测序序列的定位的百分 比正常情况下会高于 70% (Total Mapped Reads or Fragments), 其中具有多个定位的测序序列 (Multiple Mapped Reads or Fragments) 占总体的百分比通常不会超过 10% 3.3.1 Reads 与参考基因组比对情况统计 表 2 Reads 与参考基因组比对情况一览表 Sample_name fjoligo fjh2o Total reads 10217260 9908402 Total mapped 10086456 (98.72%) 9813087 (99.04%) Multiple mapped 403275 (3.95%) 565005 (5.7%) Uniquely mapped 9683181 (94.77%) 9248082 (93.34%) Read-1 4838515 (47.36%) 4626263 (46.69%) Read-2 4844666 (47.42%) 4621819 (46.65%) Reads map to '+' 4835963 (47.33%) 4623616 (46.66%) Reads map to '-' 4847218 (47.44%) 4624466 (46.67%) 3.3.2 Reads 在参考基因组不同区域的分布情况 根据基因组的基因注释信息, 对 Total mapped reads 比对到基因组上的各个部分的情况进行统计 正常情况下, 基因区的 reads 定位的百分比含量应该较高, 而定位到基因间区的测序序列可 能是因为基因组注释不完全或是非编码 RNA 或是背景噪声 图 4 Reads 在参考基因组不同区域的分布情况
3.3.3 Reads 在染色体上的密度分布情况 对 Total mapped reads 的比对到基因组 ( 分正负链 ) 的密度进行统计, 下图是随机抽取部分 reads, 展示其在染色体上的 mapping 情况 图 5 Reads 在染色体上的密度分布图最外圈是基因组 ; 中间的灰色背景区是实际的抽取的 reads 的分布情况, 红色 mapping 到正链, 蓝色到负链 ; 最里面的圆圈, 橘黄色为正链 coverage 分布, 绿色为负链 coverage 分布, 超过所有 coverage 集均值 +3 倍标准差的奇异点被舍弃 3.3.4 Reads 比对结果可视化 我们提供 RNA-seq Reads 在基因组上比对结果的 bam 格式文件, 部分物种还提供相应的参考基因组和注释文件, 并推荐使用 IGV (Integrative Genomics Viewer) 浏览器对 bam 文件进行可视化浏览 IGV 浏览器具有以下特点 :(1) 能在不同尺度下显示单个或多个读段在基因组上的位置, 包括读段在各个染色体上的分布情况和在注释的外显子 内含子 剪接接合区 基因间区的分布情况等 ;(2) 能在不同尺度下显示不同区域的读段丰度, 以反映不同区域的转录水平 ;(3) 能显示基因及其剪接异构体的注释信息 ;(4) 能显示其他注释信息 ;(5) 既可以从远程服务器端下载各种注释信息, 又可以从本地加载注释信息 3.3.5 参考序列比对分析 Q&A 图 6 IGV 浏览器界面 问 : 有参分析都需要什么文件? 答 : 相应的参考基因组 基因结构注释文件 (gtf/gff/gff3/bed 等格式, 推荐 gtf,gff) 蛋白注释文件(ptt) 及基因的 GO 注释文件的直接下载链接以及基因功能描述文件 问 : 造成 mapping rate 较低的原因可能有哪些? 答 : 当 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 对于序列拼接和定量都大有裨益
3.4 基因表达水平分析 3.4.1 基因表达水平分析 一个基因表达水平的直接体现就是其转录本的丰度情况, 转录本丰度越高, 则基因表达水平越高 在 RNA-seq 分析中, 我们可以通过定位到基因组区域或基因编码区的测序序列 (reads) 的计数来估计基因的表达水平 Reads 计数除了与基因的真实表达水平成正比外, 还与基因的长度和测序深度成正相关 为了使不同基因 不同实验间估计的基因表达水平具有可比性, 人们引入了 FPKM 的概念,FPKM(expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) 是每百万 fragments 中来自某一基因每千碱基长度的 fragments 数目, 其同时考虑了测序深度和基因长度对 fragments 计数的影响, 是目前最为常用的基因表达水平估算方法 (Trapnell, Cole, et al., 2010) 本次采用 HTSeq 软件对各样品进行基因表达水平分析, 使用的模型为 union 结果文件分别统计了不同表达水平下基因的数量以及单个基因的表达水平 一般情况下,FPKM 数值 0.1 或者 1 作为判断基因是否表达的阈值, 不同的文献所采用的阈值不同 表 3 不同表达水平区间的基因数量统计表 FPKM Interval fjoligo fjh2o 0~1 398(9.03%) 400(9.07%) 1~3 52(1.18%) 44(1.00%) 3~15 446(10.12%) 401(9.10%) 15~60 1609(36.50%) 1376(31.22%) >60 1903(43.17%) 2187(49.61%) 表 4 基因表达水平统计表 Gene_id fjoligo fjh2o ECDH10B_3754 495.825010706267 84.0202871877443 ECDH10B_4100 94.7079245431313 103.789926862447 ECDH10B_4501 29.4734205009745 28.958414276974 ECDH10B_1328 247.725964717051 178.67087835058 3.4.2 基因表达水平对比 通过所有基因的 FPKM 分布图以及 violin 图对不同实验条件下的基因表达水平进行比较 对于同一实验条件下的重复样品, 最终的 FPKM 为所有重复数据的平均值 图 7 不同实验条件下基因表达水平比对图 图一 :FPKM 分布图, 横坐标为 log10(fpkm+1), 纵坐标为基因的密度 ; 图二 :FPKM 的 violin 图, 横坐标为样品名称, 纵坐标为 log10(fpkm+1), 每个区域的 violin 图对应五个统计 量 ( 至上而下分别为最大值, 上四分位数, 中值, 下四分位数和最小值 ), 每个 violin 的宽度代表处于该表达量下的基因的多少 3.4.3 基因表达水平分析 Q&A 问 : 基因表达水平如何计算? 答 : 在 RNA-seq 技术中,FPKM(expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) 是每百万 fragments 中来自某一基因每千碱基长度的 fragments 数目, 其同时考虑了测序深度和基因长度对 fragments 计数的影响, 是目前最为常用的基因表达水平估算方法 (Trapnell, Cole, et al., 2010) 问 : 认为基因表达的阈值是多少? 为什么设置为这个阈值? 答 : 有参转录组当中, 认为 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 数目 FPKM 与 RPKM 的用途有一定的相似之处, 都是为了消除技术偏差的表达水平的表示方式 不同的是 FPKM 观察出双端 reads 中 fragment 的差异而 RPKM 关注的是 reads:rather than using read counts, approximates the relative abundance of transcripts in terms of fragments observed from an RNA-Seq experiment, which may not be represented by a single read, such as in paired-end RNA-Seq experiments
3.5 RNA-seq 整体质量评估 3.5.1 RNA-Seq 相关性检查 生物学重复是任何生物学实验所必须的, 高通量测序技术也不例外 (Hansen et al.) 生物学重复主要有两个用途 : 一个是证明所涉及的生物学实验操作是可以重复的且变异不大, 另一个是为了确保后续的差异基因分析得到更可靠的结果 样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标 相关系数越接近 1, 表明样品之间表达模式的相似度越高 Encode 计划建议皮尔逊相关系数的平方 (R 2 ) 大于 0.92( 理想的取样和实验条件下 ) 具体的项目操作中, 我们要求生物学重复样品间 R 2 至少要大于 0.8, 否则需要对样品做出合适的解释, 或者重新进行实验 图 8 RNA-Seq 相关性检查 热图 : 样品间相关系数热图 ; 散点图 ( 若样品多于 4 组, 则仅展示生物学重复之间的散点图 ): 样品间的相关系数散点图,R 2 :pearson 相关系数的平方 3.5.2 RNA-seq 整体质量评估 Q&A 问 : 样品间的相关性有何意义? 如何计算? 答 : 样品间的相关性反应了样品间的相似程度, 即不同处理或组织的样品在表达水平方面的相似度 相关系数越接近 1, 样品间的相似度越高, 样品间的差异基因也越少 生物学重复间 的样品的相关系数应大于生物学重复外的样品的相关系数 相关系数的计算方法有三种 :A. Pearson correlation; B. Spearman rank correlation; C. Kendall s τ 诺禾使用 R 语言进行 Pearson 相关 系数的计算 3.6 基因差异表达分析 3.6.1 差异表达基因列表 基因差异表达分析的输入数据为基因表达水平分析中得到的 readcount 数据, 基因差异表达分析主要分为三步 : 1) 首先对 readcount 进行标准化 (normalization); 2) 然后根据模型进行假设检验概率 (pvalue) 的计算 ; 3) 最后进行多重假设检验校正, 得到 FDR 值 ( 错误发现率 ) 针对不同情况, 会采用不同的软件进行基因差异表达的分析 分析方法如下表 : 类型软件标准化方法 pvalue 计算模型 FDR 计算方法差异基因筛选标准 有生物学重复 DESeq(Anders et al, 2010) DESeq 负二项分布 BH padj < 0.05 无生物学重复 DEGseq(Wang et al, 2010) TMM 泊松分布 BH log 2 (FoldChange) > 1&qvalue < 0.005 第 i 个基因在第 j 个样本中的 read count 值为 Kij, 则有 泊松分布 :Kij ~ P(µij) 负二项分布 :Kij ~ NB(µij,σij2) 表 5 差异基因列表 Gene Id fjoligo fjh2o log2foldchange pval p-adjusted ECDH10B_0002 133.824316207899 315.218529390383-1.236 2.5498e-25 6.9976e-25 ECDH10B_0005 184.873678209619 22.4583176155998 3.0412 2.5945e-28 7.663e-28 ECDH10B_0007 144.207237292994 51.6006583310805 1.4827 1.5102e-08 2.3821e-08 ECDH10B_0010 34.3213224757326 183.676954784727-2.42 2.4548e-31 7.7183e-31 差异基因列表主要包括的内容 : Gene_id: 基因编号 readcount_sample1: 校正后样品 1 的 readcount 值 readcount_sample2: 校正后样品 2 的 readcount 值 log2foldchange: log2(sample1/sample2) pvalue(pval): 统计学差异显著性检验指标 qvalue(padj): 校正后的 pvalue qvalue 越小, 表示基因表达差异越显著
3.6.2 差异表达基因筛选 用火山图可以推断差异基因的整体分布情况 图 9 差异基因火山图有显著性差异表达的基因用红色点 ( 上调 ) 和绿色点 ( 下调 ) 表示, 无显著性差异表达的基因用蓝色点表示 ; 横坐标代表基因在不同样本中表达倍数变化 ; 纵坐标代表基因表达量变化差异的统计学显著性 3.6.3 差异基因聚类分析 聚类分析用于判断差异基因在不同实验条件下的表达模式, 表达模式类似的基因可能具有相似的功能, 或是共同参与同一代谢过程或细胞通路 因此通过将表达模式相同或相近的基因 聚集成类, 可以用来推测未知基因的功能或已知基因的新功能 以不同实验条件下的差异基因的 FPKM 值为表达水平, 做层次聚类 (hierarchical clustering) 分析, 不同颜色的区域代表不同的聚类分组信息, 同组内的基因表达模式相近, 可能具有相似的功能或参与相同的生物学过程 除了差异基因表达量 FPKM 层次聚类分析, 我们还分别用 H-cluster K-means 和 SOM 等三种方法对差异基因的相对表达水平值 log 2 (ratios) 进行聚类 不同的聚类算法分别将差异基因分为若 干 cluster, 同一 cluster 中的基因在不同的处理条件下具有相似的表达水平变化趋势 图 10 差异基因聚类图 图一 : 整体 FPKM 层次聚类图, 将 log10(fpkm+1) 值进行归一化转换 (scale number) 并进行聚类, 红色表示高表达基因, 蓝色表示低表达基因 颜色从红到蓝, 表示 log10(fpkm+1) 从大到小 ; 图二 :log2(ratios) 折线图, 每个子图中的灰色线条表示一个 cluster 中的基因在不同实验条件下相对表达量, 蓝色线条表示这个 cluster 中的所有基因在不 同实验条件下相对表达量的平均值,x 轴表示实验条件,y 轴表示相对表达量
3.6.4 基因表达维恩图 当只有一组差异分析比较组合时, 绘制基因表达维恩图 该图展示了各 group 表达的基因个数, 以及其重叠关系 图 11 基因表达维恩图 每个大圆圈中的数字之和代表该 group 表达了的基因总个数, 圆圈交叠的部分表示 group 之间共有的表达基因, 以 FPKM>1 为基因表达的标准 3.6.5 差异表达分析 Q&A 问 : 如何判断一个基因是否是差异基因? 如果是差异基因, 如何判断该基因的表达量是上调还是下调? 答 : 差异基因筛选标准参见正文 6.1 中表格 如果基因的 log2foldchange>0, 则认为该差异基因是上调, 反之, 若 log2foldchange<0, 认为该差异基因是下调 问 : 能否用 FPKM/RPKM 进行差异分析? 答 : 在做差异分析时, 诺禾是采用 readcount 数据作为差异分析输入文件, 通过 DESeq 或者 TMM 标准化后, 进行差异分析 FPKM/RPKM 实际上也是对 readcount 进行标准化处理的一种方法, 各种标准化方法优劣势比较见下图 (Dillies, M. A. et al, 2013), 可以看出, 在进行差异分析时,DESeq 和 TMM 的标准化效果最好,FPKM/RPKM 的标准化效果较差, 所以, 不推荐使用 FPKM/RPKM 进行差异分析 问 : 如何判断差异基因在两个样品间的差异大小? 答 : 差异的显著情况可通过 log2foldchange 来判断, log2foldchange 越大, 差异倍数越大 问 : 某基因在两个样本中表达量差别很大, 却不存在与显著差异的基因列表中, 这是为何? 答 : 差异基因的筛选是基于统计学意义的, 不能直观的通过两个数值的大小判断差异基因的是否 : 首先 : 受测序深度的影响, 有些样品的测序深度较深, 可能导致该样品的 readcount 数值较高, 做差异分析的第一步就是要消除测序深度的影响, 对原始数据进行标准化处理 ; 其次 : 在差异分析过程中, 需要对 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, 能否说明一个表达, 一个不表达? 答 : 在有参项目中, 一般默认 FPKM>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 值进行了多重假设检验, 能更好地控制假阳性率 3.7 差异基因 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 被差异基因富集的概率 3.7.1 差异基因 GO 富集列表 表 6 差异基因 GO 富集列表 GO accession Description Term type Over represented p-value Corrected p-value DEG item DEG list GO:0006412 translation biological_process 1.2902e-13 3.9996e-10 76 1051 GO:0005840 ribosome cellular_component 1.9759e-11 3.0626e-08 49 1051 GO:0003735 structural constituent of ribosome molecular_function 4.5917e-11 4.7448e-08 49 1051 GO:0030529 ribonucleoprotein complex cellular_component 1.3658e-10 1.0585e-07 52 1051 结果表格详细内容如下 : GO_accession:Gene Ontology 数据库中唯一的标号信息 Description:Gene Ontology 功能的描述信息 Term_type: 该 GO 的类别 (cellular_component: 细胞组分 ;biological_process: 生物学过程 ;molecular_function: 分子功能 ) Over_represented_pValue: 富集分析统计学显著水平 Corrected_pValue: 矫正后的 P-Value, 一般情况下,Corrected_pValue < 0.05 该功能为富集项 DEG_item: 与该 GO 相关的差异基因的数目 DEG_list:GO 注释的差异基因数目 3.7.2 差异基因 GO 富集柱状图 差异基因 GO 富集柱状图, 直观的反映出在生物过程 (biological process) 细胞组分 (cellular component) 和分子功能 (molecular function) 富集的 GO term 上差异基因的个数分布情况 我们挑选了 富集最显著的 30 个 GO term 在图中展示, 如果不足 30 条, 则全部展示 图 12 GO 富集柱状图 每组两张图 ; 图一 : 纵坐标为富集的 GO term, 横坐标为该 term 中差异基因个数 不同颜色用来区分生物过程 细胞组分和分子功能, 带 * 为显著富集的 GOterm 图二 : 对图一中的 GO, 按生物过程 细胞组分和分子功能三大类别及差异基因上下调分类画的三个子图
3.7.3 差异基因 GO 富集 DAG 图 有向无环图 (Directed Acyclic Graph,DAG) 为差异基因 GO 富集分析结果的图形化展示方式 图中, 分支代表包含关系, 从上至下所定义的功能范围越来越小, 一般选取 GO 富集分析的结 果前 10 位作为有向无环图的主节点, 并通过包含关系, 将相关联的 GO Term 一起展示, 颜色的深浅代表富集程度 我们的项目中分别绘制生物过程 (biological process) 分子功能 (molecular function) 和细胞组分 (cellular component) 的 DAG 图 图 13 GO 富集有向无环图 每个节点代表一个 GO 术语, 方框代表的是富集程度为 TOP10 的 GO, 颜色的深浅代表富集程度, 颜色越深就表示富集程度越高, 每个节点上展示了该 TERM 的名称及富集分析的 p-value 3.7.4 差异基因 GO 富集分析 Q&A 问 :GO 富集分析所使用的软件是什么? 答 :GO 富集分析均采用 R 包, 富集采用的为 GOseq, 有向无环图采用的为 topgo 问 :GO 富集分析一般分析到二级, 还能继续分析吗? 答 :GO 富集分析是针对所有注释的 GO 进行统计检验, 任何等级的都有 问 :GO 富集和 GO 分类有何区别? 答 :GO 分类是将每个基因与其对应的 GO 功能联系起来, 以获取基因的 GO 注释信息, 而 GO 富集分析则是将 GO 功能相似的基因集通过统计学检验算法富集到一起, 从而方便研究具有某 一类 GO 功能的基因 3.8 差异基因 KEGG 富集分析 在生物体内, 不同基因相互协调行使其生物学功能, 通过 Pathway 显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径 KEGG(Kyoto Encyclopedia of Genes and Genomes) 是系统分析基因功能 基因组信息数据库, 它有助于研究者把基因及表达信息作为一个整体网络进行研究 作为 Pathwa 相关的主要公共数据库 (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 富集分析 3.8.1 差异基因 KEGG 富集列表 表 7 差异基因 KEGG 富集列表 Term Database ID Sample number Background number P-Value Corrected P-Value Ribosome KEGG PATHWAY ecd03010 50 78 0.00579498555878 0.486778786938 Bacterial chemotaxis KEGG PATHWAY ecd02030 16 19 0.0219788406991 0.863154067253 Pyruvate metabolism KEGG PATHWAY ecd00620 28 47 0.0594189888655 0.863154067253 Thiamine metabolism KEGG PATHWAY ecd00730 10 12 0.0670066516342 0.863154067253 结果表格详细内容如下 : Term:KEGG 通路的描述信息 ID:KEGG 数据库中通路唯一的编号信息 Sample number: 该通路下差异基因的个数 Background number: 该通路下注释基因的个数 P-value: 富集分析统计学显著水平 Corrected P-value: 矫正后的统计学显著水平,Corrected P-value < 0.05 该功能为富集项
3.8.2 差异基因 KEGG 富集散点图 散点图是 KEGG 富集分析结果的图形化展示方式 在此图中,KEGG 富集程度通过 Rich factor Qvalue 和富集到此通路上的基因个数来衡量 其中 Rich factor 指该 pathway 中富集到的差异基因个数 (Sample number) 与注释基因个数 (Background number) 的比值 Rich factor 越大, 表示富集的程度越大 Qvalue 是做过多重假设检验校正之后的 Pvalue,Qvalue 的取值范围为 [0,1], 越接近于零, 表示富集越显著 我们挑选了富集最显著的 20 条 pathway 条目在该图中进行展示, 若富集的 pathway 条目不足 20 条, 则全部展示 图 14 差异基因 KEGG 富集散点图 纵轴表示 pathway 名称, 横轴表示 Rich factor, 点的大小表示此 pathway 中差异表达基因个数多少, 而点的颜色对应于不同的 Qvalue 范围 3.8.3 富集 KEGG 通路图 为便于查看差异基因在通路图中的分布情况, 我们将差异基因标注到通路图中, 查看方法如下 : 拿到全部分析结果后, 打开 results 结果目录中 DEG_KEGGEnrichment/DEG_KEGGPath 文件夹下相应比较组合的 html 文件, 点击不同的通路名称会弹出相应的通路图 其中, 包含上调基因的 KO 节点标红色, 包含下调基因的 KO 节点标绿色, 包含上下调的标黄色 鼠标悬停于标记的 KO 节点, 弹出差异基因细节框, 标色同上, 括号中数字为 log2(fold change) 以上步骤可脱机实现, 如连接互联网, 点击各个节点, 可以连接到 KEGG 官方数据库中各个 KO 的具体信息页 3.8.4 差异基因 KEGG 富集分析 Q&A 图 15 显著富集的 KEGG pathway 代谢通路图 问 : 为什么编码同一个酶的基因, 会有的上调有的下调? 答 : 这些编号的基因存在着多个条目, 也可能包含了一个家族的多个基因, 它们间的调控机制可能尚不清楚, 反映在图上会有部分上调, 部分下调的现象, 这是比较常见的现象 问 :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
3.9 SNP 和 InDel 分析 3.9.1 SNP 和 InDel 分析 SNP 全称 Single Nucleotide Polymorphisms, 是指在基因组上由单个核苷酸变异形成的遗传标记, 其数量很多, 多态性丰富 从理论上来看每一个 SNP 位点都可以有 4 种不同的变异形式, 但实际上发生的只有两种, 即转换和颠换, 二者之比为 1:2 SNP 在 CG 序列上出现最为频繁, 而且多是 C 转换为 T, 原因是 CG 中的 C 常为甲基化的, 自发地脱氨后即成为胸腺嘧啶 一般而言, SNP 是指变异频率大于 1% 的单核苷酸变异 InDel(insertion-deletion) 是指相对于参考基因组, 样本中发生的小片段的插入缺失, 该插入缺失可能含一个或多个碱基 我们首先通过 samtools 和 picard-tools 等工具对比对结果进行染色体坐标排序 reads 去重复等处理, 然后通过变异检测软件 GATK2 分别进行 SNP Calling 和 InDel Calling, 最后进行过滤, 得到如下表形式的分析结果 其中 InDel 分析结果每列的含义和 SNP 结果是一致的, 这里不再展示 表 8 SNP 分析结果 #CHROM POS REF ALT fjoligo fjh2o Gene ID Gene Name Description CP000948.1 3864720 G A 0,23 0,23 -- -- -- CP000948.1 16764 T C 92,0 144,29 ECDH10B_0018 hokc toxic membrane protein%2c small CP000948.1 236234 T C 201,27 88,71 ECDH10B_0226 thrw trna-thr CP000948.1 504888 G A 0,1 0,1 ECDH10B_0494 ybcc DLP12 prophage%3b predicted exonuclease #CHROM:SNP 位点所在染色体 POS:SNP 位点坐标 REF: 参考序列在该位点的基因型 ALT: 该位点的其它基因型 other coloums: 每个个体该位点的基因型 ( 数字代表支持该位点的 reads 个数 ) Gene ID: SNP 所在基因 ID Gene Name: 基因名称 Description: 基因描述信息 3.9.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: 插入缺失 是指相对于参考基因组, 样本中发生的小片段的插入缺失, 该插入缺失可能含一个或多个碱基 3.10 新转录本预测 用 Rockhopper 软件 (R. McClure, et al, 2013) 将测序结果根据参照参考基因组进行组装, 并与已注释的基因模型进行比较, 发现新的转录本区域 通过 Blastx 与 nr 库作比对 ( 其中,evalue 值设为 le-5), 对新预测的转录本区域进行注释, 将能注释上的转录本区域作为具有编码潜能的新转录本区域 表 9 新转录本预测 Gene_ID Start Stop Strand NR_GI NR_ID Novel00001 101738 101856 + 377887476 EHU51953.1 Novel00002 626612 626716 + Novel00003 739872 739976 + Novel00004 748318 748961 - 转录本编号转录起始位点转录终止位点链方向 NR_GI: 比对到 NR 库的基因的 Genbank ID NR_ID: 比对到 NR 库的基因的序列号
3.11 基因结构分析 原核生物功能上相关的几个基因往往串联排列在一起, 构成操纵子 (Operon) 结构作为基因的表达单元, 受上游共同的调控区和下游转录终止信号的调控 转录时, 几个基因转录在一 条 mrna 链上, 再分别翻译成各自不同的蛋白质 原核生物基因结构和调控模式如下图所示 : 我们通过 Rockhopper 软件, 根据 reads 在基因组上的分布情况, 对操纵子 转录起始位点 (Transcription Start Site, TSS) 和转录终止位点 (Transcription Termination Site, TTS) 进行预测 然 后提取基因转录起始位点上游 700bp 序列, 通过 TDNN(Time-Delay Neural Network) 方法进行启动子预测 3.11.1 操纵子预测 表 10 操纵子预测结果 Start Stop Strand Number of Genes Genes 337 5020 + 3 ECDH10B_0002, ECDH10B_0003, ECDH10B_0004 12163 15298 + 2 ECDH10B_0014, ECDH10B_0015 16751 16960-2 ECDH10B_0017, ECDH10B_0018 17489 19620 + 2 ECDH10B_0020, ECDH10B_0021 Start: 第一个基因的起始坐标 Stop: 最后一个基因的终止坐标 Strand: 链方向 Number of Genes: 基因个数 Genes: 基因名称 3.11.2 TSS 和 TTS 预测 表 11 TSS 和 TTS 预测结果 TSS TTS Strand Gene 105 265 + ECDH10B_0001 328 2800 + ECDH10B_0002 5115 5530 + ECDH10B_0005 6468 5658 - ECDH10B_0006 TSS: 转录起始位点坐标 TTS: 转录终止位点坐标 Strand: 链方向 Genes: 基因名称 3.11.3 启动子预测 表 12 启动子预测结果 Sequence ID strand Start Position End Position Score Sequence ECDH10B_0002.head + 67 112 0.89 CAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACT ECDH10B_0003.head + 358 403 0.87 GAGTTTAACGCCGAGGGTGATGTTGCCGCTTTTATGGCGAATCTGTCACA ECDH10B_0003.head + 478 523 0.94 AATATTGATGAAGATGGCGTCTGCCGCGTGAAGATTGCCGAAGTGGATGG ECDH10B_0005.head + 121 166 0.88 CGTTTTATTGCTGCGACCAACGTGAACGATACCGTGCCACGTTTCCTGCA Sequence ID: 输入的序列 ID strand: 链方向 Start Position: 启动子起始位点坐标 End Position: 启动子终止位点坐标 Score: 精准度 Sequence: 启动子序列
3.12 UTR 分析 3.12.1 UTR 预测和长度分布 我们根据基因转录起始位点 ( 转录终止位点 ) 和翻译起始位点 ( 翻译终止位点 ) 信息, 提取 5'UTR(3'UTR) 序列, 并对其长度分布情况进行统计 针对 5 UTR, 用 RBSfinder 软件对 SD 序列进行预测 ( 其中,rbs region length 我们设为 50bps); 针对 3 UTR, 用 TransTermHP 软件对不依赖 ρ 因子的终止子进行预测 图 16 UTR 长度分布图 横轴表示 UTR 长度区间, 纵轴是不同区间 UTR 密度的统计, 红色虚线代表 UTR 长度的均值 3.12.2 5 UTR SD 序列预测 表 13 5 UTR SD 序列预测 gene_id: 基因编号 Start: 基因起始坐标 Stop: 基因终止坐标 Strand: 链方向 Pattern:SD 序列信息 Gene_ID Start Stop Strand Pattern Position ECDH10B_0008 8238 9191 + AAGAG 8223 ECDH10B_0009 9306 9893 + CGGAA 9293 ECDH10B_0014 12163 14079 + TGGAG 12150 ECDH10B_0015 14168 15298 + AGGGG 14150 Position:SD 序列起始坐标 3.12.3 3 UTR 不依赖 ρ 因子的终止子预测 表 14 3 UTR 不依赖 ρ 因子的终止子预测 Gene_ID Term_start Term_end Strand 5'_tail 5'_stem Loop 3'_stem 3'_tail ECDH10B_0009 9916 9899 - TAGCGAATAAAAAAA TCCCCCC GAGC GGGGGGA TCTCAAAACAAT ECDH10B_0016 16607 16581 - GGAGAATAACAAAAA TGGTCATCTGGA GCT TACAGGTGGCCA TTCGTGGGACAG ECDH10B_0020 18671 18690 + CAGGACGGTTTACCG GGGAGCC ATAAAC GGCTCCC TTTTCATTGTTA ECDH10B_0049 50326 50360 + AGAATTTACGGCTAG CGCCGGATGCGACGC CGGTC GCGTCTTATCCGGCC TTCCTATATCAG gene_id: 基因编号 Term_start: 终止子起始坐标 Term_end: 终止子终止坐标 strand: 链方向 5'_tail: 茎环结构 5' 尾序列 5'_stem: 茎环结构 5' 茎序列 loop: 茎环结构环序列 3'_stem: 茎环结构 3' 茎序列 3'_tail: 茎环结构 3' 尾序列 3.13 反义转录本预测 Cis-natural antisense transcripts(cis-nats) 反义转录本是由源 DNA 链相同区域转录的内生 RNA 分子, 与正义转录本存在部分收敛或分散方向的重复 据目前研究发现, 反义转录本是重要的生物起源机理, 主要通过表观遗传学上的改变, 对基因进行调控 反义转录本分为三种类型 :enclosed( 全部包含 ) convergent (3'-3' 重叠 ) 和 divergent(5'-5' 重叠 ) 对于链特异性建库的 RNA-seq 数据, 可以鉴定其反义转录本在基因组上的位置 种类以及数量等 表 15 反义转录本预测结果
3.14 srna 分析 细菌中, 长度在 50~500nt 的非编码 RNA 通常定义为小 RNA(small RNA, srna) 用 Rockhopper 软件发现新的基因间区转录本, 通过 Blastx 与 nr 库作比对, 对新预测的转录本区域进行注 释, 将注释不上的转录本作为候选的非编码 srna 通过 RNAfold 软件和 IntaRNA 软件对候选的 srna 分别进行二级结构预测和靶基因预测 3.14.1 srna 预测和长度分布 图 17 srna 长度分布图 横轴表示 srna 长度区间, 纵轴是不同区间 srna 密度的统计, 红色虚线代表 srna 长度的均值 3.14.2 srna 二级结构预测 图 18 srna 二级结构 3.14.3 srna 靶基因预测 表 16 srna 靶基因预测结果 srna_id mrna_id energy(kcal/mol) srna_position mrna_position srna00002 ECDH10B_1971-7.93484 125 -- 134 108 -- 118 srna00004 ECDH10B_1971-4.98008 23 -- 44 384 -- 408 srna00008 ECDH10B_1971-8.70759 30 -- 45 268 -- 280 srna00012 ECDH10B_1971-4.5503 42 -- 49 44 -- 51 srna_id:srna 编号 mrna_id: 靶基因编号 energy(kcal/mol): 自由能 srna_position:srna 互补位置 mrna_position: 靶基因互补位置 3.14.4 srna 表达水平分析 表 17 srna 表达水平统计结果 srna_id fjoligo fjh2o srna00042 9.47359944674179 10.3820688547943 srna00030 7.21798053085089 71.1913292900178 srna00031 43.5502780536787 26.0326502627677 srna00004 372.220778262307 590.773208382487 srna_id:srna 编号其他列 :srna 在各样品中的表达量 (fpkm)
四 参考文献 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. Busch, A., A. S. Richter, et al. (2008). IntaRNA: efficient prediction of bacterial srna targets incorporating target site accessibility and seed regions. Bioinformatics.(IntaRNA) Hofacker, I. L. and P. F. Stadler (2006). Memory efficient folding algorithms for circular RNA secondary structures. Bioinformatics.(RNAfold) Kingsford, C. L., K. Ayanbule, et al. (2007).Rapid, accurate, computational discovery of Rho-independent transcription terminators Illuminates their relationship to DNA uptake. Genome biology.(transtermhp) McClure, R., D. Balasubramanian, et al. (2013). Computational analysis of bacterial RNA-Seq data. Nucleic acids research.(rockhopper) Suzek, B. and M. Ermolaeva S. Schreiber and SL Salzberg. (2001). A probabilistic method for identifying start codons in bacterial genomes. Bioinformatics.(RBSfinder) 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) Anders, S.(2010). HTSeq: Analysing high-throughput sequencing data with Python.(HTSeq) 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: a MapReduce 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) Waibel, A. H., Hanazawa, T., Hinton, G. E., Shikano, K., Lang, K. J.. (1989).Phoneme Recognition Using Time-Delay Neural Networks IEEE Transactions on Acoustic, Speech, and Signal Processing Vol. 37 no.3, 328-339.(TDNN)