thymeleaf 排除list中的某个数据_分析展示你的RNA-seq数据，从这里开始

公司用illumina测序对建好库的RNA样品进行测序后，会得到一堆后缀为的Rawdata。然后在经过公司或者实验室人员将Rawdata进行比对后，得到了表达矩阵的数据。那么怎么对这几万个基因进行分析呢？有什么策略可以看到你想看到的东西呢？

一. 处理数据之前，我们先要了解数据类型，

通常比对完的数据有两种。

一种是经过均一化处理过后的FPKM（或者也可以是RPKM、TPM等，不同的均一化方式）。另一种数据是未经过均一化处理的count的数据。这类样本不可以进行直接比较，而是要经过标准化之后才能比较。

那么上述两类标准化的方法有什么不同呢？

首先我们要知道RNA-seq的数据为什么要标准化，RNA-seq要解决的一个关键问题就在于定量，像qPCR一样，这样不同样本才能比较，而这些标准化的方法主要想解决两个问题：

我们一个个介绍：

FPKM的计算公式如图

其中C是比对到该基因的外显子上的片断数，N是所有map至基因组的reads的碱基数，L就是该基因外显子碱基全长。

简单的来说FPKM均一化的方式考虑了总基因数和基因长度，所以不同样本的测序深度差异它是可以解决的。人们主要是用FPKM和RPKM来当统计量。不过TPM较具优势，因为该数值可以直接进行样本间的比较。它的不同在于先除以长度再除以这个的总和，这种方式算出来的不同样品的TPM总和是相同的，（也意味着TPM数值能体现出比对上某个基因的reads的比例。）

需要注意的是cuffdiff计算FPKM时会根据你的样本矩阵来调整，即是最终得到FPKM值并不是组内样品单独运算得到的FPKM值的平均值。同个矩阵里跑的若是同批处理的样本可以直接比较。

FPKM 的计算过程可以参照这里： https:// zhuanlan.zhihu.com/p/37 196518 https:// qiubio.com/archives/309 1

• 如何标准化FPKM呢？如果拿到的FPKM出现的严重的批次效应（batch effect），你需要将其去掉。

去掉batch effect的方法一个是可以用housekeeping gene 来这类不变的量来标准化， 一个可以用combat包或者limma包来进行Quantile标准化具体见

具体标准化方法可以参考这两篇
RNA-seq中的那些统计学问题（二）FPKM/RPKM之外的那些标准化方法 https://www. omicsclass.com/article/ 1113

• 标准化count的方法有许多，如R包deseq2、limma、edgeR等，而这些包的输入也只能是count，而不能是做过均一化的FPKM等。（这是由后序分析需要用到的统计学方法决定的。）

deseq2标准化的优势在于，它不仅对测序深度进行了标准化，而且有文库补偿（弥补library composition）的功能。因为不同的样本含有不同的特异高表达的基因，这些基因会对总基因数，从而对其它基因的表达量造成影响。比如，若X基因只有某一个样本A高表达，FPKM的计算时在不表达的B样本中为了文库一致，这个低表达X基因的B样本就会拉高同样本中的其它基因的表达量，而A样本中除了X基因的其它基因就会被拉低，可能导致些许阳假性的差异。

简单的来说，deseq2标准化就是找标准化因子，有如下几个步骤：

a. 取自然对数（以e为底，用对数的一个原因就是对数后均值不太容易受到异常值的影响）

b. 对数完每行的均值（相同基因，这种方法也叫几何平均数）

c. 去掉-Inf。即所有样本中有表达量为0的基因，这样做有两个好处，一个是把所有表达量为0的都去掉了，一个就是可以把某个样本高表达，有些样品不表达的基因给去掉。那么剩下的就是所有样品都有表达的管家的基因了。

d. 对数矩阵减均值。其本质上是以某个基因的平均值为参考，对每个样本中的基因X进行均一化。这步也是计算过程中出现负值的原因，因为count都是整数，原本取对数后不会出现有负数的情况。而对数后相减后就有可能出现负数，比如。这步可以理解为是考虑文库补偿的问题，是FPKM做不到的。

y=ln(x)对数函数

e. 计算每个样本对数的中位数，不用平均数是为了排除一些极端表达基因的影响。

f. 将中位数再变为对数前的数，这就是每列的标准化因子。

g. 将raw count除以每列的标准化因子，得到标准化后的矩阵。

看起来有些复杂，但其实你只要输入count，这个软件一步就能完成。你只要记得，deseq2只是一个差异分析的软件，就是类似于做方差分析的软件一样，只不过其通过log变换和中位数挑选来排除异常值的影响。

