homer -- 通过差异基因寻找转录因子

看我写的教程不如看官网(\滑稽):http://homer.ucsd.edu/homer/

1. Installation

cd $HOME/software
mkdir homer
cd homer
wget http://homer.ucsd.edu/homer/configureHomer.pl
perl configureHomer.pl -install
export PATH=$HOME/software/homer/bin/:$PATH # 最好写到.bashrc里

2. Data donwload

perl $HOME/software/homer/configureHomer.pl -list # 列出可供下载的数据包,按需下载
# perl $HOME/software/homer/configureHomer.pl -install mouse
# perl $HOME/software/homer/configureHomer.pl -install mm8
perl $HOME/software/homer/configureHomer.pl -install hg19 # 我只需要人类的数据
perl $HOME/software/homer/configureHomer.pl -install human-p # 人类启动子数据

3. Find TF

目前我只用过findMotifs.pl这个功能,可以根据给定的基因ID或者fasta文件寻找已知的或潜在的转录因子结合区域,即motif,并给出对应的TF和统计数据。

养成好习惯,记得先 --help~

findMotifs.pl --help

两种用法:

  • findMotifs.pl <input list> <promoter set> <output directory> [additoinal options]
  • findMotifs.pl targets.fa fasta motifResults/ -fasta background.fa

如果已获得差异基因,直接使用第一种用法就行了,比较简单。

参数设置:

  • -start:TSS上游多少bp作为启动子区域,带负号,默认-300
  • -end:TSS下游多少bp作为启动子区域,不用带正号,默认50
  • 物种类型
  • 输出目录

其他的好像也不用怎么设置,功能都挺齐全的昂,GO、KEGGG、MsigDB富集分析、染色体区域、蛋白质相互作用以及我不了解的东东,略fancy~

示例:

测试数据:https://github.com/JiahaoWongg/Files/blob/main/homer_geneID.txt

鉴于github访问时好时坏。。。

ENSG00000271798.1
ENSG00000228350.1
ENSG00000169903.7
ENSG00000072163.19
ENSG00000242110.8
ENSG00000162639.16
ENSG00000133131.15
ENSG00000171617.14
ENSG00000258757.1
ENSG00000145604.16
ENSG00000100297.16
ENSG00000285160.1
ENSG00000011143.17
ENSG00000276043.5
ENSG00000117480.16
ENSG00000196260.5
ENSG00000185379.20
ENSG00000186185.14
ENSG00000129173.13
ENSG00000279347.1
ENSG00000127324.9
ENSG00000076770.14
ENSG00000186283.14
ENSG00000105427.10
ENSG00000273132.1
ENSG00000226065.1
ENSG00000122642.11

geneID最好用ENSEMBL,当然如果不是软件也会给你转换。

启动子范围我选择-800~200,另外需要指定文件输出目录,这里为out

findMotifs.pl homer_geneID.txt human out -start -800 -end 200

跑的还挺慢的吧,两个小时?

结果:

$ ls out
biocyc.txt              homerMotifs.motifs10  knownResults.html           prosite.txt
biological_process.txt  homerMotifs.motifs12  knownResults.txt            reactome.txt
cellular_component.txt  homerMotifs.motifs8   lipidmaps.txt               seq.autonorm.tsv
chromosome.txt          homerResults          molecular_function.txt      smart.txt
cosmic.txt              homerResults.html     motifFindingParameters.txt  smpdb.txt
gene3d.txt              interactions.txt      msigdb.txt                  wikipathways.txt
geneOntology.html       interpro.txt          pathwayInteractionDB.txt
gwas.txt                kegg.txt              pfam.txt
homerMotifs.all.motifs  knownResults          prints.txt

这里主要看homerResults文件夹,为预测的潜在的TF结合motif,至于已存在的TF结果knownResults我还没研究。

可以用浏览器打开homerResults.html进行查看。

不过,貌似并没有提供预测结果的汇总,我这里写了一个R函数,用于提取汇总homerResults里的结果:

