如何寻找已发表文章的转录组

要寻找已发表文章的转录组，首先需要明确所需研究的基因或生物系统。接着就可以通过在线数据库查找，如GEO、ArrayExpress、SRA等，这些数据库中存储了大量公开发表的原始转录组数据。在这些数据库中，可以使用关键字搜索、分类浏览等多种方式找到与自己研究相关的数据。一旦找到了目标数据，就可以下载并进行数据处理、分析，以获取自己需要的转录组信息。同时，可以通过PubMed等文献检索工具，查找关于该基因或生物系统的相关发表论文，以更深入地了解转录组数据的意义和应用可能。

1）ROI

了解过三代测序数据分析的人，对于CCS 环形一致性序列的概念肯定不会陌生，在iso-seq中，提出了比CCS 更加灵活的一个概念：ROI

ROI , 全称 reads of insert，可以理解为插入片段，首先看下三代测序文库构建阶段的reads示意图：

对于上述的文库片段，测序产生的reads 示意图如下：

由于是一个环状分子， 随着测序反应的进行，会循环测序；如果把插入片段的正负链都测了一次，就做1个full pass;

对于CCS 而言，要求至少有2个full pass , 才能去生成CCS reads; 三代测序的特点就是读长很长，可以达到十几kb, 对于短的插入片段而言，CCS这样定义当然没有问题，但是对于全长转录本

而言，转录本长度很长，比如转录本长度1kb, 读长3kb, 此时在一个零模波导孔（ZMW）中测序的reads 就不可能达到2个full pass , 也就产生不了CCS reads, 为了解决这个问题，提高reads的利用率，提出了ROI 的概念，ROI 指的就是插入片段，上图测序reads 产生的ROI 如下：

ROI 不要求满足2个full pass， 相对CCS 而言，更加适合全长转录本的分析；

2）artifacts

这种序列是由于文库制备阶段，adapter 序列错误的将两条转录本的序列链接构成了一个环状分子，这个和adapter 浓度有关，通常这种reads 产生的比例很少，小于05%， 在后续的分析中，这部分reads 需要去除

在PCR 反应中，由于不完全延伸的产物作为了下次扩增反应的引物，导致出现嵌合体序列，直观上看，就是PCR产物来源于两条或者多条reads；

PCR 产生的嵌合体序列，在PCR 反应体系中，这种序列是不可避免的，大约有3%的比例，在后续的分析过程中，可以借助软件去除这部分reads；

FL , Full-length reads, 全长转录本

从raw data 到 ROI , 在从ROI 去除 artifacts reads 之后，我们就得到了用于后续分析的clean reads;

clean reads 就已经是转录本的序列了，我们首先看一下clean reads 当中，哪些是全长转录本；哪些不是全长转录本，这个 *** 作就是分类，classify

全长转录本的示意图如下：

对于全长转录本而言，其ROI reads 中包含5‘ primer 和 3‘ primer; 而且会出现polyA 为结构；（polyA 针对mRNA和部分lncRNA）

对于不同大小的文库，其全长转录本的比例也不同：

可以看到，文库片段越长，全长转录本的比例越低；

4） consensus transcript isoforms

一致性转录本序列，一个ZMW 产生一个转录本的reads, 肯定会有冗余的reads 出现，这是通过聚类（cluster）的方式，就全长转录本序列进行聚类，可以得到一致性的转录本序列；

数据分析流程： >

一般我们拿到10x空间转录组数据分析的结果最先看的肯定是web_summary网页报告，因为从这个结果里面我们大概就能判断你的数据好不好，不好的问题在哪里，数据到底能不能用等等。这里来详细介绍一下怎么看10x 空间转录组web_summary网页版报告。

10x 空间转录组网页版报告模板如下：

下面来详细介绍一下每块区域每个指标的含义。

Reads总体情况统计区：

Number of Reads： 样本总的测序reads数，双端测序这个是指一端的reads数，实际上算数据量需要用reads 2 读长。

Valid Barcodes： barcode校准后有效的barcode数占总的reads的比例，Space Ranger会先尝试纠正barcode序列中的序列错误，然后再进行统计。

