1. 下载Cisbp数据库中的motif文件
(1) 下载文件
CISBP
点击Bulk downloads进入到Bulk downloads界面
Bulk downloads
选择想对应的物种，然后点击Downloa Entire dataset Archive下载数据
(2) 理解文件
我下载下来的文件名字是Macaca_fascicularis_2024_09_12_9:10_am.zip
然后进行解压，解压完之后的文件如下所示
Macaca_fascicularis_2024_09_12
a. logos_all_motifs文件夹
logos_all_motifs文件夹放的是motif的图像
image.png
以M08124_2.00_fwd.png 和M09199_2.00_rev.png为例子
M08124_2.00_fwd.png：该文件表示 M08124_2.00 基序在正链（forward strand）上的图像或可视化结果。这是基序在正向方向上的表现。
M09199_2.00_rev.png：该文件表示 M09199_2.00 基序在反链（reverse strand）上的图像或可视化结果。它是基序在反向方向上的表现，通常是正链的互补序列。
我一共下载下来的是5899个motifs
cd Macaca_fascicularis_cisbp_20240912/logos_all_motifs
ls | grep rev | wc -l
# 5899
ls | grep png | wc -l
# 11798
ls | grep rev | wc -l
# 5899
b. pwms_all_motifs文件夹
pwms_all_motifs文件夹放的是motifs的pwms的文件,也是5899个文件
pwms_all_motifs
cat M09067_2.00.txt
M09067_2.00.txt
c.TF转录因子相关文件
TF_Information_all_motifs_plus.txt
TF_Information_all_motifs_plus
TF_Information_all_motifs.txt
TF_Information_all_motifs
TF_Information.txt
TF_Information
cat TF_Information_all_motifs_plus.txt | wc -l # 70015行 包含了TF和所有的motif
cat TF_Information_all_motifs.txt | wc -l # 70015行 包含了TF和所有的motif 感觉和TF_Information_all_motifs_plus.txt文件一样
cat TF_Information.txt | wc -l # 1340行 是和motif结合的转录因子的信息
这个1340行，我是从官网home页面进去后选择物种点击Go
image.png
到了这个页面点击download excel spreadsheet
image.png
(http://cisbp2.ccbr.utoronto.ca/downloads.php)
但解压后PWM.txt文件是空的，也不知道为什么
image.png
总之反正我摸索了下，得从Bulk downloads界面去下载你所需要的物种的motifs,那么解压下来就会有motif_pwm.txt文件了
2. 格式转换：将CisBP Motif文件处理为ArchR兼容的Motif PWMs格式
(1) 循环读取每个 PWM 文件构建PWMatrixLists文件
这一步我是参考这个
rm(list = ls())
# Macaca_fascicularis_cisbp_20240912/pwms_all_motifs从http://cisbp-rna.ccbr.utoronto.ca/bulk.php选择相应物种，然后download entire dataset archive
library(TFBSTools) # 参考https://blog.csdn.net/u012110870/article/details/102804617
motif_dir <- "Macaca_fascicularis_cisbp_20240912/pwms_all_motifs"
PWList <- list()
# 初始化跳过的文件计数器
skipped_files <- 0
# 循环读取每个 PWM 文件
for (file in list.files(motif_dir, pattern = ".txt", full.names = TRUE)) {
df <- read.table(file, header = TRUE, row.names = 1)
# 检查文件是否为空或是否有有效的碱基频率信息
if (nrow(df) == 0 || all(is.na(df))) {
message(paste("Skipping file:", file, "- No valid data"))
skipped_files <- skipped_files + 1 # 增加跳过文件计数
next
}
# 将数据转换为矩阵
mt <- as.matrix(df)
# 提取文件名作为 motif ID
motif_id <- substr(basename(file), 1, 6)
# 创建 PWM 对象并存储
PWList[[motif_id]] <- PWMatrix(ID = motif_id, profileMatrix = t(mt))
}
# 将 PWM 对象组合成 PWMatrixList
PWMatrixLists <- do.call(PWMatrixList, PWList) # 5899
# 打印结果
print(PWMatrixLists)
# 打印跳过的文件数量
message(paste("Total number of skipped files:", skipped_files)) # 919 # 一共是6818，减去919还剩5899
(2) 添加TF的信息
# 尝试使用 fill 参数读取文件
tf_info <- read.table(
"Macaca_fascicularis_cisbp_20240912/TF_Information_all_motifs_plus.txt",
header = TRUE,
sep = "\t",
stringsAsFactors = FALSE,
fill = TRUE # 自动填充不完整的行
)
tf_info$Motif_ID <- substr(tf_info$Motif_ID, 1, 6)
length(unique(tf_info$TF_ID)) # 1340
tf_info <- tf_info[tf_info$Motif_ID %in% names(PWMatrixLists), ]
length(unique(tf_info$TF_ID))
# [1] 846 过滤了, 有些没有PWMS文件
motif_ids <- names(PWMatrixLists)
# 将 profileMatrix 的列名设为 NULL
# tmp <- PWMatrixLists
# 循环遍历每个 PWM 对象并补充元数据
for (motif_id in motif_ids) {
# 提取对应的 motif 信息
motif_info <- tf_info %>% filter(Motif_ID == motif_id)
if (nrow(motif_info) > 0) {
# 提取相关信息
motif_name <- motif_info$TF_ID[1]
family_name <- motif_info$Family_Name[1]
species <- "Macaca fascicularis" # 统一物种信息
medline <- motif_info$PMID[1]
motif_type <- motif_info$Motif_Type[1]
collection <- motif_info$MSource_Identifier[1]
# 补充到 PWMatrix 对象中
PWMatrixLists[[motif_id]] <- PWMatrix(
ID = motif_id,
name = motif_name,
matrixClass = family_name,
profileMatrix = PWMatrixLists[[motif_id]]@profileMatrix, # 保留原始矩阵
strand = "+", # 假设为正链
tags = list(
comment = "Motif information added from TF_Information_all_motifs_plus.txt",
medline = medline,
type = motif_type,
collection = collection,
species = species
)
)
} else {
message(paste("No information found for motif:", motif_id))
}
}
# 检查更新后的 PWMatrixLists 因为pwm的colnames需要是NULL，所以多这一步
print(PWMatrixLists)
for (i in seq_along(PWMatrixLists)) {
# 将 profileMatrix 的列名设为 NULL
colnames(PWMatrixLists[[i]]@profileMatrix) <- NULL
}
(3) 去除上面列和不是1的motif 列
其实这个是因为要么原始文件行和加起来不是1，要么就是R中浮点精度的问题，是个坑
我是根据这个Issue 提示然后过滤掉这些error motifs的，我还没想到最终怎么去解决这些error motifs，先这样子
# 初始化一个空列表来存储不符合标准的 PWM 名称
error_list <- list()
# 遍历 PWMatrixLists 中的所有 PWM
for (i in seq_along(PWMatrixLists)) {
# 获取当前 PWM 的名字
pwm_name <- names(PWMatrixLists)[i]
# 获取当前 PWM 的 profileMatrix
profile_matrix <- PWMatrixLists[[i]]@profileMatrix
# 检查列的和是否为 1
colsum_check <- all.equal(colSums(as.matrix(profile_matrix)), rep(1, ncol(profile_matrix)))
# 如果不符合条件，将 PWM 名称加入到错误列表中
if (!isTRUE(colsum_check)) {
error_list[[pwm_name]] <- colsum_check
}
}
# 输出错误列表的 PWM 名称
print(names(error_list))
length(error_list)
PWMatrixLists_new <- PWMatrixLists[names(PWMatrixLists) %ni% names(error_list)] # %ni% 不在额返回true，就是提取不是error_list的PWMatrixLists
PWMatrixLists_new <- ArchR:::.summarizeJASPARMotifs(PWMatrixLists_new)
saveRDS(PWMatrixLists_new$motifs, "mcf_cisbp_PWMatrixLists_new_remove_errorlist.rds")
最终从原先的5899 motifs变成了5418 motifs
image.png
3. 在ArchR中使用addMotifAnnotations函数运行Motif PWMs文件
rm(list = ls())
setwd("~/analysis/20240914_atac_analysis")
library(ArchR)
library(motifmatchr)
library(BSgenome.Mfascicularis.Ensembl.mcf6)
addArchRThreads(threads = 12)
mcf_motifs <- readRDS("mcf_cisbp_PWMatrixLists_new_remove_errorlist.rds")
projHeme5 <- loadArchRProject(path = "Save-ProjHeme_MACS2_addpeaks/", force = FALSE, showLogo = TRUE)
print("projHeme5_addMotifAnnotations_mcf_motifs_500.rds")
projHeme5 <- addMotifAnnotations(
ArchRProj = projHeme5,
annoName = "mcf_motifs_500",
motifPWMs = mcf_motifs,
motifSet = NULL,
collection = NULL,
species = NULL
)
saveRDS(projHeme5, "projHeme5_addMotifAnnotations_mcf_motifs_allls.rds")
摸索了很久，解决了各种bug，终于总算可以自己构建成功了，可以用来跑archr中的addMotifAnnotations了，我真棒呀哈哈哈哈
