使用TPM/FPKM/RPKM进行差异分析真的可以消除系统误差吗？

[ixxmu]

https://mp.weixin.qq.com/s/L5ikz6qmCTSoO4ut4ahQvg

[ixxmu]

使用TPM/FPKM/RPKM进行差异分析真的可以消除系统误差吗？ by 生信菜鸟团

这周曾老师给我分享了一篇文章，题为《Don’t trust TPM/FPKM/RPKM too much. They don’t promise to cancel the systematic error of RNA-Seq data》，文章主要讨论了RNA-Seq数据中的系统误差不能被TPM/FPKM/RPKM消除的问题。TPM/FPKM/RPKM是一种“全局归一化”，在理论上可以消除线性系统误差。但是，作者认为这种假设对组学数据分析来说太天真了，并且强调组学数据的非线性系统误差没有生物信息学方法可以消除，并以GSE159751为例进行了相关分析。

---

背景知识

首先我们来谈谈什么是线性和非线性系统误差？为什么TPM/FPKM/RPKM在理论上可以消除线性系统误差？

线性误差

组学数据中的线性系统误差指的是由于技术或实验条件的限制，导致测量值与真实值之间存在一定的偏差，这种偏差在数学上可以用线性模型来表示。具体而言，对于一组测量数据 y，它们与真实值 x 的关系可以用以下的线性模型表示：y = Ax + b + e

其中，A 是系数矩阵，b 是偏置向量，e 是误差向量。在这个模型中，A 和 b 是已知的常数，而 e 表示测量误差，也就是线性系统误差。

消除线性系统误差的方法主要有两种：校准和标准化。

校准的基本思路是通过添加一些已知的标准样品，建立标准曲线，然后利用标准曲线对未知样品的测量结果进行修正。具体而言，我们可以将标准样品的真实值和测量值用线性回归模型拟合，得到一个系数和截距，然后利用这个系数和截距对未知样品的测量值进行校准。

标准化的基本思路是将测量数据进行线性变换，使得它们的均值为零，方差为一。这个过程通常称为Z-score标准化：

Z = (x - μ) / σ

其中，x 是原始测量数据，μ 是所有样品的均值，σ 是所有样品的标准差。将这个标准化的结果作为最终的测量结果，可以有效消除线性系统误差e。

下面是一些在转录组数据分析中可能遇到的线性误差的例子：

1. 序列长度：不同长度的转录本会影响到测序数据的比对和计数。这个误差可以通过对测序数据进行长度标准化来消除。
2. 测序深度：测序深度的不同会影响到表达量的计算。

这些误差可以通过对表达量进行归一化处理来消除，例如，TPM、FPKM、RPKM等方法。

关于这些方法我们此前也有介绍：

• RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同

• Counts FPKM RPKM TPM CPM 的转化

也欢迎大家关注我们语雀平台团队整理的“生物统计从理论到实践”，里面涉及到许多生物信息的统计学知识讲解

RPKM FPKM 和 TPM：https://www.yuque.com/biotrainee/biostat/chapter3-4

RPKM FPKM 和 TPM主要目的是去除测序数据的技术偏差：测序深度和基因长度

3. 不完全的RNA提取：对于不完全的RNA提取，可能会影响到所测定的RNA的数量，从而影响到表达量的计算。这个误差可以通过对样品进行RNA质量控制、RNA浓度检测、DNase I处理等步骤来消除。
4. PCR扩增引入的偏差：PCR扩增过程中，不同的转录本可能会被扩增的效率不同，从而引入偏差。这个误差可以通过进行PCR扩增的优化来消除。
5. RNA降解：RNA的降解会导致RNA浓度的下降和RNA质量的降低，从而影响到表达量的计算。这个误差可以通过对样品进行RNA质量控制和标准化处理来消除。

这些误差通常会引起表达量测量的偏差，但是这些误差通常是线性的，因此可以采用不同的标准化方法和归一化方法来消除或减小这些误差，从而提高数据处理的准确性和可靠性。

非线性误差

组学数据中的非线性系统误差指的是由于技术或实验条件的限制，导致测量值与真实值之间存在一定的偏差，而这种偏差无法用线性模型来完全描述。

与线性系统误差不同，非线性系统误差不仅仅是误差向量 e 的线性组合，而可能涉及到多个变量之间的复杂关系，如幂律、指数函数、对数函数等。

因此，对于组学数据中的非线性系统误差，我们需要采用更加复杂的模型来描述和修正。常见的方法包括使用高阶多项式、神经网络等非线性模型来对数据进行拟合和校正。

非线性系统误差是指在组学数据分析中可能存在的不符合线性模型的偏差。这种偏差可能由于实验设计、样本处理、RNA测量、数据分析等各个方面的因素引起，因此可能影响分析的准确性和可靠性。非线性系统误差不能通过简单的计算和全局归一化等基本方法来消除，需要更为复杂和有效的技术来处理。在组学数据分析之前，必须对可能出现的非线性系统误差有所了解并对其进行评估，以保证数据的质量和准确性。

下面是一些在转录组数据分析中可能遇到的非线性误差的例子：