Valid UMIs： 有效的UMI数占总的reads的比例。

Sequencing Saturation： 测序饱和度值，就是在当前测序深度情况下，有多少比例的捕获到的mRNA被测出来了，比如这这里的测序饱和度是74%，说有74%的mRNA基因被检测出来了，如果加大测序深度会有更多的mRNA被检测出来。

Q30 Bases in Barcode： barcode序列的Q30值

Q30 Bases in RNA Read： 捕获的mRNA序列的Q30值

Q30 Bases in UMI： UMI序列的Q30值

Mapping结果统计区：

Reads Mapped to Genome： 比对到基因组上reads的比例

Reads Mapped Confidently to Genome： 唯一比对到基因组上reads的比例，也就是我们常说的mapped uniquely reads，不过这里如果某条reds唯一比对到一个基因的exon区，同时又比对到了一处非exon区，还是算唯一比对到exon区的reads。

Reads Mapped Confidently to Intergenic Regions： 比对到唯一基因间区的reads的比例

Reads Mapped Confidently to Intronic Regions： 比对到唯一内含子区的reads的比例

Reads Mapped Confidently to Exonic Regions： 比对到唯一外显子区的reads的比例

Reads Mapped Confidently toTranscriptome： 比对到唯一基因转录组上reads的比例，这一部分会包括剪切位点的reads。 这一部分的reads就是用来对UMI进行计数统计的 。细心的朋友们可能会发现这一部分的reads比例比Reads Mapped Confidently to Exonic Regions的值要低，这是因为有些基因的exon是有overlap的，处于overlap区域的reads最终是不进入UMI计数的。

Reads Mapped Antisense to Gene： 比对到基因转录组的反义链区域的reads比例，这部分reads是没有意义的。从这里我们也可以发现 10x空间转录组建库和比对是有方向性的 。

Spot信息统计区：

Fraction Reads in Spots Under Tissue： 比对到唯一基因转录组上reads（ Reads Mapped Confidently to Transcriptome ）有多少比例覆盖在组织区域的spot上，这里是93%，那就说明只有7%的reads分布在组织之外的灰色区域的。 10x软件在这里有一个默认的阈值为50% ，认为这个比例值超过50%结果是正常的，低于50%则回到网页最上面区域提示报错信息（认为可能是透化不完全导致背景RNA过高或者是组织区域选的不合适）。从这个50%的阈值上我们也可以判断10x的这个空间转录组技术还是存在一定缺陷的，它允许接近50%的reads散落在组织以外的区域，说明组织透化这一步想让对应区域的mRNA完全都落入对应spot点里面去还是很难的。

Mean Reads per Spot： 每个spot的平均reads数， 10x这里采用的是所以测序reads总是除以组织上检测到的spot数（跟单细胞的统计方法是一样的） ，理论上来说这样统计是不合理的，因为总的reads包括没有比对上的reads、没有mapping到转录本上的reads、组织区域以外的spot上的reads，所以是 不能真实的反应每个spot上实际的reads数的 。

Median Genes per Spot： 每个spot的基因中位数

Total Genes Detected： 检测到的基因总数

Median UMI Counts per Spot： 每个spot的中位UMI数

样本信息区：

Sample ID： 样本id

Chemistry ： 试剂版本

Slide Serial Number ： Slide信号和区域

Reference Path： 参考基因组路径

Transcriptome： 基因组转录组版本

Pipeline Version： spaceranger软件版本

Analysis区域

UMI分布展示： 左边是图像上UMI的分布，右边是tsne降维可视化后的UMI的分布，鼠标放置到图像上会现在对应的位置信息和对应spot上的UMI count数。从这个图我们可以判断UMI主要分布在组织的哪些区域，哪些区域没有捕获到mRNA或捕获的mRNA特别少。

Cluster的分布展示： 左边是cluster在组织图像上的分布，右边是tsne降维可视化后的cluster的分布，鼠标放置到图像上会现在对应的位置信息和对应spot上的cluster值和该cluster占总的spot的比例。这个上cluster分群在组织上的层次关系特别明显。

