1. 准备数据

使用BatMeth2这个软件生成的结果，将不同样本之间的甲基化合并到一个文件中，最后生成文件merged.Region.mC.txt

2. 脚本

library(ggsci)
library(tidyr)
library(dplyr)
library(stringr)

meth <- read.delim("D:/槜李/结果_图/PCA/merged.Region.mC.txt") %>% #导入文件
mutate(methRatio = Methylated_Count/Total_Count) #增加一列，计算甲基率

region_meth_cpg = filter(meth) %>%
dplyr::select(-c(Methylated_Count, Total_Count)) %>%
pivot_wider(names_from = sample, values_from = methRatio) %>%
mutate(id = str_c(Chromosome, Position, Strand, Methylation_Type, sep = "_")) %>% #文件中Chromosome, Position, Strand, Methylation_Type列的名称
dplyr::select(id, starts_with('ZL')) %>%
column_to_rownames(var = "id") #这一步生成的region_meth_cpg文件格式在下面

meth_pca = pca(region_meth_cpg,
removeVar = 0.5,
metadata = sample_info) #sample_info是以sample当作列名的文件，导入这个文件时rowname = 1

rownames_to_column(meth_pca$rotated, var = "sample") %>%
left_join(sample_info, by = "sample") %>%
ggplot(aes(x = PC1, y = PC2)) +
geom_point(shape = 21, size = 8, aes(fill = group)) +
scale_fill_aaas() +
labs(fill = NULL) +
theme_bw() +
theme(axis.title = element_text(size = 16),
legend.direction = "horizontal",
legend.position = "top")

region_meth_cpg的文件格式，中间存在NA值也没有影响

3. 结果