引言

❝
pyscenic 分析除了需要我们自己的表达矩阵文件，此外还需要 转录因子名称，转录因子到 motif 的注释文件 和 基因的 motif 富集信息 三个文件。官网只提供了 人，小鼠，果蝇 这三个物种的信息文件。对于其它物种，我们进行 pyscenic 分析仍然是一个挑战。官网提供了一个如何构建 cisTarget databases 的教程，我们分享一下如何构建上面提到的三个文件。

需要 fasta 和 motif 文件：

一次性构建需要三个文件：

参数解释：

安装 create_cisTarget_databases

github：https://github.com/aertslab/create_cisTarget_databases

克隆到本地，创建环境：
git clone https://github.com/aertslab/create_cisTarget_databases

Create conda environment.

conda create -n create_cistarget_databases \
‘python=3.10’ \
‘numpy=1.21’ \
‘pandas>=1.4.1’ \
‘pyarrow>=7.0.0’ \
‘numba>=0.55.1’ \
‘python-flatbuffers’

install Cluster-Buster precompiled binary

Activate conda environment.

conda activate create_cistarget_databases

pip install numpy==1.22.4
安装 Cluster-Buster：
cd “${CONDA_PREFIX}/bin”

Download precompiled Cluster-Buster binary.

wget https://resources.aertslab.org/cistarget/programs/cbust

Make downloaded binary executable.

chmod a+x cbust
安装 UCSC 的 liftOver，bigWigAverageOverBed：

Install some UCSC tools

cd “${CONDA_PREFIX}/bin”

Download liftOver.

wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver

Download bigWigAverageOverBed.

wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/bigWigAverageOverBed

Make downloaded binaries executable.

chmod a+x liftOver bigWigAverageOverBed
下载转录因子-motif 信息

❝
JASPAR 是一个收集了转录因子结合位点 motif 的数据库。motif，即转录因子结合位点结合的序列特征，转录因子通常通过识别 motif，与基因的调控区域结合，发挥转录调控作用。该数据库包含了6 大类物种的信息，分别是Fungi（真菌）、Insecta（昆虫）、Nematoda（线虫）、Plantae（植物）、Urochordata（尾索动物）和 Vertebrata（脊椎动物）。 网址：https://jaspar.elixir.no/

❝
如果是做植物的，可以使用 PlantTFDB 数据库下载相关信息。地址：https://planttfdb.gao-lab.org/

❝
首先确定你研究的物种属于哪个大的种属里面，再去下载 非冗余的 meme 格式 的文件，左侧的 Individual PFMs 则是每个 motif 为单个文件，这里我们下载脊椎动物的：

文件内容就是 每个 motif 及对应的转录因子及 PFM 矩阵：
wget https://jaspar.elixir.no/download/data/2024/CORE/JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt

less -S JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt

MEME version 4

#

ALPHABET= ACGT

#

strands: + –

#

Background letter frequencies

A 0.25 C 0.25 G 0.25 T 0.25

#

MOTIF MA0004.1 Arnt

letter-probability matrix: alength= 4 w= 6 nsites= 20 E= 0

0.200000 0.800000 0.000000 0.000000

0.950000 0.000000 0.050000 0.000000

0.000000 1.000000 0.000000 0.000000

0.000000 0.000000 1.000000 0.000000

0.000000 0.000000 0.000000 1.000000

0.000000 0.000000 1.000000 0.000000

URL http://jaspar.genereg.net/matrix/MA0004.1

#

MOTIF MA0069.1 PAX6

letter-probability matrix: alength= 4 w= 14 nsites= 43 E= 0

0.046512 0.093023 0.093023 0.767442

0.046512 0.046512 0.000000 0.906977

0.093023 0.604651 0.023256 0.279070

0.906977 0.046512 0.023256 0.023256

0.069767 0.790698 0.023256 0.116279

0.023256 0.000000 0.953488 0.023256

0.023256 0.860465 0.093023 0.023256

0.488372 0.046512 0.046512 0.418605

0.023256 0.093023 0.023256 0.860465

0.046512 0.325581 0.581395 0.046512

0.837209 0.000000 0.139535 0.023256

