1. 什么是样本关系探索

之前我们下载了5种阶段的柑橘,每个阶段取3个样本.15个样本之间关系,组间关系和组内关系,都在本次样本关系探索范围内.

2. 样本关系探索

2.1 导入样本信息表

我们回忆一下样本信息表是什么.最简单的样本信息表,只包含样本的分组信息.如果要更复杂的话,可以加上我们要研究的内容,比如重要、年龄、花青素含量等.

#读取表,并重命名列
sample_info = read.table('24.DE_analysis/sample_info.txt',col.names = c('group','sample'))
#也可以使用rename函数
read.table('24.DE_analysis/sample_info.txt') %>%
    rename(group = V1,sample = V2) -> sample_info

2.2 评估样本相似性

主要有三种方法:相关系数、距离矩阵、PCA主成分分析.

2.2.1 相关系数

2.2.1.1 相关系数计算

评估相似性也就是评估样本之间的关系.样本内有多个统计数据,我们可以通过相关系数衡量样本之间的相似性.那么评估相似性需要什么数据呢?是表达数据.也就是导入表达矩阵(注意是标准化后的).

genes_exp = read.table(file='23.merge/genes.TMM.TPM.matrix')
genes_exp

读入后,我们就要比对它们的相似性.计算相关系数这里有三种方法Pearson相关系数、Spearman相关系数,Kendall 相关系数.这里我们主要说Pearson相关系数.

Pearson相关系数
我们可以观察到下图样本 $1$ 和样本 $2$ 构成了两个向量.

我们用python把它画出来,长这样:

上图描述了X,Y两个样本的分布,我们可以计算相关系数.

对于线性、连续的数据来说,使用最多的就是Pearson相关系数.

在实际计算样本相似性前,我们先举个例子,计算相关系数.这个可以选择算具体的相关系数,比如Spearman.

geneA = c(2,3,7,8,10)
geneB = c(3,2,7,5,3)
geneC = c(39,76,NA,65,23)
#cor 计算两个向量之间的相关系数
cor(geneA,geneB)
cor(geneA,geneB,method='spearman')

#[1] 0.4054654
#[1] 0.3077935

也可以使用Python算,我直接用上面的X,Y数值.

np.corrcoef(x,y)
#array([[1.        , 0.99942998],
#       [0.99942998, 1.        ]])

如果要计算geneC的相关系数就会返回NA,因为该样本有个表达数据没有测到.如果要计算,就必须把NA删除,也就是将该样本所在行计算时也要删除.下面是use参数使用.

complete.obs:表示只使用完整的观测.
everything:表示使用所有观察数据
all.obs: 这一个参数要求数据中不能含有NA，否则报错，只有不含有NA才能正常计算
na.or.complete: 这个参数的作用和complete.obs类似，只不过如果每一列向量都有NA，则计算相关性的矩阵并不会报错,而会给出NA.
pairwise.complete.obs: 它的作用是如果某一列向量中有NA,那么计算该向量与其他向量相关性时,去除具有NA的那一行(如果其他向量相互计算没有NA,就完整计算).

cor(geneA,geneC)#NA
cor(geneA,geneC,use='complete.obs')#-0.3976854

如果我们想像Python一样得到矩阵,可以先转为dataframe.请注意,计算时是将样本Cs3删掉了,因为use='complete.obs'.

df_test = data.frame(
  row.names = c('Cs1', 'Cs2', 'Cs3', 'Cs4', 'Cs5'),
  geneA,geneB,geneC)
df_test
cor(df_test)
cor(df_test,use='complete.obs')

我们通过上面发现,cor用于计算列之间的相关系数,如果我们的样本是按行摆放的,就需要转置.

sample_cor = cor(genes_exp)
#保留两位小数
sample_cor = round(sample_cor,digits = 2)
sample_cor

我们使用%>%这个符号,一定要先引入tidyverse

2.2.1.2 绘制相关系数热图

如果没导入要先导入.这里直接根据矩阵的值转成颜色了.我们还可以看文档优化.

pheatmap(sample_cor)

cluster_cols: 列是否聚类(相似性高的放在一起)
cluster_rows: 行是否聚类(相似性高的放在一起)
cellwidth/height: 长宽.
border_color: 边界线颜色
font_size: 字体大小
angle_col: 将横纵坐标的字体弄斜.
display_number: 显示数值
fontsize_number: 数值大小

如果嫌颜色不好看,也可使用调色板Click.

library(RColorBrewer)
#数字表示从对应字符串的颜色序列中拿出多少个颜色
color = brewer.pal(5, "Paired")
#"#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C" "#FB9A99"
#我们可以先展示下颜色
library(scales)
show_col(color)

pheatmap(sample_cor,
         color = color,
         cluster_cols = F,cluster_rows = F,cellwidth = 15,cellheight = 15,display_numbers = T)

但是颜色只拿了5个,这里却有8个数字,所以需要多拿一些.颜色数量需要比数值多一个.

2.2.2 样本相关性-距离矩阵

2.2.2.1 距离矩阵计算