这一部分主要展示亚群的top基因的信息，因为不管是单细胞还是空间转录组基本上后面都会自己另外重新分析的，所以这部分和上面的cluster分布信息意义不大。

Sequencing Saturation（测序饱和度）

对reads进行随机抽样，观察不同spot平均reads的情况下测序饱和度的分析，一直到实际的测序深度测序饱和度的值，理论上当所有转化的mRNA转录本均已测序后，饱和度接近10（100％），虚线表示测序到合理的饱和点位置，也就是说就是测序深度再高也不可能饱和度达到100%。

Median Genes per Spot（sopt点的中位基因）

也是对reads进行随机抽样，观察不同spot平均reads的情况下spot的中位基因的值，曲线最高点的斜率能反应增加测序深度能得到最大的spot的中位基因数。

总结

对于web_summary的结果我们大概重点可以从下面几个方面来看数据效果

1、总的spot数，这个其实由组织的大小而定，没有具体好坏的说法；

2、每个spot的中位基因数，中位基因数太少说明捕获效果不好，有可能透化步骤条件不够优化，当然也有可能是试剂或芯片的问题；

3、测序饱和度，每个点的UMI中位数，sopt平均reads数，饱和度、sopt平均reads数和中位UMI数都太低说明测序深度不够，需要加大测序量。

4、基因组的比对率，比对率太低有可能是样品污染；

5、组织spot上reads的比例，比对太低有可能透化时间不够导致背景RNA过高，需要优化透化条件，也有可能组织区域选的不好，这个可以通过LoupeBrowser手动选择组织区域。

所谓时序分析 (time series analysis) 在 data science 中是非常重要的一个方向。对大多数商业行为而言如果能够通过已有不同时间数据来进行预测就有可能大大提高自己的胜率。通常时间序列数据会包括趋势部分和不规则部分， 我们需要做的就是剔除不规则部分然后找到趋势所在，再进行预测。在预测过程中通常可以采用移动平均法、局部加权回归法、指数平滑法和自回归整合移动平均等方法。

生物学的时间相关数据本身预测属性和商业数据相比要弱很多。一种是单一条件的纯时间序列，主要看不同基因的表达模式，根据相似的表达谱将基因归为多个类有助于找到功能相似的基因。另一种情况是含有对照和处理的时间序列，需要再考察不同条件的差异基因。

关于时间序列转录组数据分析的工具，近三年来有两篇偏综述和测评类的文章（一个人写的）。

在这两篇文章中还是提到了一些工具，但其中有一些用到matlab（这软件贵啊），有一些年久失修或者不维护或者和最新R版本不兼容，筛筛捡捡能用的且文章里认为还不错的也就剩下三四个。

来自于 DESeq 的方法，下文中提到的 ImpluseDE2 和 MaSigPro 都使用了这种模型。

来自于 maSigPro 方法，所谓多项式回归区别常见的线性回归，会把一次特征转换成高次特征的线性组合多项式，比使用直线拟合更加准确。但是到底用几次方需要具体分析，次数过高会出现过拟合。在能够解释自变量和因变量关系的前提下，次数应该是越低越好，这也算是奥卡姆剃刀原则吧。

所谓自回归是统计上一种处理时间序列的方法，用 至 来预测本期 的表现并假设它们为线性关系。简单说就是用自己来预测自己，因为是从回归分析中线性回归发展而来只是用x预测x，所以叫自回归。

同样是来自于 DESeq 的方法，下文中提到的 ImpluseDE2 和 MaSigPro 也都使用了这种方法。

似然比检验 (likelihood ratio test,LRT) 用于比较两个模型的拟合优度进而确定哪个模型与样本数据拟合的更好。其中一个是具有一定数量项的完整模型，另一个是删掉完整模型中一部分项的简化模型。LRT 检验中，自由度等于在简化模型中减少的模型参数数目，LR近似符合卡方分布。一个相对复杂的模型与一个简单模型比较，如果可以显著地适合一个特定数据集，那么这个复杂模型的附加参数就能够用在以后的数据分析中。