0.255814 0.255814 0.302326 0.186047

0.023256 0.116279 0.069767 0.790698

0.023256 0.000000 0.395349 0.581395

URL http://jaspar.genereg.net/matrix/MA0069.1

…

❝
我们需要对下载的数据进行处理，得到每个 motif 对应一个矩阵文件，需要以.cb 结尾，文件名和 motif 保持一致：

提取motif及矩阵

grep -E “MOTIF|^[[:space:]][0-9]” JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt | sed ‘s/MOTIF />/g’ |sed ‘s/^[[:space:]]//g’ > tf_motif_pfm_maxtrix.txt

替换空格

sed -i ‘s/>MA([0-9.]*) />MA\1_/’ tf_motif_pfm_maxtrix.txt

check

less tf_motif_pfm_maxtrix.txt

>MA0004.1_Arnt

0.200000 0.800000 0.000000 0.000000

0.950000 0.000000 0.050000 0.000000

0.000000 1.000000 0.000000 0.000000

0.000000 0.000000 1.000000 0.000000

0.000000 0.000000 0.000000 1.000000

0.000000 0.000000 1.000000 0.000000

>MA0069.1_PAX6

0.046512 0.093023 0.093023 0.767442

0.046512 0.046512 0.000000 0.906977

0.093023 0.604651 0.023256 0.279070

0.906977 0.046512 0.023256 0.023256

0.069767 0.790698 0.023256 0.116279

0.023256 0.000000 0.953488 0.023256

0.023256 0.860465 0.093023 0.023256

0.488372 0.046512 0.046512 0.418605

0.023256 0.093023 0.023256 0.860465

0.046512 0.325581 0.581395 0.046512

0.837209 0.000000 0.139535 0.023256

0.255814 0.255814 0.302326 0.186047

0.023256 0.116279 0.069767 0.790698

0.023256 0.000000 0.395349 0.581395

保存每个文件：

输出保存到每个文件

mkdir motif_dir
cd motif_dir
awk ‘/^>/{if(file) close(file); filename=substr($0,2)”.cb”; print $0 > filename; file=filename; next} {print >> file}’ ../tf_motif_pfm_maxtrix.txt

ls

MA0002.3_Runx1.cb MA0515.1_Sox6.cb MA0755.2_CUX2.cb MA1116.2_RBPJ.cb MA1623.2_Stat2.cb

MA0003.5_TFAP2A.cb MA0516.3_SP2.cb MA0756.3_ONECUT2.cb MA1117.2_RELB.cb MA1624.2_Stat5a.cb

MA0004.1_Arnt.cb MA0517.2_STAT1::STAT2.cb MA0757.2_ONECUT3.cb MA1118.2_SIX1.cb MA1625.2_Stat5b.cb

MA0006.2_Ahr::Arnt.cb MA0518.2_Stat4.cb MA0758.1_E2F7.cb MA1119.2_SIX2.cb MA1627.2_Wt1.cb

MA0007.4_Ar.cb MA0519.2_Stat5a::Stat5b.cb MA0759.3_ELK3.cb MA1120.2_SOX13.cb MA1628.2_Zic1::Zic2.cb

MA0009.2_TBXT.cb MA0520.2_Stat6.cb MA0760.2_ERF.cb MA1121.2_TEAD2.cb MA1629.2_Zic2.cb

❝
准备每个 motif id 文件：
grep “>” tf_motif_pfm_maxtrix.txt|sed ‘s/>//g’ > motifs_id.txt

head motifs_id.txt

MA0004.1_Arnt

MA0069.1_PAX6

MA0071.1_RORA

MA0074.1_RXRA::VDR

MA0101.1_REL

MA0107.1_RELA

MA0115.1_NR1H2::RXRA

MA0116.1_Znf423

MA0119.1_NFIC::TLX1

MA0130.1_ZNF354C

准备启动子序列文件

❝
前面提到需要的fasta序列（regions/genes），可以是转录因子结合区域的序列，也可以自己提取每个基因上游的启动子区域的序列，这里我们手动去提取基因转录上游3K的序列。
下载猪的基因组和注释文件（https://www.ensembl.org/Sus_scrofa/Info/Index）：

