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

2021-08-18 本文已影响0人生信摆渡

看我写的教程不如看官网(\滑稽）：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)
}

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

1. Installation

2. Data donwload

3. Find TF

猜你喜欢

热点阅读