已知一个蛋白domain 找到对应的基因家族并做进化树

1)根据domain序列，直接打开pfam网站，首页上面有一个sequence search 输入序列后得到pfam号

2)之后在自己服务器上操作：找到目标物种的注释文件（一般会带有pfam号） 输入命令 grep pf*** <注释文件> |less 先查看，确定后直接输出到新文件  grep pf*** <注释文件>  |cut -f 1 > protain_transcript 

3)之后可使用软件seqkit 下载后解压后可以直接使用 需要有基因组protain的序列(fasta文件)

seqkit grep -f protain_transcript  <genome_protain.fasta>   > goal_protain_family.fa

得到该文件后使用windows版本的mega，可参考下面的比对步骤（6）

以下方法较麻烦，建议使用上述方法

前提有一个domain序列（需要做成fasta格式，直接用vim 编辑器进行编辑文档即可 。第一行为>name，第二行为序列），基因组的蛋白序列（fasta文件）

1.需要通过blastp来得到相应的基因（转录本） 。以下称呼 query fasta （已知的蛋白序列） 蛋白库（基因组的蛋白序列通过blast中的makeblastdb建索引）

1）首先对基因组蛋白序列进行建库    

-in 用于输入基因蛋白序列位置    -dbtype   建库类别为蛋白则输入prot  后面两个参数可以直接加上 -out 为输出文件位置

2）接下来进行比对

blastp -query <比对的文件> -db <建立的蛋白库索引> -evalue 1e-10  -out <输出的文件名> -outfmt <有多个输出格式，一般常用6，此处用6>  -num_alignments <用来显示前n个比分质量高的，此处选10>  -num_threads <使用的线程数，看自己文件大小及服务器承载能力而定>

3）result（12列）

重点关注前几个，此处我卡的是identity>70 ，query 长度为163，我卡的是length>105

4)此处文件应该是这样的

queryid  基因组蛋白序列中对应的id    100.000    163    0    0    1    163   17    179    4.40e-123    346

5）因为我得到的是对应的基因ID（转录本id）需要得到对应的fasta文件，此处用到软件seqkit （谷歌该软件，复制下载链接，直接在linux上使用 命令       wget   <加上复制的下载链接>   下载完成后使用   tar  -zxvf  <下载的文件名>  解压后即可使用，我此处用的是绝对路径）

使用命令     ~/biosoft/seqkit grep -f <将blast后的基因id该列拿出来> <基因组的fasta文件>    > prootain.fa

6）得到fa文件后，打开在自己电脑上下载好的mega软件（windows版的）载入文件后点击窗口左上方的aligment按钮，提示未选中文件，是否选择全部文件，点击是后，选择第一个 clustalW ，弹出窗口，点击确定后进行比对