使用AnimalTFDB3.0对动物基因组的转录因子进行注释

1. 动物转录因子数据AnimalTFDB3.0简介

动物转录因子数据库AnimalTFDB3.0对97个动物基因组的转录因子(Transcription Factor)和转录辅助因子(Transcription cofactor)进行了归纳整理。基于DNA结合结构域,将动物转录因子分成了73个基因家族,将转录辅助因子分成了83个基因家族。此外,动物转录因子分为六大类(Basic Domain Group、Zinc-Coordinating Group、Beta-Scaffold Factors、Helix-turn-helix、Other Alpha-Helix Group和Unclassified Structure),动物转录辅助因子也分为六大类(Co-activator/repressors、Chromatin Remodeling Factors、General Cofactors、Histone-modifying Enzymes、Cell Cycle和Other Cofactors)。

动物转录因子数据库AnimalTFDB3.0提供了网页工具进行转录因子分析。该网页工具一次仅支持上传不超过1000条蛋白序列。该网页工具采用HMM算法进行计算,分析速度比较快。数据库作者对大部分转录因子家族采用了PFam数据库的HMM模型,并自己构建了一些转录因子家族的HMM模型。但作者没有提供HMM数据库的下载,仅提供了97个动物转录因子和转录辅助因子蛋白序列的FASTA格式文件下载。

使用AnimalTFDB3.0数据库的网页工具仅能进行转录因子分析,且由于仅能支持不超过1000条序列进行分析而比较麻烦,不利于动物全基因组的转录因子分析。以下讲解下载AnimalTFDB3.0数据库FASTA文件,并自行编写程序进行转录因子和转录辅助因子注释。

2. 下载AnimalTFDB3.0数据库蛋白序列

首先,进入AnimalTFDB3.0数据库下载页面,采用复制粘贴方式获取97个物种的物种名,保存到文件species.txt中。

Ailuropoda melanoleuca;Anas platyrhynchos;Anolis carolinensis;Aotus nancymaae;Astyanax mexicanus;Bos taurus;Caenorhabditis elegans;Callithrix jacchus;Canis familiaris;Capra hircus;Carlito syrichta;Cavia aperea;Cavia porcellus;Cebus capucinus;Cercocebus atys;Chinchilla lanigera;Chlorocebus sabaeus;Choloepus hoffmanni;Ciona intestinalis;Ciona savignyi;Colobus angolensis palliatus;Cricetulus griseus chok1gshd;Cricetulus griseus crigri;Danio rerio;Dasypus novemcinctus;Dipodomys ordii;Drosophila melanogaster;Echinops telfairi;Equus caballus;Erinaceus europaeus;Felis catus;Ficedula albicollis;Fukomys damarensis;Gadus morhua;Gallus gallus;Gasterosteus aculeatus;Gorilla gorilla;Heterocephalus glaber female;Heterocephalus glaber male;Homo sapiens;Ictidomys tridecemlineatus;Jaculus jaculus;Latimeria chalumnae;Lepisosteus oculatus;Loxodonta africana;Macaca fascicularis;Macaca mulatta;Macaca nemestrina;Mandrillus leucophaeus;Meleagris gallopavo;Mesocricetus auratus;Microcebus murinus;Microtus ochrogaster;Monodelphis domestica;Mus caroli;Mus musculus;Mus pahari;Mus spretus;Mustela putorius furo;Myotis lucifugus;Nannospalax galili;Nomascus leucogenys;Notamacropus eugenii;Ochotona princeps;Octodon degus;Oreochromis niloticus;Ornithorhynchus anatinus;Oryctolagus cuniculus;Oryzias latipes;Otolemur garnettii;Ovis aries;Pan paniscus;Pan troglodytes;Papio anubis;Pelodiscus sinensis;Peromyscus maniculatus bairdii;Petromyzon marinus;Poecilia formosa;Pongo abelii;Procavia capensis;Propithecus coquereli;Pteropus vampyrus;Rattus norvegicus;Rhinopithecus bieti;Rhinopithecus roxellana;Saimiri boliviensis boliviensis;Sarcophilus harrisii;Sorex araneus;Sus scrofa;Taeniopygia guttata;Takifugu rubripes;Tetraodon nigroviridis;Tupaia belangeri;Tursiops truncatus;Vicugna pacos;Xenopus tropicalis;Xiphophorus maculatus