1. 批次效应：如果实验中存在多个批次，可能会导致不同批次之间的表达量存在显著差异，这是因为批次效应会影响RNA提取、样品制备和测序等步骤，从而影响表达量测量的准确性。
2. 技术因素：转录组测序技术中存在多个步骤，如RNA提取、RNA库制备和测序等，每个步骤都可能会引入一些非线性误差，例如，PCR扩增过程中可能会引入偏差等。
3. 剪接变异：基因的剪接变异会导致不同转录本长度的不同，因此在标准化方法中只考虑基因长度而不考虑剪接变异可能会引入非线性误差。
4. RNA质量：RNA质量是影响表达量测量的重要因素，如果RNA质量不好，则会导致表达量测量偏低或偏高，从而引入非线性误差。

所以作者认为组学数据的非线性系统误差没有生物信息学方法可以消除，在此基础上，任何不同来源的数据集的合并都是不可接受的。

---

小编有话说：

此外，根据全篇文章的解读以及小编后续复现的分析也可发现，作者的文章暗含了支持使用raw counts进行差异分析的立场，认为为了理论上消除系统误差而使用TPM/FPKM/RPKM是不可取的。所以稍微总结一下作者的两个立场，以便读者带着问题阅读：

• TPM/FPKM/RPKM并虽然理论上可以消除线性系统误差，但对差异分析而言结果却不如raw counts不理想，这在一定程度反映了这些均一化方法的局限性和处理带来的潜在偏差
• 在实验设计和操作环节需要尤其注意非线性误差的潜在影响，因为没有生物信息学方法可以处理非线性误差，任何不同来源的数据集的合并都是不可接受的

---

以下引用均为文章作者所述翻译内容

histogram

您可能认为必须使用TPM/FPKM/RPKM，而不是raw counts，以避免样本之间的系统性误差。但每百万reads的标准化是一种“全局标准化”，理论上只能抵消线性系统误差。然而，对于组学数据分析来说，这样的假设太天真了。让我们看一个示例，GSE159751，它由两组26个样本组成。根据GSM记录，GSE159751_RAW.tar提供的矩阵的为FPKM值。请注意，以直方图表示的FPKM的分布模式各不相同。似乎有两种类型的形状：单峰型和双峰型，与疾病状态无关。这有力地表明存在系统性错误。（Figure 1）

首先来到GEO看看这个数据集https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159751

FPKM (Fragments Per Kilobase of transcript per Million mapped reads) 是一种常用的转录组表达量的标准化方法。FPKM处理的转录组数据的分布直方图通常呈现出类似于正态分布的形式，这是由于转录本的表达量在大多数生物学实验中都是连续变量，同时在样本中存在许多不同的转录本。

虽然FPKM转录组数据的分布可以近似为正态分布，但是具体的形态和趋势仍然取决于具体的实验设计和生物体。在一些情况下，转录本表达量可能呈现出双峰分布或者其他类型的非正态分布（这里FPKM的分布模式有两种形状：单峰和双峰状，无论疾病状态如何）。

我们先尝试使用GEOquery包在线获取：

library(GEOquery) 
library(AnnoProbe)

eset <- getGEO("GSE159751", destdir=".", AnnotGPL = F, getGPL = F)
class(eset)
view(eset)

exp <- exprs(eset[[1]])  # 空值

发现为空后转为手动下载GSE159751_RAW文件，并对各个样本GSM进行合并获得FPKM的表达矩阵merge_GSM_file.txt

##### 多样本直方图-------------
mydata <- read.table("merge_GSM_file.txt", header=TRUE, sep=" ")
d <- read.table("GSE159751_RAW/GSM4838989_N1.txt",header = T,sep = '\t')
coln <- d[,1]  # gene列
mydata$gene <- coln
mydata$X <- NULL
# reshape2::merlt把长数据转换为宽数据
library(reshape2)
library(tidyverse)  # 数据分析和可视化的大型library集合，包含有dplyr、ggplot2等常用库
library(hrbrthemes)  # 提供了多个精美的主题，可以美化图表
library(viridis)  # 提供了多个平滑的配色方案，可以用来设置颜色
library(forcats)  # 针对离散变量(factors)的再次封装，提供了更多的独特编号和排序功能
long_data <- mydata %>%
  melt(id.vars = c("gene"),  # 指定保留的列
       variable.name = "sample",  # 通过多列列名新生成列的名称
       value.name = "FPKM")  # 通过多列值生成新列的名称
