基于motif序列筛选含该motif的ChIP-Seq峰及靶基因的算法需求

一、将HOMER Motif转换为正则表达式 #

HOMER输出的motif含简并碱基（如[C/T/G]），直接转换为正则的字符类格式即可：

• 去掉分隔符/，将可选碱基放在方括号[]内，比如[C/T/G]转为[CTG]，[T/C]转为[TC]
• 核心保守序列直接保留，比如CAAGG直接写在正则中

示例：若目标motif规则是第一位可选C/T/G，第二位可选T/C，核心为CAAGG，对应的正则表达式为：

motif_regex <- "[CTG][TC]AAGG"

二、筛选含目标Motif的峰序列 #

使用R的grepl()函数匹配正则表达式，从df中筛选符合条件的峰：

# 示例数据框
df <-data.frame(
  position=c( "chr1,1165251,1165500","chr1,2147751,2148000"),
  SEQ=c("CAGCACCAGGAACAGTGCGGCCGTGAGCTCAGCCTTCAAGGTCAAGCTCCTCACCTCCCTGTCAGAGCCCTGGCTCCAGCGT",
        "AAGAGGGGAAAAGACCATAGCAGAGAGCTCAGCCTGAGCTCCTCACCTCCCTGTCAGAGCCACGG")
)

# 筛选含motif的峰
filtered_peaks <- df[grepl(motif_regex, df$SEQ, ignore.case = TRUE), ]
# 输出结果与示例一致
print(filtered_peaks$position)

三、峰位置与基因列表的重叠分析 #

使用Bioconductor的GenomicRanges包高效处理基因组区间重叠：

1. 准备数据格式 #

拆分峰的位置信息为染色体、起始、终止，并转换为GRanges对象：

if (!require("GenomicRanges")) {
  BiocManager::install("GenomicRanges")
  library(GenomicRanges)
}

# 拆分峰位置
filtered_peaks_split <- do.call(rbind, strsplit(filtered_peaks$position, ","))
colnames(filtered_peaks_split) <- c("chr", "start", "end")
filtered_peaks_split <- transform(filtered_peaks_split, 
                                  start = as.numeric(start), 
                                  end = as.numeric(end))

# 转为GRanges对象
peaks_gr <- GRanges(
  seqnames = filtered_peaks_split$chr,
  ranges = IRanges(start = filtered_peaks_split$start, end = filtered_peaks_split$end)
)

2. 处理基因列表并寻找重叠 #

将基因列表（需包含chr、start、end、gene_name）转为GRanges对象，然后寻找重叠区间：

# 示例基因列表（替换为你的真实数据）
gene_list <- data.frame(
  gene_name = c("GeneX", "GeneY"),
  chr = c("chr1", "chr1"),
  start = c(1165000, 2147000),
  end = c(1166000, 2149000)
)

# 转为GRanges对象
genes_gr <- GRanges(
  seqnames = gene_list$chr,
  ranges = IRanges(start = gene_list$start, end = gene_list$end),
  gene_name = gene_list$gene_name
)

# 寻找重叠的靶基因
overlaps <- findOverlaps(peaks_gr, genes_gr)
target_genes <- gene_list[subjectHits(overlaps), ]

# 合并峰与靶基因信息
final_result <- cbind(filtered_peaks[queryHits(overlaps), ], target_genes)

关键说明 #

• 正则匹配时ignore.case=TRUE可兼容序列中的大小写碱基
• GenomicRanges支持多种重叠规则（如部分重叠、完全包含），可通过findOverlaps()的type参数调整（默认是任意重叠）
• 如果motif是反向互补的，需同时匹配正向和反向互补序列，可生成互补正则后用|拼接（如motif_regex <- "[CTG][TC]AAGG|CCTTG[GA][AC]"）

内容的提问来源于stack exchange，提问作者mehdi heidari