为了测试多个时间点的任何差异，可以使用包含时间因子的设计和时间因子在简化公式中被删除的另一个设计。对于包括对照和实验组的时间序列，可以使用包含条件因子，时间因子和两者相互作用的公式。在这种情况下，使用具有不包含相互作用项的简化模型的似然比检验将测试该条件是否在参考水平时间点（time 0）之后的任何时间点可以诱导基因表达的变化。

EBseq-HMM 采用的方法，来自于 BEseq。

这个软件最早发表在2007年，相对老一些好在目前仍然在维护，其主要目的是给时序数据进行基于 模糊聚类算法 的聚类。我们常见的聚类算法可以分为严格聚类(hard clustering)和模糊聚类(Fuzzy clustering )(也叫做宽松聚类 soft clustering)。严格聚类会将一个基因只聚到一类中，kmeans 就属于严格聚类。而模糊聚类允许同一数据属于多个不同的类，其聚类结果是一个数据对聚类中心的隶属度，0到1之间。对于分类很开的数据使用严格聚类是没问题的。但对于时序表达量数据来说，不同的类常常会有重叠，所以可以尝试宽松聚类方法。算法需要首先设定一些参数，若初始化参数不合适，可能影响聚类结果的正确性。

在使用 Mfuzz 时首先应该进行数据标准化处理 ，可以使用类似于 FPKM 或者 TPM 的表达结果也可以使用 DESeq2 矫正后的结果进行比较分析，另外不支持值为0的数据，所以需要加上 pseudocount 。除此之外，Mfuzz 接受的数据格式为 ExpressionSet，需要对矩阵进行转换。

这个包只能进行聚类，是找不了有处理对照组的差异基因的。需要注意。

运行masigpro 主要有四步：

有两点内容需要注意：对于无对照的单一时序数据处理方法；以及处理转录数据时的特殊参数。因为这个包不会对数据进行标准化，所以应该提前做好，使用 DESeq2 即可。

另外，在实际分析的时候可能会出现 glmfit: algorithm did not converge 的警告。这是由于进行 logistic 回归时，依照极大似然估计原则进行迭代求解回归系数，glm 函数默认最大迭代次数是 25，当数据不太好时 25 次迭代可能还不收敛，一方面可以增大迭代次数。但当增大迭代次数仍然不收敛就需要对数据进行异常值检验等进一步处理。通常把一些表达量极低或者极高的基因删除掉，这个问题就可以解决。

ImpulseDE2 是最近才出来的一个R包，在前面提到的综述评测文章中认为这个包找时序数据中的差异基因效果最好，它可以用来解决两类问题。

这个包中，有一个 plotHeatmap 函数，可以借助 ComplexHeatmap 对数据整体进行热图的绘制同时提取不同类的基因，也可以使用 plotGenes 看某一个基因的表达情况。

在展示的热图中会出现四部分，包括 transient and transition trajectorie，其中每一种 tarjectorie 又包括 up 和 down 两类。所谓的 transient 可以理解为时序数据在中间某一个时间点存在up 或者 down peak，即在某一个时间点存在表达的最大或者最小值；而所谓的 transient 可以理解为一个持续的变化，比如持续的升高或者持续的降低。

EBSeq-HMM 是基于 EBSeq 二次开发的工具，主要用于分析时序数据。在计算的时候首先基于负二项分布对参数进行估计，然后利用自回归隐马模型将基因的表达进行分类。比较神奇的是，最终给到的结果会标示为 Up-Up-Down-Down-Down 之类的若干 path，然后你可以选出你感兴趣的 path 进行后续分析。

因为目前做的数据是没有对照的单一时间序列数据，所以还不能体会哪一个找出的差异基因更准确些。但是如果只是想把所有的基因根据不同的时间点分为若干表达 pattern，似乎结合 Mfuzz 和 ImpulseDE2 就可以了。

当然，涉及到聚类，尤其是非监督聚类的时候通常主观因素还是较强，如果能对关键基因或者数据有一个大致的估计预判 *** 作起来会相对轻松些，如果没有，可能就需要结合不同类的生物学意义等角度来找合适的聚类数目了。

>