# gene列这里没啥用，丢了算了
long_data$gene <- NULL
# 根据sample分组信息 plot
myplot <- long_data %>%
  # 使用fct_reorder()函数调整sample的顺序
  # mutate(sample = fct_reorder(sample, FPKM)) %>%
  ggplot( aes(x=FPKM, color=sample, fill=sample)) +
  geom_histogram(alpha=0.6) +
  # 使用“scale_fill_viridis”和“scale_color_viridis”函数添加填充和颜色，以使每个文本标签具有不同的颜色
  scale_x_continuous(trans = "log10",  # log10 区域narrow down
                     labels = scales::label_number(accuracy = 1)) +  # 标签不使用科学计数法
  scale_fill_viridis(discrete=TRUE) +
  scale_color_viridis(discrete=TRUE) +
  # 使用“theme_ipsum”函数将主题设置为“ipsum”主题，这是一种较新的主题，可使图形看起来更美观
  theme_ipsum() +
  # 使用“theme”函数对其他主题元素进行微调，如删除图例，调整面板间距等
  theme(
    legend.position="none",
    panel.spacing = unit(0.1, "lines"),
    strip.text.x = element_text(size = 8)
  ) +
  # 使用“xlab”和“ylab”函数分别添加x轴和y轴标签
  xlab("the histograms of the distribution of FPKM of each sample.") +
  ylab("counts") +
  # facet_wrap()是ggplot2包中的一个函数，它用于将数据集中的不同组分别绘制在不同的小图中，以避免在同一轴上绘制过多分组而造成图像混乱的问题。这个函数可以根据数据集中的一个或多个变量，生成多个小图。每个小图都可以使用相同的绘图参数。通过facet_wrap()函数生成的多个小图组成了一张大图，方便对数据进行比较和分析
  facet_wrap(~sample)
myplot
# N8.1列其实为C8 下载的数据集错把C8写成了N8

可以发现，各样本FPKM值的分布确实如作者所说“FPKM的分布模式有两种形状：单峰和双峰状，无论疾病状态如何”

---

PCA

随后PCA结果表明，样本表达谱很可能是由于非线性偏差而聚集在一起的，而不是由于生物学背景。请看主成分分析的结果，疾病状态（蓝色和黄色）与PC1（横轴）无关（图2，上图）。所以我将样本标记为PC1得分高和低（粉色和绿色）（图2，下部）。So I marked samples as PC1 score high and low (pink and green).

rm(list=ls())
data <- read.table("merge_GSM_file.txt",header = T,sep = ' ')
genecol <- read.table("GSE159751_RAW/GSM4838989_N1.txt",header = T,sep = '\t')
genecol <- genecol[,1]
data$X <- NULL
data$gene <- genecol
save(data,file = "merge_GSMdata_withGene.RData")

rm(list = ls())
load("1.merge_GSMdata_withGene.RData")  # data 表达矩阵

options(stringsAsFactors = F) 
getOption('timeout')  # 60
options(timeout=10000)
library(FactoMineR)
library(factoextra)

group=rep(c("Non-IBD","CD"),each=13)
group_list=factor(group,levels = c("Non-IBD","CD"))
table(group_list)#检查一下组别数量
# Non-IBD      CD 
# 13      13

# PCA
data[1:4,1:4]
# N1       N2          N3          N4
# 1 4.936070  4.70440  6.96842018  7.35134879
# 2 0.285694  0.00000  0.06803618  0.03205427
# 3 9.950503 12.70891 20.55543822 17.82955841
# 4 1.425634  1.57114  1.24927397  1.94193101
# gene列转换为行名
rownames(data) <- data$gene
data$gene <- NULL
# 添加group列
data <- as.data.frame(t(data))
data$group <- group_list
data[1:4,1:4]
# ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457
# N1        4.936070      0.28569398        9.950503        1.425634
# N2        4.704400      0.00000000       12.708906        1.571140
# N3        6.968420      0.06803618       20.555438        1.249274
# N4        7.351349      0.03205427       17.829558        1.941931
table(data$group)
# Non-IBD      CD 
# 13      13 
# 绘图
dat_pca <- PCA(data[,-ncol(data)], graph = FALSE)  # 去除非数值列 pca结果
PCAplot<- fviz_pca_ind(dat_pca,
                     geom.ind = "point", # 只显示点，不显示文字
                     col.ind = data$group, # 用不同颜色表示分组
                     palette = c("#00AFBB", "#E7B800"),
                     addEllipses = T, # 是否圈起来，少于4个样圈不起来
                     legend.title = "disease state groups") + theme_bw()
PCAplot

# So I marked samples as PC1 score high and low (pink and green).
pca_group <- rep("PC1 high socres",26)
pca_group[which(dat_pca$ind$coord[,1]<0)] <- "PC1 low scores"
table(pca_group)
# PC1 high socres  PC1 low scores 
# 13              13 
PCAplot_byPCA1<- fviz_pca_ind(dat_pca,
                       geom.ind = "point", # 只显示点，不显示文字
                       col.ind = pca_group, # 用不同颜色表示分组
                       palette = c("#00AFBB", "#E7B800"),
                       addEllipses = T, # 是否圈起来，少于4个样圈不起来
                       legend.title = "PC1 scores groups") + theme_bw()
PCAplot_byPCA1

save(pca_group, group_list,file = "2.PCA.RData")

可以发现样本确实没有按照“disease state”分组聚集

---

hierarchical clustering

然后我对数据执行了分层聚类。当然，样本不是按疾病状态进行聚类，而是按预期的PC1评分类别进行聚类。更有趣的是，PC1分数分类与直方图的形状有关。（图3）

