单细胞ChIP-seq 样本聚类【TF-IDF/LSI】

2020-05-24 本文已影响0人 caokai001

缘由：

之前翻译了一篇单细胞ChIP文章单细胞ChIP-seq技术（CoBATCH）
看到有这个图，大概意思：将H3K27ac 位点及其细胞分别进行了层次聚类，使用HOMER进行每个enhancer module de nove TF, 使用GREAT进行GO富集分析。使用二项分布进行P-Value。

image.png

以为就是行是很多peak 位置，列是不同细胞的count matrix聚类热图，用pheatmap 等等自动画出来了，但是看文章的方法部分，有点晕，LSI score 进行降维聚类???下面列出图解和method 部分截图。

image.png

image

我理解图E做法是，选出20112 个common peak 区域(Tips：没搞懂作者如何选取的这些peak位点，将所有细胞bam 合并，用macs2 callpeak -broad 大概可以获得20万个peak，猜测选取q 值最小的20112 peak区域) ; 列是不同的细胞。

Step1 : 想将矩阵转换成binary matrix (文章多次提到需要将peak-cell matrix 进行binary ,可能是保证单个细胞单个位点DNA 覆盖度只能为0/1 ；和sc-RNA-seq 可能不一样)
Step2: 进行TF-IDF 进行权重校正,好像是自然语言处理里对一篇文章中，每一个字出现的概率进行校正，当这个字在多个文章中出现越不重要，当这个字在这一篇文章中，多次出现越重要。
Step3: 进行SVD 降维
Step4: 对细胞进行聚类和可视化

下面使用测试数据，完成上面操作：

### Step1：建立测试数据
dat <- matrix(sample(x = 0:2,size = 1000,replace=T,prob = c(0.5,0.3,0.2)),
  nrow = 100)
dat <- dat[apply(dat, 1, sum)!=0 ,]
dim(dat)
rownames(dat) <- paste0("peak_",1:100)
colnames(dat) <- paste0("cell_",1:10)
head(dat)

> head(dat)
       cell_1 cell_2 cell_3 cell_4 cell_5 cell_6 cell_7 cell_8 cell_9 cell_10
peak_1      1      1      0      0      0      0      0      1      0       1
peak_2      1      0      0      1      0      0      1      1      0       1
peak_3      0      0      1      0      1      0      1      0      0       0
peak_4      1      0      0      0      0      0      0      1      1       1
peak_5      1      1      0      0      0      1      1      1      1       0
peak_6      1      0      1      1      1      0      0      1      1       0

### Step2： binary matrix
dat[dat > 0] = 1
View(dat)

### Step3:TF-IDF
# 参考：https://blog.csdn.net/kMD8d5R/article/details/88564616

### 去掉0 的peak_num
library(tidyverse)
DAT <- dat
f_table <- as.data.frame(DAT) %>% mutate("peak_num"=rownames(.)) %>% 
  gather(key=cell,value=count,-peak_num)  %>% select(c(2,1,3)) %>%  
  filter(count!=0)

### Step4：加载需要包
library(pacman)
p_load(tidytext,rio,gtools)

# 参考：自然排序https://stackoverflow.com/questions/2778039/how-to-perform-natural-sorting-in-r


### Step5: 计算tf-idf
tf_idf_table <- f_table %>% 
  bind_tf_idf(term = peak_num, document = cell,n=count)
head(tf_idf_table)

### Step6 ： 转换成peak-cell weight matrix
tf_idf_table_df <- tf_idf_table %>% select(c(1,2,6)) %>% spread(key = peak_num,value =tf_idf,fill=0 ) %>% 
  as.data.frame()  
rownames(tf_idf_table_df) <- tf_idf_table_df$cell
## 转置
tf_idf_table_df <- tf_idf_table_df[,-1] %>%  t() 

## 对行列排序
tf_idf_table_df <- tf_idf_table_df[mixedorder(rownames(tf_idf_table_df) ),]
tf_idf_table_df <- tf_idf_table_df[,mixedorder(colnames(tf_idf_table_df))] 

head(tf_idf_table_df)

### Step7: SVD 降维
# 参考：R语言SVD运用， https://blog.csdn.net/u013259893/article/details/40483189
# svd 输入矩阵，或者数据框
inputdata <- tf_idf_table_df
svdresults <- svd(inputdata, nu = min(nrow(inputdata), ncol(inputdata)),
                  nv = min(nrow(inputdata), ncol(inputdata)), LINPACK = F)

# Tips : 特征值10个，细胞存在10维坐标
cell_coor <- svdresults$v  ## 10个细胞降维坐标
colnames(cell_coor) <- paste0("cell_",1:10)

### Step8: 层次聚类
# 参考：https://zhuanlan.zhihu.com/p/28264290
str(cell_coor)
# 按行进行标准化，消除量纲，计算样本间距离
scale_cell_coor <- as.matrix(cell_coor) %>% t() %>% scale() %>% t()

hc <- hclust(dist(scale_cell_coor ,method = "euclidean"),method = "ward.D2")
hc

# 画聚类树状图
plot(hc,hang = -0.01,cex=0.7)