download genome and gtf

wget https://ftp.ensembl.org/pub/release-113/gtf/sus_scrofa/Sus_scrofa.Sscrofa11.1.113.gtf.gz

wget https://ftp.ensembl.org/pub/release-113/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz
计算正负链基因上游3000bp的位置：
library(rtracklayer)
library(tidyverse)
library(Biostrings)

gtf <- import(“Sus_scrofa.Sscrofa11.1.113.gtf.gz”,format = “gtf”) %>%
as.data.frame()

get gene promoters regions

genes <- gtf %>%
dplyr::filter(type == “gene” & gene_name != “NA”) %>%
dplyr::select(seqnames, start, end, strand, gene_name) %>%
dplyr::mutate(start2 = ifelse(strand == “+”,start – 3000,end + 1),
end2 = ifelse(strand == “+”,start – 1,end + 3000))

check

head(genes)

seqnames start end strand gene_name start2 end2

1 1 226161299 226217308 – ALDH1A1 226217309 226220308

2 1 226414306 226432279 + ANXA1 226411306 226414305

3 1 227574648 227772671 + RORB 227571648 227574647

4 1 227811204 227963309 – TRPM6 227963310 227966309

5 1 228033996 228045013 – C9orf40 228045014 228048013

6 1 228054915 228092263 – CARNMT1 228092264 228095263

对基因组染色体重新命名：

load genome file

gemome <- readDNAStringSet(filepath = “Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz”)
gemome

DNAStringSet object of length 613:

width seq names

[1] 274330532 GCTTAATTTTTGTCATTTCTCACCCCTGCTCTTGA…CTTTCTTTTTTCTTGAGACTTACGTTCATTGGGGA 1 dna:primary_ass…

[2] 151935994 GAAGTGCCCTGCCAAGTGGCCTAAGTCTTTCTCTC…GGTTGGTTAGGGTTAGGGTTAGATCTTCTCTCTTT 2 dna:primary_ass…

[3] 132848913 AGAAATGTAGAAGAAATGGGTGATTTGTAATAAAA…GGTTAGGTAGGGTTAGGGTTAGGGTTAGGGTTAGG 3 dna:primary_ass…

[4] 130910915 CTCCTTTGACTCTTGGTCTAAATTAGCTCACCTGA…TAGGGCTGCTTCCTGTGGCCTATATCTCTCTTCTT 4 dna:primary_ass…

[5] 104526007 ACACGGCTCGGGGATGGGTGCGGCGGCTCCGCGTC…TGATGGTACGAACGTCGTGAAGGTACAAGGAACTA 5 dna:primary_ass…

… … …

[609] 15496 TCTACTTTTCTCCCTCAATACTATGAGTAAAGGTA…AGATGTGGTACATATACACAATGGAATATTACTCG AEMK02000270.1 dn…

[610] 15455 GTAAGTACTAATCTCTAGACTCCCTGGTGCATGTG…GAAGTTGTTTACGAAAGACAAGGGAGAATACACAA AEMK02000323.1 dn…

[611] 15107 GAGGTTTAAACTTAACAAATAACTTCAGAATGTGT…ACACACAAAAAATAAACAGAAAAAGTACACGCTTC AEMK02000417.1 dn…

[612] 15107 GACGGCCCGTGGGGTCGACCATATGTCCCGGGCGA…CCTGGGTGTCGAACAGCTGTCCCAGAGGACTCACA AEMK02000654.1 dn…

[613] 15096 CCCTCTGAACACTCACTGCCAGCCTTTTTTTTTTT…TTGTGGATTTAAATGGTGAATTCATGGCATTTAAA AEMK02000519.1 dn…

rename seqnames

