seqkit根据基因id_基因家族分析保姆级教程（分子进化）-生信小白自学之路

开始前先熟悉一下我自己

简单介绍一下，我是20级研究生，普通本科和中科院海洋所联培的学生，什么是联培呢，简单来说就是理论课程在学校上，根据学校安排，一般是一年，我们半年就上完了，后面就去海洋所工作了。正式的说就是工作，不是上学，就是打工，个人感觉，不过这并没有什么影响，丝毫不会耽误一个人进步。联培算是三赢的方式，作为学生来说要好好珍惜这个平台。

下面介绍我对基因家族的认知过程

我是小白，在这之前我什么都不知道，我很迷茫，一头雾水。

然后和海洋所的导师见了一次面，她给我布置了这个课题，导师说很简单，师兄说这就是个小课题。我？哦？那好吧，当我实践起来，？？？emmmmmmmm？？？？？？？？

真的是，不知所措。第一次开组会是在定了课题后的一周，我做了ppt就讲啊讲啊，他们说没什么内容....我也确实是花了时间展现一点东西，可是实在是没有实质性的进展。后来我仔细一看，确实没什么东西，我就发现，这和学校不一样，你得拿出实质性的进展、结果才可以。

师兄帮我安装了相关软件，虽然我不知道他在干啥，现在我觉得我也行的操作，在前段时间的我看来，哇好牛逼，大神大神。感激。真的，我感激每一个帮助过我的人。然后我的师姐大体给我讲了一下思路，最初的思路。外文网站，咱都没接触过，我就照猫画虎的模仿，尝试，现在终于算是搞明白了那些套路，明白了原理很多就好办多了。我觉得做这个最重要的是思路，思路怎么来呢，简书这个平台不错，我在这上面学了很多操作，还有就是文献，但是文献的方法讲的基本是一笔带过，讲个专业名词。咱也不懂。师姐讲的也不是那么仔细，我也不能老麻烦人家，所以基本都是我自己不断的摸索，现在我仿佛到了那种境界，豁然开朗的境界，这个东西不要和别人比，优秀的人实在是太多了，和自己比就好了，自己进步了自己开心就好了，非要和别人比，那对于别人来说可能就是小case啦。所以自己好好搞这一方面就好了。

开始正文

一、基因家族基础知识与研究思路介绍

1.基因家族简介

基因是染色体上一段可以发生转录的区域(内含子外显子启动子)。转录本才是基因的研究实体，转录本transcript(别名 剪切体)是由一条基因通过转录形成的一种或多种可供编码蛋白质的成熟的mRNA。一条基因通过内含子的不同剪接可构成不同的转录本。设计转录本实验可以研究内含子剪切机制、表观遗传、RNA编辑等，通常是考察一条基因对应的不同转录本的调节机制等。

基因家族(gene family)，是来源于同一个祖先，由一个基因通过基因重复而产生两个或更多的拷贝而构成的一组基因，它们在结构和功能上具有明显的相似性，编码相似的蛋白质产物，同一家族基因可以紧密排列在一起，形成一个基因簇，但多数时候，它们是分散在同一染色体的不同位置，或者存在于不同的染色体上的，各自具有不同的表达调控模式。

按功能划分：把一些功能类似的基因聚类，形成一个家族，例如GH家族(糖苷水解酶家族)等。按照序列相似程度划分：一般将同源的基因放在一起认为是一个家族，一般使用orthoMCL进行聚类

motif是蛋白质分子具有特定功能的或者作为一个独立结构域一部分相近的二级结构聚合体序列高度相似的序列，互为同源gene，归属于一个基因家族(拷贝数目多于1)

结构域的角度来说，具有保守结构域(某个或多个)的序列，即为某个基因家族的序列(可能同时要不具有另外的某个结构域)

2.常规的基因家族分析流程

在这些常规的生信分析后，一般的文章还会加上一些湿实验去验证，例如不同非生物条件下基因家族的表达等(PCR为主)。1.确定的研究基因家族

2.了解你研究的基因家族的特征

3.可参考收录了基因家族特征的网站

4.查找相关文献

5.数据下载

A.基因组序列信息，存储基因组序列信息的.fasta文件。还有其蛋白质序列，也是以.fasta结的文件。一般来说注释的比较好的基因组都会含有这些文件。

B.基因组基因结构注释信息。储存基因的intron，exon，CDS,gene等坐标信息的.gff3或.gtf文件。

C.基因家族隐马可夫模型，hmm文件

3.基因家族鉴定的工具hmmer：

一般寻找基因家族，都可以通过保守结构域来预测，从而找到物种的某一基因家族

在鉴定基因家族时，常用到的工具是hmmsearch，里面常用的算法有三种。一般我们使用--cut_tc算法对隐马可夫模型进行搜索，tc算法是使用pfam提供的hmm文件中trusted cutoof的值进行筛选，相对比较可靠