deseq2矫正的原理可以看原北卡罗来纳大学教堂山分校的Josh Starme的StatQuest系列视频教程，里边的统计学原理值得学习，也有人将这个系列的视频整理成了学习笔记 StatQuest学习笔记25——差异表达分析

二. 介绍完两种基本数据类型后，我们以我们用TCGA上下载的肝癌和胆管癌RNA-seq数据来举例说明一下分析过程。

我们在得到数据后，对样本的整体情况要有一个大致的判断，这样才能保证数据分析前没有问题。

1. 各样本表达的情况。用箱线图看一下，不同样品之间的表达量的均值要相对一致。若不一致，后序要用经过标准化至均值一致才能比较。

我们用TCGA上下载的肝癌和胆管癌的count数据来举个例子，可以看到count在标准化前中位数值不一。

Raw count矩阵分布

而经过deseq2七步矫正标准化过后，中位数值明显统一了。

标准化过后的count分布

2. 计算样品之间的皮尔森相关系数（PCC）。对样本的相关性作一个判断，检查同类样本的重复性，不同样本的相似程度。

PCC of FPKM（这里因为样本太多就挑选了一些样本）

3. 层次聚类（hierarchical clustering）。目的同上，生物学上相近的样品应该聚在一起，不同的相距较远。

三. 上述几个标准都符合后，我们就可以开始对数据进行分析了，首先是看你的分析目的。

RNA-seq可以做的大都是相关性研究，通过比较找到一些差异，从基因表达上给你的课题指明一定的方向，一般来说，单独做RNA-seq，有如下几个常见的目的。

1. 如果你的样本是实验组与对照组的关系，那么寻找差异基因是关键，这可以通过RNA变化来推测蛋白的差异。找差异基因的软件有很多，上面说的处理count的软件都可以。（这里需要注意的是，unchange的基因应该是去掉所有表达量都为0的那些的基因后的list。）

2. 如果你的样本不止两组，那么可以用k-means聚类来所有样品中各类基因的表达分布，研究你想要的类别。

3. 如果样本具有递进关系，比如时间上，或者病情上，那么也可以分析表达模式（Expression patterns）来找你想要研究的cluster。

我选取的样本没有明显递进关系，找的paper的图

4. 如果你想看两个样本list之间的交集，可以用venn图来找到相应的基因

一些venn的方法做法 http:// bioinformatics.psb.ugent.be /webtools/Venn/
或者是sangerbox里的venn也很好用

5. 如果你想看你的样品的亲缘关系，你可以使用上面的层次聚类，也可以使用降维分析PCA的方法来看。甚至有些时候还可以通过PCA的图观察到同类细胞的变化趋势，比如细胞可塑性。

6. 当然你也可以直接找你关注的基因通路的list。

上述3个方法都可以得到你的想研究的gene list，一个个找你想要的基因当然可以，不过差异基因有上千个时候就很麻烦了，通常我们会用基因富集通路的方法去找两组样本是否在某些通路上得到富集，从而在那个方向上设计实验。

四. 那么通路富集的方法有哪些呢？

1. 网站Metascape，方便快捷，可选数据库。Metascape

2. 网站g:Profiler，它虽然数据库分类的挺多，可视化也做的好看，但在计算p-value方面，可能算的偏大。a web server for functional enrichment analysis and conversions of gene lists

3. Y叔的R包ClusterProfiler ，是业界大神Y叔写的一个R包，除了做通路富集，还集成了许多方便使用的可视化命令。

利用clusterProfiler可视化数据

详情可以参考 https://www. jianshu.com/p/e133ab316 9fa https://www. jianshu.com/p/feaefcbdf 986 >

4. GESA软件 GESA是电脑软件，按照格式调整表达矩阵，扔到软件里，就可以看基因list是否在某些通路上富集了。 可能某些通路并不是1-2方法中富集靠前，但又想知道该通路有无富集就可以用这个办法。

GESA具体链接 https://www. jianshu.com/p/c0ea67728 907

五. List 的可视化，可视化是把数据呈现给大家的最后一步。

1. 对于gene list 可以直接用pheatmap用作图，其中注意一些细节。

2. 点图可以用Plot，也可以用ggplot2来画一些好看一些的图，比如火山图，条形图，拆线图等，或者用sangerbox里的图也很好用。

今天就先写到这里啦，终于码完了，希望我的拙见可以帮助到你。有机会再分享一些简单的代码和R语言学习小技巧。