subString <- function(strings, idx, sep = NA){

    strings = as.character(strings)
    if(is.na(sep)){
        res = as.character(lapply(strings, function(x) paste(strsplit(x, "")[[1]][idx], collapse = "")))
    } else{
        res = sapply(strsplit(strings, sep), function(x) x[idx])
    }
    return(res)
}

summaryHomer <- function(outFolder){

    homerFolder = paste0(outFolder, "/homerResults")
    xFiles = list.files(homerFolder, ".motif$")
    xFiles = xFiles[-grep("similar", xFiles)]
    xFiles = xFiles[-grep("RV", xFiles)]
    xFiles = xFiles[order(as.numeric(gsub("\\.", "", gsub("motif", "", xFiles))))]
    texts  = sapply(paste0(homerFolder, "/", xFiles), readLines)
    chunks = sapply(texts, function(x) strsplit(x[1], "[\t]"))

    motif = sapply(chunks, function(x) subString(x[1], 2, ">"))
    match = sapply(chunks, function(x) subString(subString(x[2], 2, "BestGuess:"),  1, "/"))
    score = sapply(chunks, function(x) rev(strsplit(x[2], "[()]")[[1]])[1])
    count = sapply(chunks, function(x) subString(x[6], 3, "[T:()]"))
    ratio = sapply(chunks, function(x) subString(x[6], 2, "[()]"))
    p_value = sapply(chunks, function(x) subString(x[6], 2, "P:"))

    xresT = data.frame(motif, 
                       match, 
                       score = as.numeric(score), 
                       count = as.numeric(count),
                       ratio_perc = as.numeric(gsub("%", "", ratio)), 
                       p_value = as.numeric(p_value)
                       )
    rownames(xresT) = gsub(".motif", "", basename(rownames(xresT)))
    return(xresT)
}

只需要提供homer的输出目录,我们来使用一下:

> summaryHomer(out)
               motif              match score count ratio p_value
