有参转录组分析

Posted 2023-04-24

参考技术A • 总体上来说HISAT利用了BWA和Bowtie的算法，同时解决了mRNA中不存在内含子难以比对的问题，比对上代主流RNA-seq比对工具能快50倍，同时需求更少的内存<8G（这就意味着你可以在PC上跑数据），20个样本，每个样本一亿reads的估计，用一台电脑一天之内能够跑完。使用者可以提供精确的基因注释来提高在已知基因区域的准确性，但这是可选项。

• Transcript assembly and quantification with StringTie

• RNA-seq的分析依赖于精准的对于基因isoform的重建以及对于基因相对丰度的预测。继承于Cufflinks，StringTie相对于之前开发的工具更为准确，需求内存和耗时也更少。

• 使用者一样可以使用注释文件来帮助StringTie运行，对于低丰度的数据比较有帮助。此时StringTie依旧会对非注释区域进行转录本的组装，所以注释文件在这里是可选选项。

首先添加环境变量，如果使用conda的话需要 source activate

hisat2支持对gtf文件构建索引，并且可以向其中添加外显子，SNP等信息，进一步完善索引

随后利用 ballgown 或者 featurecount 得到表达矩阵进行下游分析

对于可变剪切一般使用 rMATs 进行统计，rMATS是一款对RNA-Seq数据进行差异可变剪切分析的软件。其通过rMATS统计模型对不同样本（有生物学重复的）进行可变剪切事件的表达定量，然后以likelihood-ratio test计算P value来表示两组样品在IncLevel（Inclusion Level）水平上的差异(从公式上来看，IncLevel跟PSI的定义也是类似的)，lncLevel并利用Benjamini Hochberg算法对p value进行校正得FDR值。rMATS可识别的可变剪切事件有5种，分别是skipped exon (SE)外显子跳跃，alternative 5′ splice site (A5SS)第一个外显子可变剪切，alternative 3′ splice site (A3SS)最后一个外显子可变剪切，mutually exclusive exons (MXE)外显子选择性跳跃和 retained intron (RI)内含子滞留

随后对结果进行可视化，

简单的转录组差异基因表达分析 -- DESeq2

参考技术A

经典的转录组差异分析通常会使用到三个工具 limma/voom , edgeR 和 DESeq2 。今天我们就通过一个小规模的转录组测序数据来演示 DESeq2 的简单流程。

对于 DESeq2 的分析流程而言，我们需要输入的数据包括：

​ 下面就以 mobData 中的数据为例简单介绍 DESeq2 的分析流程

​ 由于 mobData 中的行名没有提供基因的ID，我们也不是为了探究生物学问题，就以 mobData 的行数作为其ID

DESeqDataSet 是 DESeq2 流程中储存read counts和中间统计分析数据的对象，之后的分析都建立在该对象之上进行。

​ 在进行差异分析之前，需要对样本数据的表达矩阵进行预处理，包括：

​ 通过PCA结果来看各组样本分组情况还是不错的，但hclust的聚类结果反映的分组就略微有点混杂了，可能要聚类计算的距离函数选用不当有关。

使用 DESeq() 函数进行差异分析时，该函数干了以下三件事：

​ counts() 可以提取 DESeq object 中的表达矩阵，而 results() 可以提取差异分析的结果，其中包括了：

样本间的均值, log2 fold changes, standard errors, test statistics, p-values and adjusted p-values.

​ 使用 results() 函数时需要指明进行比较的样本，这里用 contrast=c("group_list","MM","WW") 提取 MM 组和 WW 组进行差异分析的结果。如果想要比较 WM 组和 WW 组，只要改变 contrast=c("group_list","WM","WW") 即可。

​ 检查结果中是否包含 NA 值

​ 这里 padj 中有1142个 NA 值是因为使用 results() 提取差异分析结果时，大于 alpha 值（这里是0.1）的矫正后p-value都会被当做是 NA 。因此，我们将这些 padj 值都设为 1

排序后以 log2FoldChange 绝对值大于1， padj 小于0.05为条件筛选显著的差异表达基因

​ 至此，便筛选出了 217个 在 MM 组和 WW 组之间的显著差异表达基因。至于后续的可视化分析则是因课题而异了，等以后有空了再补坑吧！

​ 有同学可能注意到，虽然我们的样本有多个组，但在差异分析时进行的还是pairwise的分析，为什么我们不可以三个组一起分析呢？学过ANOVA的同学应该都知道，ANOVA就是可以应对这种多组差异分析的情况。但要注意的是，ANOVA只可以告诉我们对于某个基因在这三组中是否存在差异，想要找出是哪一组有其他组别有差异还是需要进行pairwise t-test之类的分析。所以，在这里我们两组两组地进行分析正是出于这个考虑，并且有更方便我们解释差异分析的结果，说明在A组基因的表达量相对于B组的是上调还是下调。另外，本文的差异分析还是处于单因子水平（只有一个变量），至于多因子的差异分析以后研究透了再和大家进行分享。

完。