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