显然我们也可以使用距离衡量样本的相关性.dist()返回一个距离矩阵,可以计算行之间的距离.

sample_dist = dist(t(genes_exp))
sample_dist

2.2.2.2 层次聚类

计算完后,需要对样本之间进行聚类.需要使用hclust函数.但是我们看不到直接结果.关于层次聚类,可以看这个博客.不停迭代距离最近的点,然后构建一棵树,从实现过程来看,有点像最小生成树.


hclust(sample_dist)
#Call:
#hclust(d = sample_dist)
#
#Cluster method   : complete 
#Distance         : euclidean 
#Number of objects: 15

需要把它绘制出来.

plot(hclust(sample_dist))

2.2.3 主成分分析

2.2.3.1 主成分计算

就是在机器学习里学习的降维分析.从下面定义来说,是从一系列观测值按照不同的权重变成新的观测值,新的观测值的数量由我们决定.

主成分分析（英语：Principal components analysis,缩写：PCA）是一种统计分析、简化数据集的方法。它利用正交变换来对一系列可能相关的变量的观测值进行线性变换，从而投影为一系列线性不相关变量的值，这些不相关变量称为主成分（Principal Components）.

首先需要安装package.

BiocManager::install('PCAtools')
library(PCAtools)

在进行主成分分析之前,也需要对gene进行过滤.我们可以先去除一些表达量在不同样本中变化不大的基因.同时PCA分析要求metadata,文档描述为元数据,也就是样本信息表,它的行名必须等于mat的列名.主成分分析用到的是表达矩阵.

mat: 仅包含数字数据的数据矩阵或数据框.默认情况下,变量应位于行中,样本位于列中.
metadata: 包含元数据的数据矩阵或数据框.这些数据将存储在生成的 pca 对象中.必须遵守 rownames(metadata) == colnames(mat)

pca的返回值是一个pca类型的变量.loads为主成分负载,即线性变换中的加权数.具体如何,其实在官方文档有.

library(tidyverse)
sample_info_with_row = column_to_rownames(sample_info,var = 'sample')
pca_res = ppca(genes_exp,removeVar = 0.3,metadata= sample_info_with_row)
# removeVar: 去除30%的观测量,最好不超过0.5
# metadata: 样本信息表, 并且必须rownames(metadata) == colnames(mat)

建议看这个,补一个PCA原理.

2.2.3.2 PCA绘图

但是这输出结果不是我们想要的输出结果,我们可以将每个PCA的可解释度特征可视化.我个人理解这里是每个PCA成分为最后投影方向负责的比例.

screeplot(pca_res)

绘制PC1和PC2的二维图.把PCA1做横坐标,PCA2做纵坐标.一般来说,都是选择这两个做坐标轴.但如果要研究的东西对样本差异性可解释度低,那么就可能需要换PCA做坐标系.总的来说,PCA1区分开来不同的组,如果我们要研究不同组之间的差异,就需要让PCA1做坐标系.

biplot(pca_res,x='PC1',y='PC2')

下图看出Cs1等各个组都聚在一起,而Cs4和Cs5之间比较接近,所以它们之间的组间距离比较小.
如果需要按组分配颜色,就可以使用colby参数,还有legendPosition表示图例位置

biplot(pca_res,x='PC1',y='PC2',colby='group',legendPosition='top')

hline/vline: 添加两条横纵线,值是参照线的坐标,主要表示做个参照.
encircle: 主要是是否展示圆环,这里我个人理解是将相近的点用圆圈包起来.
ellipse: 主要是设置是否展示95%置信椭圆区间

有时候我们也会对genes_exp取 $log_2$ 然后再进行PCA分析.这是为了降低PCA特征之间的差异.但是有时候我们存在基因表达量为0的情况,这就不能直接取 $Log_2$ ,通常我们会再加 $1$ .但取不取 $log_2$ 看我们个人.如果 $log_2$ 效果不行的话,也可以尝试标准化与中心化.

3. 样本关系探索-目的

很多时候我们探索样本相关性,是为了锁定关键样本.比如我们探究油脂含量表达基因与生长周期是否有关.就是在高油和低油的作物,分别在不同生长周期采取样本,然后进行样本相关性分析(比如聚类分析),如果高油和低油的基因在10天后表达量就出现了差异,就需要着重分析10天左右(因为也可能是在第6、8等天)的差异.

4. appendix: PCA

PCA简单来说,是将x维的数据降到k维.在空间里找到一个坐标系,x维的数据经过映射变成k维,而且原数据的信息要尽可能的保存.这 $k$ 个维度我们称为主成分 $1$ ,主成分 $2$ 等等.以 $2$ 维数据举例最好理解,我们将黑点映射到黑色的直线上,同时距离(红色线长度)也要最小、映射到黑色线的点分散越开越好(也就是投影后方差最大).

PCA的关键就是找到新坐标系.找坐标系的关键在于原点和角度.原点很明显是均值,主要是找角度.但其实这个角度就是协方差矩阵的特征向量.建议这个看视频Click.

17. 样本关系探索