CHIP-seq流程学习笔记（12）-对特定gene使用bedtools getfasta后，用meme软件进行motif分析

最新推荐文章于 2025-05-30 16:29:19 发布

原创

最新推荐文章于 2025-05-30 16:29:19 发布 · 4.4k 阅读

18 ·

CC 4.0 BY-SA版权

这篇博客介绍了如何结合RNA-seq和ChIP-seq数据，利用bedtools提取基因启动子序列，并通过meme软件进行motif分析。博主分享了从gtf文件获取基因TSS信息、使用bedtools getfasta提取序列以及安装和使用meme的过程，强调了meme在de novo motif discovery和motif enrichment中的应用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

刚写题目时发现，现在的分析应该是RNA-seq与ChIP-seq杂合在一起的分析了，为了省事还接着在这个系列写吧。最近有站友发私信建议说把代码跑出来的结果放上，方便学习。这里说明一下：关于学习笔记我尽量把图片放上来，对于分析的记录，因为涉及到实验室的课题，为了避免不必要的麻烦就不放了，本意这个记录就是自己的试验记录，记录下代码过程，为以后学生发表文章参考。

1. 根据gene_list从hg19_genes.gtf文件中提取gene_TSS信息

也可以使用biomaRt包进行注释，中间出错尚未解决，有机会再尝试。

# 提取待分析的gene_list
# 读入待处理文件
my_file <- read.csv("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/significant_different_genes/significant_pvalue_different_genes_genecount_up.csv")
# 根据需要提取gene_symbol
gene_symbol <- my_file$X

# 处理基因组注释文件
# 读入特定基因组文件
hg19_genes_gtf <- read.csv("G:/hg19_genes.gtf", header = FALSE, sep = "\t")
# 对读入的数据某一列进行拆分，此处以空格为例
library(stringr)
# 设置空的"hg19_genes_gtf_V1"向量，其行数与待处理数据行数一致
hg19_genes_gtf_V1<-rep(NA,nrow(hg19_genes_gtf))
hg19_genes_gtf_V2<-rep(NA,nrow(hg19_genes_gtf))
hg19_genes_gtf_V3<-rep(NA,nrow(hg19_genes_gtf))
hg19_genes_gtf_V4<-rep(NA,nrow(hg19_genes_gtf))
# 利用for循环，对hg19_genes_gtf数据框中的第9列进行拆分提取，每一个拆分元素储存至一个变量中
for (i in 1:nrow(hg19_genes_gtf)){
   
   
  hg19_genes_gtf_V1[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[1]
  hg19_genes_gtf_V2[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[2]
  hg19_genes_gtf_V3[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[3]
  hg19_genes_gtf_V4[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[4]
}
# 对原数据框中的特定序列重新赋值
hg19_genes_gtf <- cbind(hg19_genes_gtf, hg19_genes_gtf_V1, hg19_genes_gtf_V2, hg19_genes_gtf_V3, hg19_genes_gtf_V4)

# 原始gtf文件中字符串末尾有分号故需要删除
# 此处以结尾最后一个字符为例，采用nchar()函数对字符串进行计数后减一，来代表最后一个字符。

# 设置空的"hg19_genes_gtf_V1"向量，其行数与待处理数据行数一致
# 此处NA可以用0代替
hg19_genes_gtf_V5<-rep(NA,nrow(hg19_genes_gtf))
# 利用for循环，对hg19_genes_gtf数据框中的第13列中的分号进行剔除，采用去掉末尾字符的方法
for (i in 1:nrow(hg19_genes_gtf)){
   
   
  hg19_genes_gtf_V5[i] <- substring(hg19_genes_gtf[i,13], 1, (nchar(as.character(hg19_genes_gtf[i,13]))-1))
  }
# 将编辑后的数值合并至原数据框中的特定序列
hg19_genes_gtf <- cbind(hg19_genes_gtf, hg19_genes_gtf_V5)

# 根据目的提取特定列进行保存
hg19_genes_gtf_done <- hg19_genes_gtf[, c(1, 4, 5, 7, 13)]

# 对得到的文件进行输出保存
write.csv(hg19_genes_gtf_done, "G:/hg19_genes_start_end.gtf", quote = FALSE, row.names = FALSE)

##利用match函数提取基因所在的行数并重新命名为row.NO文件