然后,根据物种名,编写下载数据库文件的脚本(每个物种下载4个数据文件)并批量下载数据。

$ mkdir /opt/biosoft/AnimalTFDB3.0
$ cd /opt/biosoft/AnimalTFDB3.0

生成文件species.list,包含97个物种名,每行一个物种名。根据http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/download网>页中的信息,使用vim编辑和perl单行命令抓取全部的物种名信息。
$ perl -p -i -e 's/;/\n/g;' species.txt

$ mkdir data_of_97_species
$ cd data_of_97_species
$ perl -e 'while (<>) { chomp; s/\s+/_/g; print "wget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF_protein.fa\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_Cof_protein.fa\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF_cofactors\n"; }' ../species.list > download.list
$ ParaFly -c download.list -CPU 2

Caenorhabditis_elegans的FASTA文件有问题,存在序列名一致的多条序列,推荐进行修改。
$ perl -p -i -e 's/>(\S+)\t(\S+)/>$2\t$1/' data_of_97_species/Caenorhabditis_elegans_*_protein.fa
$ cd ..

分别合并转录因子数据库和转录辅助因子数据库的FASTA文件。

得到转录因子序列数据库,FASTA文件头部仅包含序列名和基因家族信息。
$ cat data_of_97_species/*TF_protein.fa > AnimalTFDB3.0_TF_protein.fasta
$ perl -p -i -e 's/^(>\S+)\t.*?\t.*?\t/$1\t/' AnimalTFDB3.0_TF_protein.fasta


得到转录辅助因子数据库,FASTA文件头部仅包含序列名和基因家族信息。
$ cat data_of_97_species/*_Cof_protein.fa > AnimalTFDB3.0_Cof_protein.fasta
由于Cof_protein.fa文件中不包含转录辅助因子基因家族信息,需要编写程序读取TF_cofactors文件信息并将基因家族信息写入FASTA文件
$ perl -e 'while (<>) { @_ = split /\t/; $gf{$_[2]} = $_[3]; } open IN, "AnimalTFDB3.0_Cof_protein.fasta"; while (<IN>) { if (m/^>(\S+)\t(\S+)/) { if (exists $gf{$1}) {print ">$1\t$gf{$1}\n";} elsif (exists $gf{$2}) {print ">$2\t$gf{$2}\n";} else {print STDERR "No type: $_\n";} } else { print; } }' data_of_97_species/*_TF_cofactors > aa; mv aa AnimalTFDB3.0_Cof_protein.fasta

3. 使用DIAMOND将数据库所有蛋白和数据库自身进行比对,计算每个基因家族的BLAST阈值

使用ALL_VS_ALL的方式分别对转录因子数据库和转录辅助因子数据库进行DIAMOND比对(evalue 1e-5、identity 30%且coverage 30%)。根据DIAMOND比对结果,得到每个转录因子家族中每个基因的比对信息(evalue、hitNum),编写程序进而计算每个基因家族的BLAST阈值。

$ diamond --in AnimalTFDB3.0_TF_protein.fasta --db AnimalTFDB3.0_TF_protein
$ diamond blastp --db AnimalTFDB3.0_TF_protein --query AnimalTFDB3.0_TF_protein.fasta --out TF.tab --outfmt 6 --sensitive --max-target-seqs 50000 --evalue 1e-5 --id 30 --subject-cover 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1

编写程序geneFamily_BLAST_threshold_of_sensitivity.pl,输入DIAMOND结果文件和数据库FASTA文件,得到各个sensitivity级别上每个基因家族的evalue和hitNum阈值。

4. 将全基因组蛋白序列和数据库进行比对,再根据阈值文件进行准确注释。

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据