1. 下載Cisbp數(shù)據(jù)庫(kù)中的motif文件
(1) 下載文件
CISBP
點(diǎn)擊Bulk downloads進(jìn)入到Bulk downloads界面
Bulk downloads
選擇想對(duì)應(yīng)的物種，然后點(diǎn)擊Downloa Entire dataset Archive下載數(shù)據(jù)
(2) 理解文件
我下載下來(lái)的文件名字是Macaca_fascicularis_2024_09_12_9:10_am.zip
然后進(jìn)行解壓凤巨，解壓完之后的文件如下所示
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）上的圖像或可視化結(jié)果敢茁。這是基序在正向方向上的表現(xiàn)留美。
M09199_2.00_rev.png：該文件表示 M09199_2.00 基序在反鏈（reverse strand）上的圖像或可視化結(jié)果。它是基序在反向方向上的表現(xiàn)逢倍，通常是正鏈的互補(bǔ)序列景图。
我一共下載下來(lái)的是5899個(gè)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個(gè)文件
pwms_all_motifs
cat M09067_2.00.txt
M09067_2.00.txt
c.TF轉(zhuǎn)錄因子相關(guān)文件
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 感覺(jué)和TF_Information_all_motifs_plus.txt文件一樣
cat TF_Information.txt | wc -l # 1340行 是和motif結(jié)合的轉(zhuǎn)錄因子的信息
這個(gè)1340行亮蒋，我是從官網(wǎng)home頁(yè)面進(jìn)去后選擇物種點(diǎn)擊Go
image.png
到了這個(gè)頁(yè)面點(diǎn)擊download excel spreadsheet
image.png
(http://cisbp2.ccbr.utoronto.ca/downloads.php)
但解壓后PWM.txt文件是空的慎玖，也不知道為什么
image.png
總之反正我摸索了下趁怔，得從Bulk downloads界面去下載你所需要的物種的motifs,那么解壓下來(lái)就會(huì)有motif_pwm.txt文件了
2. 格式轉(zhuǎn)換：將CisBP Motif文件處理為ArchR兼容的Motif PWMs格式
(1) 循環(huán)讀取每個(gè) PWM 文件構(gòu)建PWMatrixLists文件
這一步我是參考這個(gè)
rm(list = ls())
# Macaca_fascicularis_cisbp_20240912/pwms_all_motifs從http://cisbp-rna.ccbr.utoronto.ca/bulk.php選擇相應(yīng)物種薪前，然后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()
# 初始化跳過(guò)的文件計(jì)數(shù)器
skipped_files <- 0
# 循環(huán)讀取每個(gè) 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 # 增加跳過(guò)文件計(jì)數(shù)
next
}
# 將數(shù)據(jù)轉(zhuǎn)換為矩陣
mt <- as.matrix(df)
# 提取文件名作為 motif ID
motif_id <- substr(basename(file), 1, 6)
# 創(chuàng)建 PWM 對(duì)象并存儲(chǔ)
PWList[[motif_id]] <- PWMatrix(ID = motif_id, profileMatrix = t(mt))
}
# 將 PWM 對(duì)象組合成 PWMatrixList
PWMatrixLists <- do.call(PWMatrixList, PWList) # 5899
# 打印結(jié)果
print(PWMatrixLists)
# 打印跳過(guò)的文件數(shù)量
message(paste("Total number of skipped files:", skipped_files)) # 919 # 一共是6818序六，減去919還剩5899
(2) 添加TF的信息
# 嘗試使用 fill 參數(shù)讀取文件
tf_info <- read.table(
"Macaca_fascicularis_cisbp_20240912/TF_Information_all_motifs_plus.txt",
header = TRUE,
sep = "\t",
stringsAsFactors = FALSE,
fill = TRUE # 自動(dòng)填充不完整的行
)
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 過(guò)濾了, 有些沒(méi)有PWMS文件
motif_ids <- names(PWMatrixLists)
# 將 profileMatrix 的列名設(shè)為 NULL
# tmp <- PWMatrixLists
# 循環(huán)遍歷每個(gè) PWM 對(duì)象并補(bǔ)充元數(shù)據(jù)
for (motif_id in motif_ids) {
# 提取對(duì)應(yīng)的 motif 信息
motif_info <- tf_info %>% filter(Motif_ID == motif_id)
if (nrow(motif_info) > 0) {
# 提取相關(guān)信息
motif_name <- motif_info$TF_ID[1]
family_name <- motif_info$Family_Name[1]
species <- "Macaca fascicularis" # 統(tǒng)一物種信息
medline <- motif_info$PMID[1]
motif_type <- motif_info$Motif_Type[1]
collection <- motif_info$MSource_Identifier[1]
# 補(bǔ)充到 PWMatrix 對(duì)象中
PWMatrixLists[[motif_id]] <- PWMatrix(
ID = motif_id,
name = motif_name,
matrixClass = family_name,
profileMatrix = PWMatrixLists[[motif_id]]@profileMatrix, # 保留原始矩陣
strand = "+", # 假設(shè)為正鏈
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 因?yàn)閜wm的colnames需要是NULL例诀，所以多這一步
print(PWMatrixLists)
for (i in seq_along(PWMatrixLists)) {
# 將 profileMatrix 的列名設(shè)為 NULL
colnames(PWMatrixLists[[i]]@profileMatrix) <- NULL
}
(3) 去除上面列和不是1的motif 列
其實(shí)這個(gè)是因?yàn)橐丛嘉募泻图悠饋?lái)不是1裁着，要么就是R中浮點(diǎn)精度的問(wèn)題拱她，是個(gè)坑
我是根據(jù)這個(gè)Issue 提示然后過(guò)濾掉這些error motifs的，我還沒(méi)想到最終怎么去解決這些error motifs桶雀，先這樣子
# 初始化一個(gè)空列表來(lái)存儲(chǔ)不符合標(biāo)準(zhǔn)的 PWM 名稱
error_list <- list()
# 遍歷 PWMatrixLists 中的所有 PWM
for (i in seq_along(PWMatrixLists)) {
# 獲取當(dāng)前 PWM 的名字
pwm_name <- names(PWMatrixLists)[i]
# 獲取當(dāng)前 PWM 的 profileMatrix
profile_matrix <- PWMatrixLists[[i]]@profileMatrix
# 檢查列的和是否為 1
colsum_check <- all.equal(colSums(as.matrix(profile_matrix)), rep(1, ncol(profile_matrix)))
# 如果不符合條件唬复，將 PWM 名稱加入到錯(cuò)誤列表中
if (!isTRUE(colsum_check)) {
error_list[[pwm_name]] <- colsum_check
}
}
# 輸出錯(cuò)誤列表的 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函數(shù)運(yùn)行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乍恐，終于總算可以自己構(gòu)建成功了测砂，可以用來(lái)跑archr中的addMotifAnnotations了，我真棒呀哈哈哈哈