names(gemome) <- sapply(strsplit(names(gemome),split = ” “),”[“,1)

head

head(names(gemome),25)

[1] “1” “2” “3” “4” “5” “6”

[7] “7” “8” “9” “10” “11” “12”

[13] “13” “14” “15” “16” “17” “18”

[19] “X” “Y” “MT” “AEMK02000452.1” “AEMK02000698.1” “AEMK02000361.1”

[25] “AEMK02000598.1”

循环提取启动子序列:

extract promter sequences

x = 1

lapply(1:nrow(genes),function(x){
print(x)
tmp <- genes[x,]

# try extract seq
out <- tryCatch(
seq <- gemome[[tmp$seqnames]][tmp$start2:tmp$end2],
error = function(e){ return(NULL) }
)

# check
if(!is.null(out)){
# strand
if(tmp$strand == “-“){
out <- reverseComplement(out)
}else{
out
}

}else{
DNAString()
}

}) %>% DNAStringSet() -> seq_list

assign names and process unique gene name

names(seq_list) <- make.unique(genes$gene_name)

seq_list

DNAStringSet object of length 17633:

width seq names

[1] 3000 AGTGATCAGAAGCAAAGGGATATGTCCATTGCTGGG…CATGCACATAAAAAGGAACAAATAAAGCTAAGTGCT ALDH1A1

[2] 3000 TAATAATAACTATTACTCAAATATTTTTATGCAATG…AAACTCTATAAAATGAGAAGCCCAAGGCTCCACTTC ANXA1

[3] 3000 GACATTAAGAGTTCTTATATTTTTCTGGCTTTCCCA…CGGGGCGATGATTGGCGGCGACGCTGACAGGAGGCG RORB

[4] 3000 CATCCAACTAAGTTTTTAATGAATGCTCGTTATGTT…AACAAGAACAATAAAATATGTCTGTGGTTATGATGA TRPM6

[5] 3000 GCAAAAGAATCAAACTGGACTACTTTCTCATACCAC…TATCGGTGCTCGGCCTCCCGCGCCGGGCGGGCTCCG C9orf40

… … …

[17629] 3000 CAGTATCCCTATTTTGAAAAGAAAGAAAAGGTGTGG…TGTGTGTATAAATACATGAATGAAAATAAAATAGTC U6

[17630] 3000 CACCTACTGACAAATATATCAGAACTAGCAAGCATA…TTTTTCACAGAACTAGAACAAAAATTTTAAAATTTG U6

[17631] 3000 TACGAATGAACTCGGAAGTATTTTGTGGTATCCTTT…AGGTAGCACACACCCAGGCGATCTTGCAAGTAAGAG ssc-mir-9785

[17632] 3000 GTTTCGGGTTTTCTCCACCCCCTCTCCCTCCTGTTT…TCATTCTGTCGCGGGCCCCTCGCTGCTCCTCCCCCA ssc-mir-6782

[17633] 3000 ACGTCTGCCCTATCAACTTTCGATGGTAGTCGCCGT…AGTGTGCCGTGCCGCCATTCCAAACCAAGCACGTAC 5_8S_rRNA

filter seq length

fasta_filtered <- seq_list[width(seq_list) > 1]

output fasta format

writeXStringSet(fasta_filtered,filepath = “Sus_scrofa_3kpromoter.fasta”,format=”fasta”)
IGV检查负链基因ALDH1A1启动子序列前后各一部分位置：

IGV检查正链基因ANXA1启动子序列前后各一部分位置：

以上结果表明我们提取的启动子序列没有问题。
构建 cisTarget databases

❝
前面准备好了文件可以开始运行 create_cistarget_motif_databases.py 了。
python create_cisTarget_databases/create_cistarget_motif_databases.py 、
-f Sus_scrofa_3kpromoter.fasta 、
-M motif_dir/ -m motifs_id.txt 、
-t 15 -o Sus_scrofa

最终生成了 Sus_scrofa.motifs_vs_regions.scores.feather， Sus_scrofa.regions_vs_motifs.rankings.feather， Sus_scrofa.regions_vs_motifs.scores.feather 三个文件：

其它两个文件

❝
转录因子到 motif 的注释文件 如何准备官网并没有说明，github上也有人提了这个问题，作者建议使用同源基因去替换官网的示例文件里面的内容，如果没用就删除那一行。比如你可以使用 OrthiFinder 软件去找自己物种和人的同源基因，再处理官网的示例文件。感觉对于其它物种的支持确实不太好，建议作者出一个手把手教程去做这么一个分析还是很必要的。
官网motif2TF文件示例：

结尾

❝
路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。