1.导入三张表（表达矩阵、样本信息表、基因信息表）
1.1 表达矩阵
都是有两种导入方法，这里只交代代码导入
read.table("path/to/you/genes.TMM.EXPR.matrix", header = T, row.names = 1)
#header = T表示存在表头，row.names = 1表示以第1行作为表头
1.2 样本信息表
awk -F "_" '{print $0"\t"$1"\t"$2}' name.lst > sample_info.txt
#awk命令生成样品信息表
sample_info <- read.delim("path/to/your/sample_info.txt")
#read.delim导入样本信息表
需要注意的是，样本信息表中的行名一定一定要与表达矩阵中的列名完全一致，否则在PCA分析时就会报错：Error in pca(gene_exp, metadata = sample_info, removeVar = 0.1) :
'colnames(mat)' is not identical to 'rownames(metadata)'
1.3 基因信息表
图形界面导入
import datasetfrom text(readr)
Browse
First row as Names(F)
Delimiter(Tab)
comment(#)
代码导入
library(readr)
library(tidyverse)
gene_info <- read_delim("path/to/the/query_seqs.fa.emapper.annotations",
delim = "\t", escape_double = FALSE,
col_names = FALSE, comment = "#", trim_ws = TRUE) %>%
select(Gene_Id = X1, Gene_Symbol = X6, GO = X7, Ko = X9, Pathway = X10, COG = X21, Gene_Name = X22)
代码解读：
readr包是为了增强数据读取能力，tidyverse包用于数据处理和可视化
readr包中的read_delim函数，用于从文本文件中读取数据并创建数据框
delim = "\t"：指定数据文件中的字段分隔符为制表符（Tab）
escape_double = FALSE：不对双引号进行转义
col_names = FALSE：指定文件没有列名，因为数据文件中的第一行没有列名
comment = "#"：指定以 "#" 开头的行为注释，这些行将被忽略
trim_ws = TRUE：在读取数据时删除字段周围的空白
%>%：R语言中的管道符
select()函数：是tidyverse包中的一个函数，用于选择和重新命名数据框中的列
2.样本相关性计算
2.1有关代码
sample_cor <- round(cor(gene_exp, method = 'spearman'), digits = 2)
library(pheatmap)
pheatmap(sample_cor)
cor(gene_exp, method = 'spearman'):cor()函数计算样本相关性，计算方法为spearman
round(x, digits = 2)：round()函数进行四舍五入处理，结果保留两位有效数字
library(pheatmap)：加载pheatmap包
pheatmap(sample_cor)：可视化处理
2.2 几种有关函数
pearson:适用于线性相关分析，可用于计算两样本间相关性
spearman:适用于等级变量但无线性关系的数据，如肿瘤的分期
kendall:适用于离散型、分类变量
3. 样本聚类
3.1 计算距离矩阵
gene_dist <- dist(t(gene_exp))
代码解读
t(gene_exp):t()是转置符，将行变成列，将列变成行
dist():dist()函数计算距离矩阵并将结果保存到gene_dist向量中，若未指明其他参数，则默认矩阵模型为欧几里得(euclidean)矩阵。
3.2 层次聚类
gene_hc <- hclust(gene_dist)
plot(gene_hc)
代码解读
通过hclust()函数计算上一步距离矩阵的结果从而得到层次聚类结果，并通过plot()函数可视化。
4.主成分分析
4.1加载分析包
library(PCAtools)
4.2 主成分分析
p <- pca(gene_exp, metadata = sample_info, removeVar = 0.1)
代码解读
gene_exp：是前面导入的基因信息表
metadata：就是样本信息表
removeVar = 0.1表示过滤波动值小于10%的基因
pca_loadings <- p$loadings
pca_rotated <- p$rotated
代码解读
PC1=ag1+ag2+ag3+....+xgx:PC1与每个基因都有关系，但关系有大有小，a就是贡献值的大小也称作主成分复杂(loadings)
screeplot(p)
代码解读
绘出来的图是主成分对样本的解释度，不同的主成分能对样本拉开较大差异。注意：PCA的数目由样本数决定，有多少个样本数就有多少个PCA。
biplot(p, x = 'PC1',
y = 'PC3',
colby = 'strain',
colkey = c('KID' = 'red', 'BLO' = 'yellow'),
shape = 'stage',
legendPosition = "bottom")
plotloadings(p)
代码解读
p是主成分分析中保存的向量p，x表示X轴，一般是由PC1作为x轴，y表示Y轴，一般是由PC1其他作为x轴，colby表示对样本信息表中某列赋予颜色，colkey表示对某列的不同组别赋予不同颜色（我们的研究没得必要区分颜色），shape表示对某列不同值赋予不同的形状，legendPosition 表示将图列放在哪边（上：top下：bottom左：left右：right）。
根据自己的研究可以不断更换y轴的输入，我们研究的是不同时期的差异，那就是意味着我们需要不同形状能得到区分的PC，找出这个PC阶段后，使用plotloadings(p)函数查看哪些基因对PC贡献大？
5.保存与加载数据
save(gene_exp, gene_info, sample_info, file = 'path/to/your/rna-seq/prepare.rdata')
load(file = 'path/to/your/rna-seq/prepare.rdata')
save和load是一对命令，在rstudio中每次退出会询问是否要保存数据，但在命令行中可能就不会询问，如果后续还需要继续分析，则需要重新加载向量，这时我们就先将结果保存到prepare.rdata中，再使用时直接加载即可。
6.完整流程的R脚本
#导入数据
##表达矩阵
gene_exp <- read.table("D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/genes.TMM.EXPR.matrix",
header = T, row.names = 1)
##样本信息表
sample_info <- read.delim("D:/share/R_data/data/rnaseq-apple/sample_info.txt", header = T,
row.names=1)
##基因信息表
library(readr)
library(tidyverse)
gene_info <- read_delim("D:/share/R_data/data/rnaseq-apple/query_seqs.fa.emapper.annotations",
delim = "\t", escape_double = FALSE,
col_names = FALSE, comment = "#", trim_ws = TRUE) %>%
select(Gene_Id = X1,
Gene_Symbol = X6,
GO = X7,
Ko = X9,
Pathway = X10,
COG = X21,
Gene_Name = X22,)
#样本相关性分析
sample_cor <- round(cor(gene_exp, method = 'spearman'), digits = 2)
library(pheatmap)
pheatmap(sample_cor)
#样本聚类分析
##计算距离矩阵
gene_dist <- dist(t(gene_exp))
##层次聚类
gene_hc <- hclust(gene_dist)
plot(gene_hc)
#主成分分析
library(PCAtools)
##主成分分析（PCA）
p <- pca(gene_exp, metadata = sample_info, removeVar = 0.1)
pca_loadings <- p$loadings
pca_rotated <- p$rotated
screeplot(p)
biplot(p, x = 'PC1',
y = 'PC3',
colby = 'strain',
colkey = c('KID' = 'red', 'BLO' = 'yellow'),
shape = 'stage',
legendPosition = "bottom")
plotloadings(p)
#保存与加载数据
save(gene_exp, gene_info, sample_info, file = 'path/to/your/rna-seq/prepare.rdata')
load(file = 'path/to/your/rna-seq/prepare.rdata')
网友评论