两个基因相关性计算出来跟文章不一样有什么内幕吗？

学员提问如下：

学员的绘图结果 vs 文献中的结果

主要差别还是坐标轴范围差异太大了，点的分布稍微不太一样，相关性值差异也比较大，到底是什么原因呢！！！

数据背景

做这个图还是得要了解一下文章背景，数据为：.cgi?acc=GSE120643。总共28个样本：

GSM3406869 WT 0hpf rep1
GSM3406870 WT 0hpf rep2
GSM3406871 WT 2hpf rep1
GSM3406872 WT 2hpf rep2
GSM3406873 WT 3hpf rep1
GSM3406874 WT 3hpf rep2
GSM3406875 WT 4hpf rep1
GSM3406876 WT 4hpf rep2
GSM3406877 WT 6hpf rep1
GSM3406878 WT 6hpf rep2
GSM3406879 MO 3hpf rep1
GSM3406880 MO 3hpf rep2
GSM3406881 MO 4hpf rep1
GSM3406882 MO 4hpf rep2
GSM3406883 MO 6hpf rep1
GSM3406884 MO 6hpf rep2
GSM3428297 MOpabpc1a 4hpf rep1
GSM3428298 MOpabpc1a 4hpf rep2
GSM3428299 MOpabpc1a 6hpf rep1
GSM3428300 MOpabpc1a 6hpf rep2
GSM3732410 ribo-amanitin_2hpf_rep1
GSM3732411 ribo-amanitin_2hpf_rep2
GSM3732412 ribo-amanitin_4hpf_rep1
GSM3732413 ribo-amanitin_4hpf_rep2
GSM3732414 ribo-amanitin_6hpf_rep1
GSM3732415 ribo-amanitin_6hpf_rep2
GSM3732416 ribo-WT_4hpf_rep1
GSM3732417 ribo-WT_4hpf_rep2

实验设计如下

文献研究对象为斑马鱼的胚胎发育，专业术语 The maternal-to-zygotic transition (MZT)：

母体-合子转换（MZT）是指在胚胎发育早期，从母体提供的遗传信息和调控机制向合子自身基因组激活和自主调控的转变过程。这一过程是胚胎发育的关键阶段，确保胚胎能够从母体的“启动”状态过渡到自主发育的状态。

• Zygotic period
• Cleavage period
• Blastula period
• Gastrula period

HPF（Hours Post Fertilization）：指的是从受精卵形成开始计算的时间，以小时为单位。例如，24 HPF 表示受精后24小时。

文章做了非常多的实验：

YBX1 morphants embryos：指的是通过 morpholino 抗义寡核苷酸（morpholino antisense oligonucleotides, MO） 技术敲低（knockdown）YBX1 基因表达的胚胎。Morpholino 是一种化学合成的寡核苷酸，能够特异性结合 mRNA 的特定序列，阻止其翻译或导致 mRNA 降解，从而实现基因敲低。

Control embryos：对照胚胎，通常是没有经过基因敲低处理的正常胚胎，用于与 YBX1 morphants 进行比较，以评估基因敲低对胚胎转录组的影响。

通过 RNA-seq 比较 YBX1 morphants 和 control embryos 的转录组差异，研究 YBX1 基因在胚胎发育中的功能。

Ribo-amanitin 是一种从毒鹅膏菌（Amanita phalloides）中提取的毒素，属于 α-amanitin 的一种形式，具有强烈的毒性。它主要通过抑制真核生物的 RNA 聚合酶 II 的活性来发挥作用，从而干扰细胞内的转录过程。

• Wide type (Danio rerio)：正常野生型胚胎不同时期
• ybx1 MO (Danio rerio)：YBX1 基因敲除后不同时期
• pabpc1a MO (Danio rerio)：pabpc1a 基因敲除后不同时期
• Wide type (Danio rerio)：amanitin处理

复现的图背景

这个图使用的是 Wide type (Danio rerio)：正常野生型胚胎不同时期 转录组测序 数据中 ybx1 与 pabpc1a这两个基因在10个样本中的相关性散点图，值为RPKM，数据文章提供在 附件 mmc2.xlsx 的sheet h. RPKM。

