搜索
转录组结果解读北京诺禾致源科技股份有限公司转录调控研究部OUTLINE简介实验部分生物信息分析概述1转录组是指特定组织或细胞在某个时间或某个状态下转录出来的所有RNA的总和，主要包括mRNA和非编码RNA。转录组研究是研究基因功能和结构的基础，对生物体的发育和疾病的发生具有重要作用。RNA-seq技术流程主要包含两个部分，建库测序和数据分析。实验部分（RNA检测、建库、测序)）2•琼脂糖凝胶电泳：分析样品RNA完整性及是否存在杂质污染。•NanoPhotometerspectrophotometer：检测RNA纯度（OD260/280及OD260/230比值）。•Agilent2100bioanalyzer：精确检测RNA完整性。链特异性文库优势：相同数据量下可获取更多有效信息；能获得更精准的基因定量、定位与注释信息5➢1、一般动物样品会有三条带：28S、18S、5S，如果提取过程经过过柱处理或者利用CTAB+LiCl方法提取，5S可能较暗或者没有。➢昆虫或者软体动物等样品只有1条比较明显的带，例如：牡蛎、果蝇、螨虫、蝗虫、蚊、蚕等➢2、植物样品有三条带：25S、18S、5S，有些特殊物种或部位可能本身含条带比较多，如果条带清晰，也可初步判定合格➢3.原核生物中主要有5S、16S、23SrRNA叶片小鼠蚊动物植物原核RIN5RIN7RIN8RIN9RIN4RIN6RIN10RIN2RIN1RIN值范围示意图问与答文献要求OD260/OD230≥1.8，OD260/OD230如果小于2.0，则表明样品中被碳水化合物、盐类或有机溶剂污染；OD260/OD230合格的标准是多少呢？答：OD260/OD230≥2.0，且OD260/OD280≥2.0这说明RNA提取结果是相当好的，一般在1.8-2.1之间就说明RNA结果十分好，但是nanodrop的灵敏度没有2100好，因此我们主要根据2100检测结果来判定RNA是否合格，一般只要RIN值和RNA总量达到我们的判定标准的话，我们就会判为合格。生物信息分析（有参转录组）3有参转录组生物信息分析（医口转录组）3医口转录组生物信息分析31、数据质控1、1测序数据说明见结果文件：QC生物信息分析31、2测序数据过滤对原始数据进行过滤：去除带接头(adapter)的reads；•去除含N(N表示无法确定碱基信息)的reads；•去除低质量reads(Qphred=20的碱基数占整个read长度的50％以上的reads)。每个样本的测序数据过滤情况，如图所示。生物信息分析31、3测序错误率分布测序错误率分布存在以下两个特征：•测序错误率随着测序序列长度的增加而升高。原因：测序过程中化学试剂的消耗导致的，为Illumina高通量测序平台所具有的特征。•前6个碱基具有较高的测序错误率，反转录所需的随机引物和RNA模版的不完全结合单个碱基位置的测序错误率应该低于1%，最高不超过6%生物信息分析31、4GC含量分布GC含量在物种间存在一定特异性，反转录使用的6bp随机引物，前几位碱基在核苷酸组成上有一定偏好性，产生正常波动，随后则趋于稳定。普通建库方法，由于序列的随机性打断和双链互补等原则，理论上测序读段在每个位置的GC及AT含量应分别相等，且在整个测序过程基本稳定不变，呈水平线。而对于链特异性建库而言，由于只保留了单链信息，可能会出现AT分离或GC分离现象。生物信息分析31.5数据质量汇总•Q20：Phred数值大于20的碱基占总碱基的百分比•Q30：Phred数值大于30的碱基占总碱基的百分比生物信息分析（有参/医口）32、1比对率统计（HISAT2比对软件)）HISAT2的算法主要分为三个部分：•将测序序列整段比对到基因组单外显子•将测序序列分段比对到基因组的两个外显子上•将测序序列分段比对到基因组三个以上（含三个）外显子mapping高于70%(TotalReads)：参考基因组组装的较为完善所测物种与参考基因组一致相关实验不存在污染。mapping时用的是read全长，还是头尾有处理？测序得到的read1和read2的各个碱基全都是样本的序列2、比对分析生物信息分析（有参/医口）32、2数据质量汇总表格见结果文件：QC/5.Stat/align_pct.xls生物信息分析（有参/医口）32.3比对区域分布1、注释较为完善的物种比对到外显子区域的比例很高。2、比对到内含子：可能来源于前体mRNA、可变剪接事件滞留的内含子。3、比对到基因间区，可能来源于ncRNA、少许DNA片段污染，基因注释还不够完善。所有样本的测序reads在基因组区域分布情况如图所示，见结果文件：QC/4.Region。生物信息分析（有参/医口）32.4比对可视化显示reads在各染色体上的分布及在基因组中注释的外显子、内含子、基因间区等功能区域的分布，如下图所示，使用说明文档IGVQuickStart染色体的长度越长，map到的reads数目越多吗？——不一定，前提是基因在染色体上分布的密度是相同的，且这些基因的转录水平相当生物信息分析（有参/医口）33新基因预测3.1新转录本组装（StringTie软件）相对于cufflinks等软件，StringTie有以下优势：(1)拼接出更完整的转录本；(2)拼接出更准确的转录本；(3)更好的估计转录本的表达水平;(4)拼接速度更快。新基因预测的意义？非模式物种，其基因注释信息通常不是很完善，新基因预测可挖掘该物种新的基因或转录本。见结果文件2.Assemble生物信息分析（有参/医口）3StringTie预测的新转录本结构注释GTF格式内容，见结果文件：Assemble/*_novel.gtf。•V1seqname：染色体编号•V2source：注释的来源，这里的StringTie是指该转录本是由StringTie软件组装所得•V3feature：注释信息的结构类型，如gene、transcript、exon等•V9attributes：包含众多属性的列表，主要为基因编号、转录本编号等信息生物信息分析（有参/医口）33.2新转录本注释我们会对新转录本进行Pfam、SUPERFAMILY、GO、KEGG等数据库注释，其中Pfam结果如下表所示，其余数据库注释信息见结果文件：Assemble/*_novel_gene.xls。生物信息分析（有参/医口）34定量分析（featureCounts）见结果文件3.Quant生物信息分析34.1基因表达分布RNA-seq的基因表达值一般不用readcount来表示，而是用FPKM，FPKM先后对测序深度和基因长度进行了校正。有参：FPKM大于1，认为表达无参：FPKM大于0.3，认为表达FPKMintervals代表基因表达水平分为5个区间：低表达（0-1）、中表达（1-3）、中高表达（3-15）、高表达（15-60）、超高表达（60）生物信息分析34.2样本间相关性生物学重复主要有两个用途：1、实验可以重复；2、获得更可靠的差异结果（样本选择合理）相关性系数越高，其表达模式越为接近，样本相关性热图如下图所示，见结果文件：Quant/correlation.svg。生物信息分析（有参/医口）34.3主成分分析对所有样本的基因表达值（FPKM）进行PCA分析，如图所示。理想条件下，PCA图中，组间样本应该分散，组内样本应该聚在一起，见结果文件：Quant/pca.svg。问与答为什么要做生物学重复1、发表文章的需要2、不同个体、不同的处理、甚至某一瞬间样品的基因表达也是存在差异的，RNA-Seq数据可能会表现出比预期的更高的假阳性变异性，通常会通过生物学重复来屏蔽掉生物学内部变异大的不稳定的差异，得到真正的处理间的差异。3、在分析方面，增加生物学重复主要是为了减少生物学重复之间的噪音对分析结果的影响，简言之，如果组内差异大于组间差异，这种情况下得到的差异基因假阳性会高一些。生物信息分析（有参/医口）35差异分析差异分析主要分为三个步骤。•首先对原始的readcount进行标准化（normalization），主要是对测序深度的校正。•然后统计学模型进行假设检验概率（pvalue）的计算•最后进行多重假设检验校正，得到FDR值（错误发现率)。大部分差异分析软件（DESeq，DESeq2和edgeR）用原始的readcount作为输入文件，这些软件自身对会readcount做一些校正（主要是测序深度），而FPKM是校正后的表达值，所以用FPKM做差异分析相当于做了两次校正，是不合理的见结果文件4.Differential生物信息分析（有参/医口）35.1差异基因列表每个比较组合的差异显著性分析如下表所示，见结果文件：Differential/1.deglist。如何判断差异基因在两个样品间的差异大小？答：padj越小，差异越显著。也可通过|log2Foldchange|来判断差异的大小情况，|log2Foldchange|越大，差异倍数越大生物信息分析（有参/医口）35.2差异基因统计每个比较组合的差异基因（包括上调和下调）数目统计以及筛选差异的标准如下表所示，见结果文件：Differential/1.deglist/diff_stat.xls。FDR=padj=correctedpval生物信息分析3结果文件：Differential/1.deglist/{比较组合}/_volcano.png。问与答某基因在两个样本中表达量差别很大，却不存在与显著差异的基因列表中，这是为何？——差异基因的筛选是基于统计学意义的，不能直观的通过两个数值的大小判断差异基因的是否；在有重复的项目中，如果重复较差，组内差异情况会屏蔽掉部分组间的差异。所以会导致差异基因过少，再次：在计算完pvalue以后，需要对pvalue进行多重假设检验校正，得到padj，来减少假阳性。使得部分通过pvalue阀值的基因，无法通过padj的阀值。差异基因筛选条件最大能设的阈值是多少？答：最大可设阀值没有定论，一般等级比较高的文章卡的阀值都比较的严格。在一些PLOSONE等文章里面，可能卡的值比较的松，有的文章会在无重复中，只卡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，能否说明一个表达，一个不表达？——这里的readcount值仍为校正后的readcount，由于软件显示的原因，实际的readcount并不一定是0，在无参项目中，一般默认rpkm0.3时，基因表达；在有参项目中，一般默认rpkm1时，基因表达。不推荐看readcount的值看判断表达与否。能否提取部分基因来做差异分析？答：不能。差异分析是基于整体来做的。差异分析软件的作者推荐用全部readcount进行差异分析，若使用部分基因做分析，会毁坏掉数据整体的特点，如测序深度、reads分布特征。每次选定差异基因筛选阈值后，软件都会重新计算pval、padj值，同一阈值下两次计算pval/padj值可能会有不同，尤其是处于阈值附近的pval/padj值，所以差异基因筛选结果会出现微弱的不同，对数据整体的影响不大，可正常使用。为什么两组比较的readcounts值不一样，即AvsB中A的readcount值和AvsC中的A的readcount值不一样？差异基因是根据两两组对比得来的，而不