1. 读入数据：

rm(list=ls())
options(stringsAsFactors = F)
RPM <- read.table("GSE103788_filtered_RPM_genes.tsv.gz",header = T)
RC <- read.table("GSE103788_raw_counts_genes.tsv.gz",header = T)
rownames(RPM) <- RPM[,1]
RPM <- RPM[,-1]
rownames(RC) <- RC[,1]
RC <- RC[,-1]
RC[1:3,1:3]
#                         PH_WT_tumors_r1 PH_WT_tumors_r2 PH_WT_tumors_r3
#ENSMUSG00000000001_Gnai3            1827            2772            2149
#ENSMUSG00000000003_Pbsn                0               0               0
#ENSMUSG00000000028_Cdc45              44              88              72
RPM[1:3,1:3]
#                         PH_WT_tumors_r1 PH_WT_tumors_r2 PH_WT_tumors_r3
#ENSMUSG00000000001_Gnai3      111.459868      149.657246      120.341735
#ENSMUSG00000000028_Cdc45        2.684310        4.751024        4.031924
#ENSMUSG00000000031_H19          1.159134       26.994453        2.183959

2.表达矩阵探究

参考：RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同

Counts值

对给定的基因组参考区域，计算比对上的read数，又称为raw count（RC）。
计数结果的差异的影响因素：落在参考区域上下限的read是否需要被统计，按照什么样的标准进行统计。

RPM (Reads per million mapped reads)

RPM方法：10^6标准化了测序深度的影响，但没有考虑转录本的长度的影响。
RPM适合于产生的read读数不受基因长度影响的测序方法，比如miRNA-seq测序，miRNA的长度一般在20-24个碱基之间。

先来研究一下两个表达矩阵

dim(RC)
#[1] 41128    12
dim(RPM)
#[1] 17311    12

参考公众号文章给出的公式，利用raw_counts矩阵构建RPM矩阵，然后跟文章提供的RPM矩阵比较

rr <- apply(RC,2,function(x){
  tmp <- x*10E6/sum(x)
  return(tmp)
})
dim(rr)
#[1] 41128    12
rr[rownames(RPM)[1:5],1:3]
#                         PH_WT_tumors_r1 PH_WT_tumors_r2 PH_WT_tumors_r3
#ENSMUSG00000000001_Gnai3     1114.598680     1496.572460     1203.417347
#ENSMUSG00000000028_Cdc45       26.843099       47.510237       40.319241
#ENSMUSG00000000031_H19         11.591338      269.944527       21.839589
#ENSMUSG00000000037_Scml2        3.050352        5.398891        6.159884
#ENSMUSG00000000049_Apoh     19231.250248    12456.860165    21013.604440
RPM[1:5,1:3]
#                         PH_WT_tumors_r1 PH_WT_tumors_r2 PH_WT_tumors_r3
#ENSMUSG00000000001_Gnai3     111.4598680     149.6572460     120.3417347
#ENSMUSG00000000028_Cdc45       2.6843099       4.7510237       4.0319241
#ENSMUSG00000000031_H19         1.1591338      26.9944527       2.1839589
#ENSMUSG00000000037_Scml2       0.3050352       0.5398891       0.6159884
#ENSMUSG00000000049_Apoh     1923.1250248    1245.6860165    2101.3604440

肉眼可见，自己计算得到的RPM矩阵比作者给出的RPM矩阵大10倍。。。

3. 样本分组探究

3.1 相关性热图

RPM_M <- cor(RPM)
pheatmap::pheatmap(RPM_M,main = 'cor(RPM)')

PH_WT和L_TAA基本上是分开了，
除了PH_WT_tumors_r1和PH_WT_tumors_r3乱入到了L_TAA分组

3.2 样本PCA分析

colnames(RC)
# [1] "PH_WT_tumors_r1"  "PH_WT_tumors_r2"  "PH_WT_tumors_r3"  "PH_WT_notreat_r1"
# [5] "PH_WT_notreat_r2" "PH_WT_notreat_r3" "L_TAA_24_r1"      "L_TAA_24_r2"     
# [9] "L_TAA_48_r1"      "L_TAA_48_r2"      "L_TAA_control_r1" "L_TAA_control_r2"
gl_1 <- c(rep("PH_WT_tumors",3),
                rep("PH_WT_notreat",3),
                rep("L_TAA",4),
                rep("L_TAA_control",2))
gl_2 <-  c(rep("PH_WT",6),rep("L_TAA",6))
dat <- t(RPM)
dat <- as.data.frame(dat)
dat <- cbind(dat,gl_1)
suppressMessages(library(FactoMineR))
suppressMessages(library(factoextra))
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca,
             # 只显示点，不显示文字
             geom.ind = "point", 
             # 用不同颜色表示分组
             col.ind = dat$gl_1, 
             #palette = c("#00AFBB", "#E7B800"),
             # 是否圈起来
             addEllipses = TRUE,
             legend.title = "Groups"
)

PCA - 4groups.png

后3组的样本数过少，没有办法圈起来

dat <- t(RPM)
dat <- as.data.frame(dat)
dat <- cbind(dat,gl_2)

dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca,
             # 只显示点，不显示文字
             geom.ind = "point", 
             # 用不同颜色表示分组
             col.ind = dat$gl_2, 
             palette = c("#00AFBB", "#E7B800"),
             # 是否圈起来
             addEllipses = TRUE,
             legend.title = "Groups"
)