从这里开始，我们的结果开始和作者出现差别，看着作者的图3，私以为，虽然样本不是按疾病状态进行聚类，但也不能完全说按预期的PC1评分类别进行聚类吧，聚类结果也不是很完美

一开始我犯懒打算用SD标准差取差异基因简单看看聚类，但根据PC1的聚类结果都没作者那样好，这里还是贴出来，后面我还是做了一遍差异分析用差异基因热图聚类（虽然根据PC1得分分组的聚类结果还是没作者那样好）

# the R Graph Gallery：
# The heatmap() function is natively provided in R. It produces high quality matrix and offers statistical tools to normalize input data, run clustering algorithm and visualize the result with dendrograms. It is one of the very rare case where I prefer base R to ggplot2.
# ggplot2 also allows to build heatmaps thanks to geom_tile(). However, I personally prefer the heatmap() function above since only it offers option for normalization, clustering and Dendrogram.
# 所以这里我们不用ggplot2使用base
rm(list=ls())
load("1.merge_GSMdata_withGene.RData")
load("2.PCA.RData")

head(data)
rownames(data) <- data$gene
data$gene <- NULL

# > dim(data_mat)
# [1] 42072    26

data$sd <-  apply(data,1,sd)
chosen_data <- data[data$sd>2, ]  # 选择sd标准差大于2的
dim(chosen_data)
# [1] 7134   27
chosen_data$sd <- NULL

############### plot sd>2-------------------------
data_mat <- as.matrix(chosen_data)
# col_DiseaseState <- RColorBrewer::brewer.pal(2, "Set1")[group_list]
col_PC1 <- RColorBrewer::brewer.pal(2,"Set2")[factor(pca_group)]
# colSide <- c(col_DiseaseState,col_PC1)
# heatmap(data_mat,ColSideColors = col_DiseaseState)
# 太大了显示不出来
heatmap(data_mat,
        ColSideColors = col_PC1,
        xlab = "Disease State Groups",labRow = "",  # 去除横坐标
        hclustfun = hclust,  # default
        margins = c(5,8))  # 避免和后面添加的legend重叠  numeric vector of length 2 containing the margins (see par(mar = *)) for column and row names, respectively.

# Create a vector of colors and their corresponding labels
PC1_labels <- c("low","high")
PC1_color <- c("#FC8D62", "#66C2A5")  # 根据col_PC1
# Create a legend for the color vector
legend("topright", legend = PC1_labels, fill = PC1_color, 
       bty = "n", cex = 0.5, title = "PC1 groups")

###### 选择sd标准差大于50的-------------------
chosen_data1 <- data[data$sd>50, ]  # 411
chosen_data1$sd <- NULL
data_mat1 <- as.matrix(chosen_data1)
# col_DiseaseState <- RColorBrewer::brewer.pal(2, "Set1")[group_list]
col_PC1 <- RColorBrewer::brewer.pal(2,"Set2")[factor(pca_group)]
# colSide <- c(col_DiseaseState,col_PC1)
# heatmap(data_mat1,ColSideColors = col_DiseaseState,
#         hclustfun = hclust)  # 默认
heatmap(data_mat1,
        ColSideColors = col_PC1,
        xlab = "Disease State Groups",labRow = "",  # 去除横坐标
        hclustfun = hclust,  # default
        margins = c(5,8))  # 避免和后面添加的legend重叠  numeric vector of length 2 containing the margins (see par(mar = *)) for column and row names, respectively.

# Create a legend for the color vector
legend("topright", legend = PC1_labels, fill = PC1_color, 
       bty = "n", cex = 0.5, title = "PC1 groups")

##### 选择sd标准差前100----------
library(tidyverse)
chosen_data2 <- tail(arrange(data,sd), 100)
chosen_data2$sd <- NULL
data_mat2 <- as.matrix(chosen_data2)
# col_DiseaseState <- RColorBrewer::brewer.pal(2, "Set1")[group_list]
col_PC1 <- RColorBrewer::brewer.pal(2,"Set2")[factor(pca_group)]
# colSide <- c(col_DiseaseState,col_PC1)
heatmap(data_mat2,
        ColSideColors = col_PC1,
        xlab = "Disease State Groups",labRow = "",  # 去除横坐标
        hclustfun = hclust,  # default
        margins = c(5,8))  # 避免和后面添加的legend重叠  numeric vector of length 2 containing the margins (see par(mar = *)) for column and row names, respectively.

# Create a legend for the color vector
legend("topright", legend = PC1_labels, fill = PC1_color, 
       bty = "n", cex = 0.5, title = "PC1 groups")

SD>2SD>50SD top100

差异分析：

rm(list=ls())
load("1.merge_GSMdata_withGene.RData")
load("2.PCA.RData")
######将FPKM值转换为TPM值--------------
expMatrix <- data[,c(1:26)]
rownames(expMatrix) <- data$gene
fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
tpms <- apply(expMatrix,2,fpkmToTpm) #每列的基因
tpms[1:3,]
colSums(tpms) #都是每一百万（深度一样）

#####差异分析----------------
group_list
#表达矩阵数据校正
exprSet <- tpms
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