Fig.4D Scatterplots showing the normalized RPKM of ybx1 and pabpc1a from two replicates. The Pearson correlation coefficient (R) and p values are shown in the top left corner.

画一下图看看

理论上 这个图有了表达矩阵很容易就做出来了：

rm(list=ls())
library(readxl)

# 读取数据
data <- read_excel("mmc2.xlsx", sheet = "h. RPKM", skip = 1)
data <- data[,1:11]
data <- as.data.frame(data)
colnames(data)
data[1, ]
colnames(data) <- paste0(colnames(data),"_", data[1,])
colnames(data) <- c("ID","WT-0hpf_Rep1","WT-0hpf_Rep2","WT-2hpf_Rep1","WT-2hpf_Rep2",
                    "WT-3hpf_Rep1","WT-3hpf_Rep2","WT-4hpf_Rep1","WT-4hpf_Rep2",
                    "WT-6hpf_Rep1","WT-6hpf_Rep2")
data <- data[-1, ]
head(data)
data[1,2]
class(data[1,2])

# 转为数值
df <- type.convert(data)
str(df)
head(df)

## 添加symbol
library(org.Dr.eg.db) #斑马鱼org.Dr.eg.db
keytypes(org.Dr.eg.db)
library(clusterProfiler)
id2symbol <- bitr(df$ID, 
                  fromType = "ENSEMBL", 
                  toType = "SYMBOL", 
                  OrgDb = org.Dr.eg.db)
head(id2symbol)

df <- merge(id2symbol,df, by.x="ENSEMBL",by.y="ID",all.y=T)
head(df)

提取两个基因的表达

# (D) Scatterplots showing the normalized RPKM of ybx1 and pabpc1a from two replicates. 
# The Pearson correlation coefficient (R) and p values are shown in the top left corner.
# 从rpkm_raw2中把ybx1、pabpc1a相关的10个样本中的表达量复制粘贴下来
library(ggplot2)
#数据
df1 <- df[df$SYMBOL %in% c("ybx1", "pabpc1a"), ]
rownames(df1) <- df1$SYMBOL
df1 <- log2(df1[,-c(1,2)])
df1

data <- t(df1)
str(data)
data

# 计算相关性
cor_result <- cor.test(as.numeric(data[,1]), as.numeric(data[,2]), method = "pearson")
cor_result
R_value <- round(cor_result$estimate, 2)
p_value <- signif(cor_result$p.value, 3)

ggplot2绘图

#画散点图
ggplot(data, aes(x = ybx1, y = pabpc1a)) +
  geom_point(color = "brown", size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", color = "black", se = FALSE) +
  annotate("text", x = 7.2, y = 8.8, label = paste0("R = ", R_value, "\n p = ", p_value), size = 5, hjust = 0) +
  labs(x = expression(Log[2]~"(RPKM) ybx1"),
       y = expression(Log[2]~"(RPKM) pabpc1a")) +
  #scale_x_continuous(limits = c(0, 5)) +
  #scale_y_continuous(limits = c(0, 5)) +
  theme_bw(base_size = 14)

结果如上，跟文章中的坐标范围也完全不一样，R相关性差别也比较大！

看了学员的代码，她除了 log2的时候+1对数据有一些影响外，还做了坐标轴范围限定（可能是为了想跟文献对比吧！）。

那为什么数值差别这么大呢？用的也是文献提供的rpkm值。

找到问题

我反复查看了文章的数据分析方法，发现里面用了一个词语被我忽略了：Scatterplots showing the normalized RPKM of ybx1 and pabpc1a from two replicates。

也就是说文献用的标准化后的rpkm值，标准化方法。文章在做RNA-seq测序的时候加了ERCC内参，然后对RPKM进行了矫正：

问题来了：ERCC内参矫正效果如何？是转录组测序必须的吗？

我们下期来看看 如何做 ERCC内参矫正以及 矫正前后的差异比较~（未完待续）