转录组ref流程工作手册一、Reference流程生物学原理1.1实验流程图一：转录组实验流程当我们得到样品时，必须对其测序，才能得到分析所需的数据。测序根本过 程：提取样品总RNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA假设为 原核生物，那么用试剂盒去除rRNA后进入下一步)。参加fragmentation buffer 将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物(random hexamers) 合成第一条cDNA链，然后参加缓冲液、dNTPs、RNase H和DNA polymerase I 合成第二条cDNA链，在经过QiaQuick PCR试剂盒纯化并加EB缓冲液洗脱之 后做末端修复并连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后 进行PCR扩增，使用建好的测序文库进行测序。得到RNA的序列后，乂可以找到它的参考序列物种本身的基因、基因组)3发现新转录本现有数据库中对转录本的注释可能还不全面,通过高通量测序我们能检测到新的转录本 (Mortazavi, 2021 #103。我们首先从潜在gene model中挑选出长度大于I5()bp且平均覆盖度 大于2的gene model,再从中找出位于基因间区域一个基因3,端下游200bp到下一个基因 5,端上游2(X)bp之间的区域)的潜在gene model作为候选的新转录本。4基因结构以及Reads在基因组上分布的精确图形该分析主要是以图形方式概括给出Reads在基因组各个位置的分布情况，以及该位置基因的 分布情况。我们画出Reads在最长的25条染色体上的分布图，该图为SVG矢量图，如果你的 浏览器不支持SVG,请安装SVGView插件。5基因差异表达分析5.1基因表达量基因表这最的计算使用 RPKM 法(Reads Per Kb per Million reads) Mortazavi, 2021#103,其计算公式为：RPKS1106CAZ/10-设RPKM(A)为基因A的表达量，那么C为唯一比对到基因A的reads数，N为唯一比 对到基因组的总reads数，L为基因A编码区的碱基数。RPKM法能消除基因长度和测序量 差异对计算基因表达的影响，计算得到的基因表达量可直接用于比拟不同样品间的基因表达 差异。如果一个基因存在多个转录本，那么用该基因的最长转录本计算其测序覆盖度和表达 量。5.2差异分析差异表达分析找出在不同样本间存在差异表达的基因，并对差异表达基因做GO功能 分析和KEGG Pathway分析。参照Audic S.等人发表在Genome Research上的基于测序的差异基因检测方法Audic, 1997 #8该文献己被引用超过五百次)，我们开发了严格的算法筛选两样本间的差异表达 基因。假设观测到基因A对应的reads数为x,在一个大文库中，每个基因的表达量只占所有 基因表达量的一小局部，在这种情况下，p(x)的分布服从泊松分布：骸为基因A的真实转剥X!,样本一中唯一比对到基因组的总reads数为N1,样本二中唯一比对到基因组的总reads数 为N2,样本一中唯一比对到基因A的总reads数为x,样本二中唯一比对到基因A的总reads 数为y,那么基因A在两样本中表达量相等的概率可由以下公式计算：2二的对r-0(如果套力心0.5)x-0(x+力！0，或 2x(l-p(i|x)1-0夕。| 对=(=)，一 、N x!，!(l + )g，D，M然后，我们对差异检验的p value作多重假设检验校正，通过控制FDR(False Discovery Rate) 来决定p value的域值。假设挑选了 R个差异表达基因，其中S个是真正有差异表达的基因， 另外V个是其实没有差异表达的基因，为假阳性结果。希望错误比例Q=V/R平均而言不 能超过某个可以容忍的值，比方1 %,那么在统计时预先设定FDR不能超过0.01 (Benjamini, Yekutieli, 2001)。在得到差异检验的FDR值同时，我们根据基因的表达量(RPKM值)计 算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，那么说明表达差异 越显著。在我们的分析中，差异表达基因定义为FDRW0.001且倍数差异在2倍以上的基因。 得到差异表达基因之后，我们对差异表达基因做GO功能分析和KEGG Pathway分析。GO功能分析一方面给出差异表达基因的GO功能分类注释；另一方面给出差异表达基因的 GO功能显著性富集分析。GO功能分类注释给出具有某个GO功能的基因列表及基因数目统计。GO功能显著性富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO功能条 目，从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库(： geneontology.org/)的各个 term 映射，计算每个 term 的基因数目，然后应用超儿何检验，找出与整个基因组背景相比，在差异表达基因中显著富集 的GO条F1,其计算公式为m-1 MN-MP_1 V z其中，N为所有基因中具有GO注释的基因数目；n为N中差异表达基因的数目：M为所有 基因中注释为某特定GO term的基因数目；m为注释为某特定GO term的差异表达基因数 Fl o计算得到的p value通过Bonferroni校正之后，以corrected p value0.05为阈值，满足此 条件的GO term定义为在差异表达基因中显著富集的GO term。通过GO功能显著性富集分 析能确定差异表达基因行使的主要生物学功能。我们的GO功能分析同时整合了表达模式聚类分析，研究人员能方便地看到具有某功能的 所有差异基因的表达模式。例，immune response为在差异表达基因中最显著富集的一个GO term表2)。图十三显示了参与immune response的差异基因的表达模式。表2在差异表达基因中显著富集的GO-termTerms from the Process Ontology with p-value as good or better than 0.05Gene Ontology termCluster frequencyGenome frequency of useCorrected P-valueExpression Profileimmune response view aenes82 out of 807 genes,10.2%663 out of 13525 genes, 4.9%2.74e-07View Resultimmune system process view aenes100 out of 807 genes, 12 4%921 out of 13525 genes, 6.8%3.77e-06View Resultresponse to virus view qenes21 out of 807 genes,26%105 out of 13525 genes, 0.8%0.00138View Resultreaulation of apoptosis view aenes63 out of 807 genes, 7.8%583 out of 13525 genes. 4.3%0 00508View Resultrelation of programmed cell death view Qenes63 out of 807 genes,7.8%592 out of 13525 genes, 4.4%0.00840View Result&requlation of cell death view aenes63 out of 807 genes,7 8%593 out of 13525 genes, 4 4%0.00888View Resultrequlation of cell death view qenes63 out of 807 genes,7 8%593 out of 13525 genes, 44%0.00888View Result3800mllp Vs 0m 5500m Vs 0m 3800mDown vs 0m图十三参与i mmune response的差异基因表达模式聚类图KEGG Pathway 分析在生物体内，不同基因相互协调行使其生物学功能，基于Pathway的分析有助于更进一步了 解基因的生物学功能。KEGG是有关Pathway的主要公共数据Kanehisa, 2021 #96, Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背 景相比，在差异表达基因中显著性富集的Pathway。该分析的计算公式同GO功能显著性富 集分析，在这里N为所有基因中具有Pathway注释的基因数目；n%N中差异表达基因的 数目；M为所有基因中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的 差异表达基因数IR。FDR0.05的Pathway定义为在差异表达基因中显著富集的Pathway。 通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。 结果如表3所示。表3 pathway显著性富集分析列表#PathwayDEGs with pathway annotation (2085)All genes with pathway annotation (8986)PxalueQvaluePathway ID1M etaboHc pathways307 (14.72%)1081 (12.03%)1.354119e-050.002911356koOHOO2Proteasome23 (1.1%)48 (0.53%)0.00014825700.015937627ko030503B cell receptor signaling pathwav29(1.39%)70(0.78%)0.00050853410.036444944ko046624Apoptosis34 (1.63%)89 (0.996)0.0010184710.045737882ko042105Hematopoietic cell lineage31 (1 49%)80(0 89%)0.0012719050 045737882ko046406Primary immunodcficiencY16(0.77%)33 (0.37%)0.0012764060.045737882ko053407Glvcosdpho$phaiidviinositoKGPD-anchoi13(0.62%)25 (028%)0.0016188250.049721054ko00563biosynthests8N-Glyc 邳 biosynthesis18(0 86#)40 (0.45%)0.0019011400 051093137ko005109Huntingtons dkease60(2.88%)1