#判断数据是否需要转换
exprSet <- log2(exprSet+1)
#差异分析：
dat <- exprSet
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
# 获取差异表达基因
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg)
# logFC AveExpr      t   P.Value adj.P.Val     B
# ENSG00000185499  3.422   2.750  7.371 9.587e-08  0.002646 7.704
# ENSG00000119535  2.734   2.272  7.190 1.478e-07  0.002646 7.317
# ENSG00000128203  1.349   1.800  7.081 1.922e-07  0.002646 7.082
# ENSG00000101470 -2.496   3.279 -6.969 2.516e-07  0.002646 6.841
# ENSG00000211670  4.152   3.780  6.873 3.179e-07  0.002675 6.631
# ENSG00000010932 -3.515   3.148 -6.659 5.371e-07  0.003722 6.158
save(dat,deg,file = '3.2.deg.Rdata')

# load("3.2.deg.Rdata")
## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
if(T){
  logFC_t=2
  deg$g=ifelse(deg$P.Value>0.05,'stable',
               ifelse( deg$logFC > logFC_t,'UP',
                       ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
  )
  table(deg$g)
  # DOWN stable     UP 
  # 91  41858    123
  head(deg)
}

#####差异基因样本聚类热图------------

sigGenes <- rownames(deg[deg$g!="stable",])  # 214
# 用的是上一步通过差异分析得到的显著基因和tpm的结果
dat_deg <- dat[match(sigGenes,rownames(dat)),]
dat_deg[1:4,1:4]
# N1     N2    N3     N4
# ENSG00000185499 0.5743 0.7200 1.368 0.7253
# ENSG00000119535 0.7383 1.4715 1.012 0.3869
# ENSG00000101470 4.3641 4.5919 5.006 5.0930
# ENSG00000211670 0.4918 0.4594 1.199 1.2891
group_list
# 数据框形式分组
# disease
group_disease <- data.frame(group=group_list)
rownames(group_disease)=colnames(dat_deg)
library(pheatmap)
pheatmap(dat_deg,
         scale = "row",  # 标准化 非归一化
         show_colnames =T,show_rownames = F, 
         cluster_cols = T, 
         annotation_col=group_disease,
         main = "DEG_disease_cluster_rowScale")

# PC1
pca_group
group_pc1 <- data.frame(group=pca_group)
rownames(group_pc1)=colnames(dat_deg)
library(pheatmap)
pheatmap(dat_deg,
         scale = "row",
         show_colnames =T,show_rownames = F, 
         cluster_cols = T, 
         annotation_col=group_pc1,
         main = "DEG_PC1_cluster_rowScale")

我们的聚类结果可以看出，虽然“disease state”分组的聚类效果也不是很好，但比PC1要好

作者认为“disease state”分组的聚类效果不好是因为存在非线性系统差异而生物学背景被覆盖掉

单独看根据“disease state”的PCA和热图聚类的效果确实可以认为存在非线性系统差异，但我们没有按作者预期那样得到根据PC1分组聚类的结果

您可能希望应用更复杂的规范化技术来强制形状相似。所以我使用了分位数quantile标准化并再次聚类。然而，结果没有改变：PC1高样本和低样本再次聚集在一起。（图4）从大规模的微阵列数据分析中已经知道，使分布形状相似并不能纠正非线性系统误差，这对RNA-Seq也是有效的。你必须明白，没有生物信息学方法可以消除组学数据的非线性系统误差。更糟糕的是，这种复杂的归一化技术向研究人员隐瞒了数据质量方面的问题。（请了解RMA或Z分数标准化的情况。）

我们在上面的差异分析的过程在对TPM矩阵进行了log2，并在绘制热图进行聚类的过程中对每行进行了scale标准化

需要注意的是我们一开始进行差异分析时就画了boxplot查看数据分布，并去除组间差异，这应该与作者所说的“强制分布形状相似”一致

其实我看作者又画的这个图，PC1分组聚类效果感觉还是没作者说的那么好

接着作者为了进一步确定非线性系统误差的来源，从头获取fastq文件用其他软件做了一遍：

我给你看另一个实验。我下载了FASTQ文件，并用其他算法对它们进行了重新量化；fastp、HISAT2和StringTie。你可以预期，不同的量化会产生不同形状的直方图。但是，PC1得分高和低的样本仍然被再次聚类。（图5）因此，FASTQ文件中一定存在内在差异，并且它主导了所有后续的分析工作。

数据集上传作者使用的软件：于是，作者将非线性系统误差的来源锁定到了获取fastq文件等实验操作、测序工作上

DE NOVO

获取数据

下载 fastq 数据（推荐）从 EBI 数据库直接下载 fastq 格式数据，直接越过 SRA，而且不用再做格式转换https://sra-explorer.info/

质控

hisat2比对

stringtie定量

这个软件我之前没用过，这里写的详细点，转录组专辑之前在上游分析推文中使用的都是FeatureCount

参考文章：

• 学徒作业-hisat2+stringtie+ballgown流程   https://mp.weixin.qq.com/s?src=11&timestamp=1682431784&ver=4490&signature=J9bCWUQrWood3U9A81VgVa2WRUfBl76OdAUhsTS1GfIQ0ciq1dU7Tz3tbWVulXstT8fysenGBi2gzqMr3XGz3SiFI0kQmIug1tK8kV0xQKWPFsfIQrayk*qsjCoy1WVv&new=1

• RNA-seq : Hisat2+Stringtie+DESeq2   https://mp.weixin.qq.com/s?__biz=MzkyMTI1MTYxNA==&mid=2247484051&idx=1&sn=aba8f000f86dc9bfb6e330d0591ab302&scene=21&token=1159231349&lang=zh_CN#wechat_redirect

Stringtie和FeatureCount都是用于转录组定量的工具，但它们使用不同的定量策略。

在FeatureCount的情况下，它直接将读数映射到预先存在的注释（如参考基因组或基因注释），并将读数分配给相应的特征，然后可以用于量化。

另一方面，Stringtie使用从头开始的方法从读取中组装转录本并生成新的注释文件。然后使用该组装的转录组来估计来自同一组读数的表达水平。进行合并步骤以合并多个样本并生成一致转录组，该转录组可用于定量所有样本中的读数，从而提高准确性和再现性。

这两种方法各有优缺点。FeatureCount的分析速度更快、更直接，但它只依赖于预先存在的注释，可能会错过以前没有注释的新的或替代的转录本。Stringtie速度较慢，但可以捕获新的或替代的转录物，并提供了对转录物丰度的改进评估，尤其是对于低表达或重叠的转录物。

当然stringtie其实也可以直接定量，如果不需要发现新的转录本和基因的话，可直接基于 bam 文件走定量步骤，我们这里还是使用从头开的方法学习一下这个软件。从头开始的方法和软件具体命令参数可参考上面给出的文章。

组装

合并

定量

我们这里下载的是py3的版本：prepDE.py (python3) : http://ccb.jhu.edu/software/stringtie/dl/prepDE.py3

这里可以发现，我们只用了样本名.gtf文件，而没有用gene_abundances.tsv文件

打开看看就可以发现，样本名.gtf文件是包括了转录水平和基因水平的，而gene_abundances.tsv只有基因水平。prepDE.py 脚本可以直接从 gtf 文件中提取 count 值，而不需要使用 gene_abundances.tsv 文件。这样可以省去合并 gene_abundances.tsv 表格的步骤，简化了数据处理流程。

---

结束后在当前目录生成 gene_count_matrix.csv和 transcript_count_matrix.csv文件，我们要用的是基因水平上文件 gene_count_matrix.csv

到了这一步我们其实就得到了熟悉的raw counts矩阵了，打开看会发现gene_id列有许多“MSTRG”标志符，是因为我们这里使用的是从头开始的方法

MSTRG是StringTie软件生成的新基因的命名前缀，表示未在基因组注释中发现的新基因。在gene_count_matrix.csv文件中出现大量的MSTRG表示这些新基因在各个样本中都有被表达并被计算在了count矩阵中

MSTRG 是 StringTie 生成的转录本（transcript）的名称，当 StringTie 无法将转录本与已知的基因注释匹配时，会生成一个新的 MSTRG 转录本。在生成 count 矩阵时，StringTie 会将 MSTRG 转录本与已知的基因注释进行合并，并将其与已知的基因一起计算。因此，gene_count_matrix.csv 中出现的 MSTRG 转录本实际上是 StringTie 生成的新转录本和已知基因的合并结果

获取已知基因看看有多少：

48946个，和GEO给的数量42072个差不多

提取 FPKM/TPM 或 coverage 结果

我们在这里只提取FPKM 

需要借助一个perl脚本：https://github.com/griffithlab/rnaseq_tutorial/blob/master/scripts/stringtie_expression_matrix.pl

发现这个脚本--result-dirs必须要求gtf和基因丰度tsv文件在以样本命名的文件夹下

手动调整：

大致这样改文件名

这中间我有个地方看错了，来来回回改名字有点乱就不贴了。反正最后就是改成这个样：

获取逗号分隔样本文件名：运行脚本提取FPKM:

怪事出现了，脚本说它找不到transcript_id，但文件里明明是有的

查看脚本所在仓库issue：https://github.com/griffithlab/rnaseq_tutorial/issues/46https://github.com/griffithlab/rnaseq_tutorial/issues/48

发现存在相同问题：

由于我本人不会perl所以我打开脚本大致看了看没看懂哪里有问题

所以接下来我们使用raw counts矩阵进行分析

• 转录组差异分析FPKM与count处理差别大吗

---

raw counts 差异分析

PCA

# 用counts进行差异分析PC1试一试
rm(list=ls())
data <- read.table("known_gene_count_matrix.csv",sep = ',',header = T)  # 表达矩阵
rownames(data) <- data$gene_id
data$gene_id <- NULL

options(stringsAsFactors = F) 
getOption('timeout')  # 60
options(timeout=10000)
library(FactoMineR)
library(factoextra)

group=rep(c("Non-IBD","CD"),each=13)
group_list=factor(group,levels = c("Non-IBD","CD"))
table(group_list)#检查一下组别数量
# Non-IBD      CD 
# 13      13

# PCA
data[1:4,1:4]
# SRR12858287 SRR12858288 SRR12858289 SRR12858290
# ENSG00000288531                    0          11          11           0
# ENSG00000230368|FAM41C             8           0           0           0
# ENSG00000283040                    0           0           0           0
# ENSG00000234711|TUBB8P11           1           2           0           0

# 添加group列
data <- as.data.frame(t(data))
data$group <- group_list
data[1:4,1:4]
# ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457
# N1        4.936070      0.28569398        9.950503        1.425634
# N2        4.704400      0.00000000       12.708906        1.571140
# N3        6.968420      0.06803618       20.555438        1.249274
# N4        7.351349      0.03205427       17.829558        1.941931
table(data$group)
# Non-IBD      CD 
# 13      13 
# 绘图
dat_pca <- PCA(data[,-ncol(data)], graph = FALSE)  # 去除非数值列 pca结果
PCAplot<- fviz_pca_ind(dat_pca,
                       geom.ind = "point", # 只显示点，不显示文字
                       col.ind = data$group, # 用不同颜色表示分组
                       palette = c("#00AFBB", "#E7B800"),
                       addEllipses = T, # 是否圈起来，少于4个样圈不起来
                       legend.title = "disease state groups") + theme_bw()
PCAplot

# So I marked samples as PC1 score high and low (pink and green).
pca_group <- rep("PC1 high socres",26)
pca_group[which(dat_pca$ind$coord[,1]<0)] <- "PC1 low scores"
table(pca_group)
# PC1 high socres  PC1 low scores 
# 14              12 
PCAplot_byPCA1<- fviz_pca_ind(dat_pca,
                              geom.ind = "point", # 只显示点，不显示文字
                              col.ind = pca_group, # 用不同颜色表示分组
                              palette = c("#00AFBB", "#E7B800"),
                              addEllipses = T, # 是否圈起来，少于4个样圈不起来
                              legend.title = "PC1 scores groups") + theme_bw()
PCAplot_byPCA1

pca_group
mypca_group <- pca_group
save(mypca_group, group_list,file = "3.4.myPCA.RData")

其实可以明显看出来counts的PCA比FPKM的好一些

hierarchical clustering

rm(list=ls())
mydata <- read.table("known_gene_count_matrix.csv",sep = ',',header = T)
rownames(mydata) <- mydata$gene_id
mydata$gene_id <- NULL

load("1.merge_GSMdata_withGene.RData")
load("2.PCA.RData")
group_list
colnames(mydata) <- c(paste0("N",seq(1,13)),
                      paste0("C",seq(1,13)))
#表达矩阵数据校正
exprSet <- mydata
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
# 标准化 这里直接log2 
# ###limma----voom差异分析
# library(limma)
# #1
# exp=filter_count
# Group=group_list
# #2
# dge <- edgeR::DGEList(counts=exp)  # 需要借助edgeR包DGEList函数创建对象
# dge <- edgeR::calcNormFactors(dge)  # 以及对count标准化
# # 无比较矩阵
# #3
# design <- model.matrix(~Group)  # 所以类似edgeR的流程，也要单独设计方案
# #4
# v <- voom(dge,design, normalize="quantile")
exprSet <- log2(exprSet+1)
#差异分析：
dat <- exprSet
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
# 获取差异表达基因
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg)
# logFC AveExpr      t   P.Value adj.P.Val     B
# MSTRG.5743|MT1HL1 -4.968   4.318 -7.135 2.065e-07   0.01011 6.880
# MSTRG.22181|MT1G  -2.631  12.369 -6.594 7.494e-07   0.01297 5.580
# MSTRG.22176|MT1H  -3.797   9.012 -6.570 7.952e-07   0.01297 5.521
# MSTRG.22176|MT2A  -1.617  12.495 -6.400 1.201e-06   0.01469 5.105
# MSTRG.21460|ZP2   -3.606   2.371 -6.204 1.939e-06   0.01898 4.623
# MSTRG.7474|PRXL2A -1.477  11.128 -6.114 2.420e-06   0.01974 4.400

## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
if(T){
  logFC_t=2
  deg$g=ifelse(deg$P.Value>0.05,'stable',
               ifelse( deg$logFC > logFC_t,'UP',
                       ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
  )
  table(deg$g)
  # DOWN stable     UP 
  # 426  48232    288 
  head(deg)
}
# logFC AveExpr      t   P.Value adj.P.Val     B      g
# MSTRG.5743|MT1HL1 -4.968   4.318 -7.135 2.065e-07   0.01011 6.880   DOWN
# MSTRG.22181|MT1G  -2.631  12.369 -6.594 7.494e-07   0.01297 5.580   DOWN
# MSTRG.22176|MT1H  -3.797   9.012 -6.570 7.952e-07   0.01297 5.521   DOWN
# MSTRG.22176|MT2A  -1.617  12.495 -6.400 1.201e-06   0.01469 5.105 stable
# MSTRG.21460|ZP2   -3.606   2.371 -6.204 1.939e-06   0.01898 4.623   DOWN
# MSTRG.7474|PRXL2A -1.477  11.128 -6.114 2.420e-06   0.01974 4.400 stable

#####差异基因样本聚类热图------------

sigGenes <- rownames(deg[deg$g!="stable",])  # 714
# 用的是上一步通过差异分析得到的显著基因和tpm的结果
dat_deg <- dat[match(sigGenes,rownames(dat)),]
dat_deg[1:4,1:4]
# N1     N2     N3     N4
# MSTRG.5743|MT1HL1  6.238  7.942  7.723  8.553
# MSTRG.22181|MT1G  13.401 14.159 14.425 14.035
# MSTRG.22176|MT1H  10.439 11.601 11.553 11.601
# MSTRG.21460|ZP2    3.976  4.509  5.038  5.989
group_list
# 数据框形式分组
# disease
group_disease <- data.frame(group=group_list)
rownames(group_disease)=colnames(dat_deg)
library(pheatmap)
pheatmap(dat_deg,
         scale = "row",  # 标准化 非归一化
         show_colnames =T,show_rownames = F, 
         cluster_cols = T, 
         annotation_col=group_disease,
         main = "DEG_disease_cluster_rowScale")

# PC1
# pca_group
group_pc1 <- data.frame(group=pca_group)
rownames(group_pc1)=colnames(dat_deg)
library(pheatmap)
pheatmap(dat_deg,
         scale = "row",
         show_colnames =T,show_rownames = F, 
         cluster_cols = T, 
         annotation_col=group_pc1,
         main = "DEG_raw_PC1_cluster_rowScale")

save(mydata,file="3.3.mydata_deg.RData")
####
# 3.4
# 用counts进行差异分析PC1试一试
####
load("3.4.myPCA.RData")
mypca_group
group_pc1 <- data.frame(group=mypca_group)
rownames(group_pc1)=colnames(dat_deg)
group_pc1
library(pheatmap)
pheatmap(dat_deg,
         scale = "row",
         show_colnames =T,show_rownames = F, 
         cluster_cols = T, 
         annotation_col=group_pc1,
         main = "DEG_my_PC1_cluster_rowScale")

FPKM_PC1group

可以发现重新用其他上游分析软件获得的counts的聚类结果中确实也存在非线性系统误差，没有完全根据”state disease“分组聚类

但也没如作者所说的那样根据PC1得分分组聚类

---

总结

针对在一开始我们归纳的作者的两个观点，我们根据上面的复现结果进行讨论：

TPM/FPKM/RPKM并虽然理论上可以消除线性系统误差，但对差异分析而言结果却不如raw counts不理想，这在一定程度反映了这些均一化方法的局限性和处理带来的潜在偏差

我们的结果中counts相比FPKM拥有更好的PCA效果（除了FPKM处理，还可能的影响：我们自己进行的转录组上游分析软件所用参数，如质控可能更严格等），FPKM虽然理论上可以消除线性系统误差但是无法消除非线性系统误差，所以如果选择FPKM而非counts并不能达到理想的去除系统误差效果，并且这样的处理可能会掩盖掉数据的质量问题，目前主流也是根据counts来进行差异分析

在实验设计和操作环节需要尤其注意非线性误差的潜在影响，因为没有生物信息学方法可以处理非线性误差，任何不同来源的数据集的合并都是不可接受的

在我们的分析中可以发现该数据集确实存在非线性系统误差，并且经过自行使用其他软件进行上游分析得到counts进行差异分析，可以确定早在fastq文件中就存在这种误差，非线性系统误差来自实验设计、操作、测序等环节。

作者认为没有任何生物信息学方法可以在处理的过程中去除这些误差（这里做了分布形状的校准）从而进一步提出任何不同来源的数据集的合并都是不可接受的

此前表达量芯片专辑有篇推文介绍了如何整合不同数据集的结果老文新看，今天来看看两个数据集的整合分析，涉及到三种方法：

• 直接limma::removeBatchEffect()移除批次效应，合并数据集

• 分别对两个数据集进行差异分析，然后取它们的交集

• 在第二种方法基础上，RRA方法对多个排好序的差异基因集，进行求交集的同时还考虑一下它们的排序情况

这里主要谈论的是第一种方法的可行性

小编认为，①是否可以在“移除批次效应合并数据集”之前单独参照本文对非线性系统误差进行查看呢？如果各数据集实验分组聚类效果很好，是不是就可以合并呢？合并后再次查看实验分组聚类效果。

②对单个数据集，如果PCA的效果不好，转录组专辑此前也介绍了一些处理方法“亡羊补牢”：PCA效果不行，我们可以试试这样补救下，那么进行了补救（筛选聚类合格样本数据数据后）能否再进行多个数据集合并呢？

同时，我们的复现过程中存在一些跟作者不一致的地方也值得深入探讨，比如为什么作者预测认为存在非线性系统误差的数据集应该在PC1得分分组上有很好的聚类效果？欢迎大家在评论区讨论。