前面測(cè)試了常見的人類和小鼠的數(shù)據(jù)集，但是我自己是做植物的。如果想要運(yùn)行pySCENIC，就需要去建立自己的數(shù)據(jù)庫(kù)。關(guān)于這塊網(wǎng)上的資料，我自己也是遇到了一些問題，現(xiàn)在把自己的流程與大家分享。我們先拿植物里面的模式生物擬南芥為例子來試下，然后再試下其它物種的。
官網(wǎng)的教程：
https://github.com/aertslab/create_cisTarget_databases
但是官網(wǎng)的敘述很簡(jiǎn)單，或者基本沒什么敘述。我們先來查看一下幫助文件：
python create_cistarget_motif_databases.py -h
usage: create_cistarget_motif_databases.py [-h] -f FASTA_FILENAME [-F ORIGINAL_SPECIES_FASTA_FILENAME]
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? -M MOTIFS_DIR -m MOTIFS_LIST_FILENAME
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [-5 MOTIF_MD5_TO_MOTIF_ID_FILENAME] -o DB_PREFIX
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [-c CLUSTER_BUSTER_PATH] [-t NBR_THREADS]
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [-p CURRENT_PART NBR_TOTAL_PARTS]
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [-g EXTRACT_GENE_ID_FROM_REGION_ID_REGEX_REPLACE]
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [-b BG_PADDING] [--min MIN_NBR_MOTIFS] [--max MAX_NBR_MOTIFS]
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [-l] [-s SEED] [-r SSH_COMMAND]
從幫助文件來看，有幾個(gè)必須的文件：
-f FASTA_FILENAME, --fasta FASTA_FILENAME FASTA filename which contains the regions/genes to score with Cluster-Busterforeach motif. When creating a cisTarget species database from regions/genes lifted over from a different species, provide the original FASTA fileforthat species to -F.? ? #這應(yīng)該是我們關(guān)心的區(qū)域或者基因的序列文件
-M MOTIFS_DIR, --motifs_dir MOTIFS_DIR Path to directory with Cluster-Buster motifs.
?-m MOTIFS_LIST_FILENAME, --motifs MOTIFS_LIST_FILENAME Filename with list of motif IDs or motif MD5 names to be scored from directory specified by"--motifs_dir".
#這2個(gè)文件應(yīng)該是對(duì)應(yīng)的。-m指的應(yīng)該是motif的ID文件。-M指的是文件夾的路徑，這個(gè)文件夾中包含每個(gè)Cluster-Buster 相關(guān)的motif。
因?yàn)槲乙郧耙膊粚ｉT研究motif相關(guān)的，對(duì)于這塊所說的Cluster-Buster格式或者相關(guān)也不太了解。所以又專門查看了一下Cluster-Buster。
========Cluster-Buster========
官網(wǎng)文件如下：
https://github.com/weng-lab/cluster-buster
官網(wǎng)給出的用法。從官網(wǎng)的用法來說，應(yīng)該是根據(jù)motif的matrix文件去鑒定給定的序列包含不包含相應(yīng)的motif以及詳細(xì)信息。
從輸入文件的要求來看，matrix文件需要的格式是>ID加上每個(gè)位置的概率矩陣信息。
所以下面我們需要獲得TF binding相關(guān)motif的matrix。
=====下載擬南芥轉(zhuǎn)錄因子相關(guān)的binding motif=====
這些轉(zhuǎn)錄因子相關(guān)的motif研究算是比較多的，其中可以去下載的網(wǎng)站有：JASPAR（這個(gè)是包含物種比較全的）等。但是研究植物，我喜歡去專門植物相關(guān)的網(wǎng)站去下載，我是從PlantTFDB（http://planttfdb.gao-lab.org/index.php）上面下載的。
從這里下載的TF binding motif相關(guān)的信息。
#Ath_TF_binding_motifs_information.txt? ?文件里面包含所有motif的基本信息
Ath_TF_binding_motifs.meme? 可以明顯看出，文件里面包含每個(gè)motif相關(guān)的matrix信息。（注：根據(jù)cbust的要求，我們需要處理一下這個(gè)文件，只保留ID和矩陣信息。）
下面就是我處理之后的文件Ath.TF.matrix.txt，只保留了ID和矩陣信息。（為了方面后續(xù)查看，我把名字也加上去了，放到了ID后面）
下面我們就可以去執(zhí)行相關(guān)的文件了，比如我獲得了擬南芥所有基因的promoter 3K的序列，就可以去查看那些基因promoter包含相應(yīng)的motif了。
cbust? -f 1 Ath.TF.matrix.txt Ath.promoter3K.V2.fa
輸出結(jié)果如下圖所示：
剛開始，我理解的create_cistarget_motif_databases.py中 -M文件夾中的輸入文件是cbust的的輸出結(jié)果。所以我對(duì)于每個(gè)motif都運(yùn)行了一下cbust，生成結(jié)果放在了motif_cb_format文件夾中，獲得了默認(rèn)格式的每個(gè)motif的cb文件。
========運(yùn)行create_cistarget_motif_databases.py====
python create_cistarget_motif_databases.py -f Ath.promoter3K.V2.fa -M motif_cb_format/ -m Ath.TF.txt -o Athaliana -t 30
但是報(bào)錯(cuò)了。
看報(bào)錯(cuò)文件，才發(fā)現(xiàn)原來create_cistarget_motif_databases.py程序中自己要去運(yùn)行cbust，不是基于我運(yùn)行的結(jié)果去解讀。
所以就把motif_cb_format文件夾中的文件換成每個(gè)motif的矩陣文件了。
所以重新運(yùn)行剛才那個(gè)命令，就成功運(yùn)行了：
運(yùn)行輸入文件格式和詳情：
Ath.promoter3K.V2.fa? #就是每個(gè)基因promoter 3K的序列
Ath.TF.txt? ?#TF相關(guān)motif的ID，這個(gè)與文件中motif_cb_format名字一一對(duì)應(yīng)。
motif_cb_format? 文件夾，與Ath.TF.txt 中ID一一對(duì)應(yīng)
最終就產(chǎn)生了我們想要的擬南芥相關(guān)的文件。