motif1  TCAGGTGAAGGG         ZNF467(Zf) 0.635     4 0.07%    1e-8
motif2  YSGTGGAYTRYT            ZNF354C 0.746     3 0.04%    1e-7
motif3    TTCTVAATGC             BARHL1 0.558     8 2.98%    1e-7
motif4      TATTKTCC             NFATC2 0.772     9 6.70%    1e-5
motif5  CCTAACTGTTGG          BMYB(HTH) 0.640     2 0.02%    1e-5
motif6      TCKGTAGA               OSR1 0.704    10 9.54%    1e-5
motif7      AGCGGGAC          E2F6(E2F) 0.765     7 3.79%    1e-5
motif8      TGGCGATC               MZF1 0.707     8 6.11%    1e-4
motif9  TGAGCMGRGATA          GATA3(Zf) 0.605     3 0.26%    1e-4
motif10   GCGAGCTTGC              Nr2e3 0.613     4 0.85%    1e-4
motif11     ACCAGATT             TWIST1 0.723     3 0.40%    1e-4
motif12 CYGSYGCWWAMC    PB0112.1_E2F2_2 0.577     2 0.12%    1e-3
motif13 ACTGTGCCATTC    AR-halfsite(NR) 0.599     2 0.15%    1e-3
motif14   CCGTCATCAT                YY2 0.665     1 0.00%    1e-3
motif15   TCCTTATTCG              SOX14 0.620     1 0.00%    1e-3
motif16   GTGATCGGTA               CUX2 0.677     1 0.00%    1e-3
motif17   ATTGGCGTAT Oct2(POU,Homeobox) 0.703     1 0.00%    1e-3
motif18   CGATTGAATG              ZNF24 0.705     1 0.00%    1e-3
motif19   CGATGGCAGT           HIC1(Zf) 0.731     1 0.00%    1e-3
motif20   GACTGTATAG             ZBTB32 0.684     1 0.00%    1e-3
motif21   TTAGTAGCAC    PB0155.1_Osr2_2 0.705     1 0.00%    1e-3
motif22   AGTACGCGCT              HIF1A 0.640     1 0.00%    1e-3
motif23   CGAAGATGAT   POL008.1_DCE_S_I 0.608     1 0.00%    1e-3
motif24   ACTTATATGG                SRF 0.728     1 0.00%    1e-3
motif25   TTTATCGCAT               CDX2 0.758     1 0.00%    1e-3
motif26   TTTCTGTACG             ZBTB32 0.727     1 0.00%    1e-3
motif27   GGGTCGATCA   PB0030.1_Hnf4a_1 0.621     1 0.00%    1e-3
motif28   CATAAATTCG   PB0171.1_Sox18_2 0.706     1 0.00%    1e-3
motif29 CCTCAGTGTGTC             ZSCAN4 0.581     1 0.00%    1e-3
motif30 TAACCAGGATGT         SPDEF(ETS) 0.781     1 0.00%    1e-3
motif31 CTTGGGCAGTAC    PB0133.1_Hic1_2 0.679     1 0.00%    1e-3
motif32 GGGGGTTGGGTC               KLF4 0.681     1 0.00%    1e-3
motif33 TCCAAAAGAGCT             ZBTB26 0.613     1 0.00%    1e-3
motif34 AGCTGTGCTTTA              Gfi1b 0.732     1 0.00%    1e-3
motif35 CGTGAATGAAGC    PB0028.1_Hbp1_1 0.674     1 0.00%    1e-3
motif36 GTTTATTCTGGG              Foxq1 0.650     1 0.00%    1e-3
motif37 TCCCCACTGGAG               EBF1 0.679     1 0.00%    1e-3
motif38 TTTAGAGGAGGC  PB0203.1_Zfp691_2 0.629     1 0.00%    1e-3
motif39 GCTGGCATGTCT           HIC1(Zf) 0.724     1 0.00%    1e-3
motif40 GGCATGGGGCTC              CTCFL 0.671     1 0.00%    1e-3
motif41 CTGGGGCGCAGT         ZNF692(Zf) 0.648     1 0.00%    1e-3
motif42 AGAGGACGGCCG            YY1(Zf) 0.588     1 0.00%    1e-3
motif43 AAGCTTGGGAGG              Nr2e3 0.741     1 0.00%    1e-3
motif44 TCGCCCACGAGA           Egr2(Zf) 0.683     1 0.00%    1e-3
  • motif
    预测的motif序列,正链
  • match
    预测的TF
  • score
    匹配度,1为完全匹配
  • count
    你输入的基因中包含该motif的基因个数
  • ratio
    你输入的基因中包含该motif的基因个数占总输入基因个数的比例
  • p_value
    置信度

接下来,敲除?过表达?


后来又加了提取汇总knownResults里的结果的函数:

summaryHomerKnown <- function(outFolder){

    knownFolder = paste0(outFolder, "/knownResults")
    xFiles = list.files(knownFolder, ".motif$")
    xFiles = xFiles[order(as.numeric(gsub("\\.motif", "", gsub("known", "", xFiles))))]
    texts  = sapply(paste0(knownFolder, "/", xFiles), readLines)
    chunks = sapply(texts, function(x) strsplit(x[1], "[\t]"))

    motif = sapply(chunks, function(x) subString(x[1], 2, ">"))
    TF    = sapply(chunks, function(x) subString(x[2], 1, "/"))
    count = sapply(chunks, function(x) subString(x[6], 3, "[T:()]"))
    ratio = sapply(chunks, function(x) subString(x[6], 2, "[()]"))
    p_value = sapply(chunks, function(x) subString(x[6], 2, "P:"))

    xresT = data.frame(motif, 
                       TF, 
                       count = as.numeric(count),
                       ratio_perc = as.numeric(gsub("%", "", ratio)), 
                       p_value = as.numeric(p_value)
                       )
    rownames(xresT) = gsub("\\.motif", "", basename(rownames(xresT)))
    return(xresT)
}
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,902评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,037评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,978评论 0 332
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,867评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,763评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,104评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,565评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,236评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,379评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,313评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,363评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,034评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,637评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,719评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,952评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,371评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,948评论 2 341

推荐阅读更多精彩内容