前面测试了常见的人类和小鼠的数据集,但是我自己是做植物的。如果想要运行pySCENIC,就需要去建立自己的数据库。关于这块网上的资料,我自己也是遇到了一些问题,现在把自己的流程与大家分享。我们先拿植物里面的模式生物拟南芥为例子来试下,然后再试下其它物种的。
官网的教程:
https://github.com/aertslab/create_cisTarget_databases
但是官网的叙述很简单,或者基本没什么叙述。我们先来查看一下帮助文件:
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]
从帮助文件来看,有几个必须的文件:
-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. #这应该是我们关心的区域或者基因的序列文件
-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个文件应该是对应的。-m指的应该是motif的ID文件。-M指的是文件夹的路径,这个文件夹中包含每个Cluster-Buster 相关的motif。
因为我以前也不专门研究motif相关的,对于这块所说的Cluster-Buster格式或者相关也不太了解。所以又专门查看了一下Cluster-Buster。
========Cluster-Buster========
官网文件如下:
https://github.com/weng-lab/cluster-buster
官网给出的用法。从官网的用法来说,应该是根据motif的matrix文件去鉴定给定的序列包含不包含相应的motif以及详细信息。
从输入文件的要求来看,matrix文件需要的格式是>ID加上每个位置的概率矩阵信息。
所以下面我们需要获得TF binding相关motif的matrix。
=====下载拟南芥转录因子相关的binding motif=====
这些转录因子相关的motif研究算是比较多的,其中可以去下载的网站有:JASPAR(这个是包含物种比较全的)等。但是研究植物,我喜欢去专门植物相关的网站去下载,我是从PlantTFDB(http://planttfdb.gao-lab.org/index.php)上面下载的。
从这里下载的TF binding motif相关的信息。
#Ath_TF_binding_motifs_information.txt 文件里面包含所有motif的基本信息
Ath_TF_binding_motifs.meme 可以明显看出,文件里面包含每个motif相关的matrix信息。(注:根据cbust的要求,我们需要处理一下这个文件,只保留ID和矩阵信息。)
下面就是我处理之后的文件Ath.TF.matrix.txt,只保留了ID和矩阵信息。(为了方面后续查看,我把名字也加上去了,放到了ID后面)
下面我们就可以去执行相关的文件了,比如我获得了拟南芥所有基因的promoter 3K的序列,就可以去查看那些基因promoter包含相应的motif了。
cbust -f 1 Ath.TF.matrix.txt Ath.promoter3K.V2.fa
输出结果如下图所示:
刚开始,我理解的create_cistarget_motif_databases.py中 -M文件夹中的输入文件是cbust的的输出结果。所以我对于每个motif都运行了一下cbust,生成结果放在了motif_cb_format文件夹中,获得了默认格式的每个motif的cb文件。
========运行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
但是报错了。
看报错文件,才发现原来create_cistarget_motif_databases.py程序中自己要去运行cbust,不是基于我运行的结果去解读。
所以就把motif_cb_format文件夹中的文件换成每个motif的矩阵文件了。
所以重新运行刚才那个命令,就成功运行了:
运行输入文件格式和详情:
Ath.promoter3K.V2.fa #就是每个基因promoter 3K的序列
Ath.TF.txt #TF相关motif的ID,这个与文件中motif_cb_format名字一一对应。
motif_cb_format 文件夹,与Ath.TF.txt 中ID一一对应
最终就产生了我们想要的拟南芥相关的文件。