MingLi-929
commented
1 year ago 1. 测试heatmap中对矩阵的不同处理方式
1. 使用tidyr::pivot_longer()直接处理矩阵
2. 使用Matrix将矩阵转为dgRMatrix的稀疏矩阵---->提前生成全为0的矩阵----->根据稀疏矩阵得到行列index替换全为0的长矩阵(这个长矩阵与直接pivot_longer()转化的长矩阵相同)
3. 使用Matrix，将矩阵转为dgRMatrix的稀疏矩阵---->根据行列index直接生成长矩阵(这个长矩阵会比直接pivot_longer()转化的长矩阵减少，减少数值为0的)使用microbenchmark对三种方式进行测试
library("microbenchmark")
library(Matrix)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
files <- getSampleFiles()
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
mt <- getTagMatrix(peak = files[[5]],TxDb = txdb,
upstream = 1000,downstream = 1000,
type = "start_site",by = "gene",
weightCol = "V5")
mt <- t(apply(mt, 1, function(x) x/max(x)))
ii <- order(rowSums(mt))
mt <- mt[ii,]
colnames(mt) <- seq_len(dim(mt)[2])
rownames(mt) <- seq_len(dim(mt)[1])
## 第一种直接转化
mt_1<-mt
mt_1 <- mt_1 %>% as.data.frame() %>%
rownames_to_column("sample_ID") %>%
pivot_longer(-c(sample_ID),names_to = "coordinate",
values_to = "values")
mt_1$coordinate <- as.numeric(mt_1$coordinate)
mt_1$sample_ID <- as.numeric(mt_1$sample_ID)
## 第二种生成全0矩阵后进行替换
mt_2 <- as(mt, 'dgRMatrix')
mt_2_df <- data.frame(sample_ID=rep(seq_len(mt_2@Dim[1]),rep(mt_2@Dim[2],mt_2@Dim[1])),
coordinate=rep(seq_len(mt_2@Dim[2]),mt_2@Dim[1]),
values.1=rep(0,mt_2@Dim[1]*mt_2@Dim[2]))
col_pos <- mt_2@j+1
row_pos <- findInterval(seq(mt_2@x)-1,mt_2@p[-1])+1
df <- data.frame(sample_ID=row_pos,coordinate=col_pos,values=mt_sp@x)
merge_df <- merge(mt_2_df,df,all=TRUE,)[,c(1,2,4)]
merge_df[is.na(merge_df)] <- 0
## 第三种由稀疏矩阵直接生成不带有0的长矩阵
mt_3 <- as(mt, 'dgRMatrix')
col_pos <- mt_3@j+1
row_pos <- findInterval(seq(mt_3@x)-1,mt_3@p[-1])+1
mt_3_sp <- data.frame(sample_ID=row_pos,coordinate=col_pos,values=mt_3@x)
进行测试
# 对运行速度进行测试
func1 <- function(){
mt_1<-mt
mt_1 <- mt_1 %>% as.data.frame() %>%
rownames_to_column("sample_ID") %>%
pivot_longer(-c(sample_ID),names_to = "coordinate",
values_to = "values")
mt_1$coordinate <- as.numeric(mt_1$coordinate)
mt_1$sample_ID <- as.numeric(mt_1$sample_ID)
}
func2 <- function(){
mt_2 <- as(mt, 'dgRMatrix')
mt_2_df <- data.frame(sample_ID=rep(seq_len(mt_2@Dim[1]),rep(mt_2@Dim[2],mt_2@Dim[1])),
coordinate=rep(seq_len(mt_2@Dim[2]),mt_2@Dim[1]),
values.1=rep(0,mt_2@Dim[1]*mt_2@Dim[2]))
col_pos <- mt_2@j+1
row_pos <- findInterval(seq(mt_2@x)-1,mt_2@p[-1])+1
df <- data.frame(sample_ID=row_pos,coordinate=col_pos,values=mt_sp@x)
merge_df <- merge(mt_2_df,df,all=TRUE,)[,c(1,2,4)]
merge_df[is.na(merge_df)] <- 0
}
func3 <- function(){
mt_3 <- as(mt, 'dgRMatrix')
col_pos <- mt_3@j+1
row_pos <- findInterval(seq(mt_3@x)-1,mt_3@p[-1])+1
mt_3_sp <- data.frame(sample_ID=row_pos,coordinate=col_pos,values=mt_3@x)
}
bench_res <- microbenchmark("use tidyr"=func1(),
"use sparse matrix build long matrix with 0"=func2(),
"use sparse matrix build long matrix"=func3())
得到的测试结果
> bench_res
Unit: milliseconds
expr min lq mean
use tidyr 448.3222 462.9606 496.77029
use sparse matrix build long matrix with 0 4312.1044 4455.9427 4640.16852
use sparse matrix build long matrix 25.3467 29.2169 39.21166
median uq max neval
467.52880 541.4086 720.7547 100
4539.81870 4695.4377 5655.6305 100
30.37515 32.6257 171.1232 100得到如下结论
- 使用稀疏矩阵的index构建的长矩阵的速度最快，比，使用tidyr::pivot_longer()直接转化为长矩阵的速度的速度要快，快十倍。
- 先使用稀疏矩阵再转化为带有0的长矩阵的速度，比，直接使用tidyr::pivot_longer()的速度更慢，因此可以舍去这种方法。使用稀疏矩阵就直接构建不带有0的长矩阵。
2. 测试"处理矩阵"和"画图"的时间
根据上面可以知道
- 使用tidyr::pivot_longer()直接转化矩阵和直接使用稀疏矩阵的index两种方法是速度最快的
- 因此下面要对使用处理矩阵和进行画图分开测试
bench_res_tidyr <- microbenchmark(matrix=func1(),
plot=func4(mt_1))
bench_res_sparse <- microbenchmark(matrix=func3(),
plot=func4(mt_3_sp))测试结果如下所示
> bench_res_tidyr
Unit: milliseconds
expr min lq mean median uq max neval
matrix 466.1613 474.6662 507.455975 485.5752 524.41290 1134.3673 100
plot 2.3244 2.4761 2.907646 2.8150 2.98015 8.0449 100
> bench_res_sparse
Unit: milliseconds
expr min lq mean median uq max neval
matrix 26.3806 27.7055 34.810995 29.20545 35.9749 189.8263 100
plot 2.3634 2.5708 2.927292 2.89770 3.0372 5.6479 100
- 使用常规方法，画图时间占2.907646/(507.455975 + 2.907646)=0.57%
- 使用稀疏矩阵的方法，画图时间占2.927292/(34.810995 + 2.927292)=7.76%
可以得出结论：
-
进行绘制热图时候，准备矩阵的时间都在90%以上
-
使用稀疏矩阵进行处理矩阵能够大幅度减少运行时间
3. 使用稀疏矩阵直接构建的长矩阵绘图的差别
使用tidyr::pivot_longer()直接转化矩阵得到的热图
使用稀疏矩阵构建长矩阵的绘图
总结： 两者区别就在于，使用稀疏矩阵直接构建的长矩阵绘图没有给0配置颜色，直接显示空白。
GuangchuangYu
This pr change the way of plotting heatmap from base system to ggplot system. It adds the function of gradient color by giving
paletteparameter, which is the palette of gradient color.There are two ways to plot the heatmap. The first is giving a created matrix to
tagHeatmap()fucntion.peakHeatmap()to do all things in one steps.If users wanna to speed up the process, a binning method can be used by specifing the
nbin.Since this method use ggplot system, users can use ggplot method to change
a list of peaks can also be plotted.
plot body region
peakHeatmap(peak = files[[4]], TxDb = txdb, upstream = NULL, downstream = NULL, by = "gene", type = "body", nbin = 800)
extend body region with rel object
peakHeatmap(peak = files[[4]], TxDb = txdb, upstream = rel(0.2), downstream = rel(0.2), by = "gene", type = "body", nbin = 800)
extend body region with actual number
peakHeatmap(peak = files[[4]], TxDb = txdb, upstream = 1000, downstream = 1000, by = "gene", type = "body", nbin = 800)
txdb1 <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb2 <- unlist(fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene))
txdb <- list(gene19 = txdb1, geneX = txdb2) peakHeatmap_multiple_Sets(peak = files[[4]], upstream = 1000,downstream = 1000, by = c("gene19","geneX"), type = "start_site", TxDb = txdb,nbin = 800)
sample <- list(files[[4]],files[[5]]) peakHeatmap_multiple_Sets(peak = files[4:5], upstream = 1000, downstream = 1000, by = c("gene19","geneX"),type = "start_site", TxDb = txdb,nbin = 800)
txdb1 <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb2 <- unlist(fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene))[1:10000,]
txdb <- list(gene19 = txdb1, geneX = txdb2) peakHeatmap_multiple_Sets(peak = files[[4]], upstream = 1000,downstream = 1000, by = c("gene19","geneX"), type = "start_site", TxDb = txdb,nbin = 800)
peak_Profile_Heatmap(peak = files[[4]], upstream = 1000, downstream = 1000, by = "gene", type = "start_site", TxDb = txdb3, nbin = 800)
peak_Profile_Heatmap(peak = files[[4]], upstream = 1000, downstream = 1000, by = c("gene19","geneX"), type = "start_site", TxDb = txdb,nbin = 800)
peak_Profile_Heatmap(peak = files[4:5], upstream = 1000, downstream = 1000, by = c("gene19","geneX"), type = "start_site", TxDb = txdb,nbin = 800)