二、基因家族分析|基因家族成员鉴定(hmm模型＆同源blast)

1 基因家族成员的鉴定步骤详解确定研究的基因家族

家族成员的基本特征确定(参考已有物种)

参考序列集合的准备

目标物种序列和注释信息的下载或准备

双向Blast比对获取可能的成员

基于保守结构域进行进一步筛选

双向Blast比对获取可能的成员

方法一：基于hmm模型的鉴定方法

准备数据下载研究物种基因组fasta文件、注释文件gtf/gff3文件

BIR.hmm    #PF00653.22       这里是举例

目标基因家族搜索与筛选hmmsearch --cut_tc --domtblout 123.out BIR.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa.gz

#过滤筛选得到E-value小于1*10-20,先拿到序列号

grep -v "#" BIR.out|awk '($7 + 0) < 1E-20'|cut -f1 -d  " "|sort -u > BIR_qua_id.txt

#再根据序列号，从Arabidopsis_thaliana.TAIR10.pep.all.fa.gz中提取序列

less Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | /data1/spider/ytbiosoft/seqkit grep -f BIR_qua_id.txt > BIR_qua.fa

多序列比对，构建目标物种的NB-ARC基因家族的hmm模型#对筛选出来的序列用clustalw进行多序列比对

/data/shaofeng/clustalw/clustalw

弹出clustalw的操作界面，以下展示具体输入过程：

选择1(输入待比对序列)→ 输入待比对序列的文件名：BIR_qua.fa → 选择2(开始进行序列比对)→选择9(选择输出比对结构的格式为aligned)→ 按enter键 → 选择1(选择比对模式为全局比对)→ 指定输出的比对结果的文件名称:BIR_qua.aln → 回车后开始比对 → 输入一个树文件名(new GUIDE TREE file):BIR_qua.dnd (最后才能得到BIR.aln，否则BIR.aln为空)

#使用hmmbuild对这些置信的序列进行隐马尔可夫模型的构建，即构建更加准确的hmm模型来尽可能的预测目标物种中BIR基因家族中所有的成员。

hmmbuild BIR_qua.hmm  BIR_qua.aln

hmmsearch --cut_tc --domtblout BIR.second.out BIR_qua.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa

利用目标物种的hmm模型再次筛选目标物种中符合要求的序列#再次对这些基因进行过滤和提取

grep -v "#" NBS-ARC.second.out|awk '($7 + 0) < 1E-03' | cut -f1 -d " "|sort -u >final.NBS.list

less Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | /data1/spider/ytbiosoft/seqkit grep -f final.NBS.list > final_NBS-ARC_qua.fa

方法二：基于同源比对blast的鉴定方法

下载NCBI 中所有动物存在于Ref-seq中的IAP序列(Ref-seq一般被认为是比较置信的动物基因序列)

将下载的蛋白序列存放至ref.nbs.plant.fa文本文档中，上传至服务器

比对并筛选目标物种中符合要求的序列#用makeblastdb建立blast数据库

makeblastdb -in ref.nbs.plant.fa -dbtype prot -out blastdb

#用blastp进行序列搜索，得到每个序列的相似序列

blastp -num_threads 20 -db blastdb -query Arabidopsis_thaliana.TAIR10.pep.all.fa -outfmt 7 -seg yes > blastp.out &

#筛选identity大于75%的序列

cat blastp.out |awk '$3>75' |cut -f1 |sort -u > blastp_result_id.list

将上述两种方法得到gene id合并取交集，找出两种方法共有的基因家族成员，使结果更可信.comm -12 blastp_result_id.list final.NBS.list > common.list

less Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | /data1/spider/ytbiosoft/seqkit grep -f common.list > final_searh_NBS-ARC_qua.fa

最后，还可以通过一些网上的保守结构域搜索网页，进一步对所找出的结果进行验证，比如：这些工具都可再次验证所搜寻的蛋白质序列是不是含有基因家族对应domain。在查看保守结构域后，如果该区域含有IAP所对应的保守结构域，例如BIR区域等，该蛋白质序列可以保留进行后续的分析。如果在该区域没有找到对应的保守区域，为了分析的严谨性，需进行进一步的排查来确定是否要去掉该序列。这种情况一般分为两种情况，第一种就是注释无误，该序列确实丢失了对应的保守结构域，需要去掉。第二种情况就是注释有误，该序列的结构域可能没有被完整的保留下来，这种情况应该截取该序列的上下游重新注释分析。

总结及注意事项

只有一个domain，hmmer很快，但是可能结果很多，例如MAPK、MAPKK、MAPKKK等，它们的domain都为pkinase，分族是根据进化树分支结果，这时需要结合blast结果验证。

两个及以上domain，需要利用检索取两个较少结果的交集，可结合blast结果验证。