超级增强子系列9: ggseqlogo进行Motif文件可视化

library(ggseqlogo)
library(ggplot2)
library(dplyr)
library(purrr)
library(gridExtra)
read_motif_matrix <- function(file) {
  ext <- tolower(tools::file_ext(file))
  if (ext == "transfac") {
    # ----------------------------
    # TRANSFAC 格式
    # ----------------------------
    lines <- readLines(file)
    pos_lines <- lines[grepl("^[[:space:]]*[0-9]+", lines)]
    if (length(pos_lines) == 0) stop("No position data in TRANSFAC file: ", file)
    df <- read.table(text = pos_lines, header = FALSE, sep = "", stringsAsFactors = FALSE)
    if (ncol(df) < 5) stop("Insufficient columns in TRANSFAC file: ", file)
    mat <- t(as.matrix(df[, 2:5]))  # 转置为 4×L
    rownames(mat) <- c("A", "C", "G", "T")
    return(mat)
  } else if (ext == "meme") {
    # ----------------------------
    # MEME 格式
    # ----------------------------
      lines <- readLines(file)
      lines <- lines[!grepl("^\\s*$", lines) & !grepl("[A-Za-z]", lines)]
      mat_long <- as.matrix(read.table(text = lines, sep = "", header = FALSE))
      mat <- t(mat_long)  # L×4 → 4×L
      rownames(mat) <- c("A", "C", "G", "T")
  } else if  (ext == "jaspar") {
    # ----------------------------
    # JASPAR 格式（修正版）
    # ----------------------------
    lines <- readLines(file)
    # 跳过第一行（header）
    data_lines <- lines[2:5]
    df <- read.table(text = data_lines, sep = "", header = FALSE, stringsAsFactors = FALSE)
    # 第一列是碱基（A/C/G/T），其余为计数；跳过任何非数字列（如含 "[" 的列）
    # 找出从第2列开始、能转为 numeric 的列
    numeric_cols <- sapply(df[-1], function(x) {
      !any(is.na(suppressWarnings(as.numeric(as.character(x)))))
    })
    if (!all(numeric_cols)) {
      # 如果有非数字列（如 "["），只保留能转 numeric 的列
      count_cols <- which(numeric_cols) + 1  # +1 因为 df[-1]
    } else {
      count_cols <- 2:ncol(df)
    }
    mat <- as.matrix(df[, count_cols])
    rownames(mat) <- c("A", "C", "G", "T")
  } else if (ext == "pfm") {
    # ----------------------------
    # PFM 格式（原始）
    # ----------------------------
    mat <- as.matrix(read.table(file, skip = 1, fill = TRUE))
    if (nrow(mat) != 4) {
      mat <- as.matrix(read.table(file, fill = TRUE))
    }
    if (nrow(mat) != 4) stop("PFM file must result in 4 rows (A,C,G,T): ", file)
    rownames(mat) <- c("A", "C", "G", "T")
  } else if (ext == "txt") {
    # ----------------------------
    # TXT 格式（PWM，带 A:/C: 行名）
    # ----------------------------
    mat_df <- read.table(
      file,
      skip = 1,
      row.names = 1,
      sep = "\t",
      check.names = FALSE,
      stringsAsFactors = FALSE,
      fill = TRUE
    )
    mat <- as.matrix(mat_df)
    cleaned_rownames <- gsub("[^ACGT]", "", rownames(mat))
    expected <- c("A", "C", "G", "T")
    if (!setequal(cleaned_rownames, expected)) {
      stop("File ", file, " missing valid rows. Got: ", paste(cleaned_rownames, collapse = ", "))
    }
    rownames(mat) <- cleaned_rownames
    mat <- mat[expected, , drop = FALSE]
  } else {
    stop("Unsupported file extension: '", ext, "' for file: ", basename(file))
  }
  if (any(is.na(mat))) warning("NA values introduced in matrix from: ", file)
  return(mat)
}
# 匹配所有支持的格式
all_files <- list.files(
  pattern = "\\.(pfm|txt|transfac|meme|jaspar)$",
  full.names = TRUE,
  ignore.case = TRUE
)
# 批量读取
motif_names <- tools::file_path_sans_ext(basename(all_files))
motif_list <- map(all_files, read_motif_matrix)
# 可视化
plot_list <- map2(motif_list, motif_names, ~{
  ggseqlogo(.x, method = "prob") +
    theme_classic() +
    labs(title = .y, x = "Position", y = "Probability") +
    theme(plot.title = element_text(size = 9, face = "bold"),
          axis.text = element_text(size = 7))
})
# 5. 创建PDF文件
pdf("all_motifs_col.pdf", width = 8.5, height = 11)
# 每页2列，每列6个图（共12个图）
ncol <- 2
nrow <- 6
plots_per_page <- ncol * nrow
# 分页输出
for (i in seq(1, length(plot_list), by = plots_per_page)) {
  end_idx <- min(i + plots_per_page - 1, length(plot_list))
  # 获取当前页的图
  page_plots <- plot_list[i:end_idx]
  # 如果最后一页图不够，用空白填充
  if (length(page_plots) < plots_per_page) {
    blank_count <- plots_per_page - length(page_plots)
    for (j in 1:blank_count) {
      page_plots[[length(page_plots) + 1]] <- ggplot() + theme_void()
    }
  }
  # 排列当前页
  grid.arrange(
    grobs = page_plots,
    ncol = ncol,
    nrow = nrow,
    top = paste("Page", ceiling(i/plots_per_page), 
                "- Showing motifs", i, "to", end_idx)
  )
}
dev.off()
cat("✓ 已创建 all_motifs_col.pdf\n")
cat("✓ 包含", length(plot_list), "个motif\n")
cat("✓ 布局: 每页", ncol, "列 ×", nrow, "行 =", plots_per_page, "个motif\n")

 输出图片如下