1. 下载Cisbp数据库中的motif文件

(1) 下载文件

CISBP数据库

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了，我真棒呀哈哈哈哈