基因家族分析一(NBS基因家族鉴定)

https://www.bilibili.com/video/av85069377?p=2

大致流程

一：准备文件

1. 拟南芥的基因组，CDS，蛋白，GFF以及HMM文件。

GCF_000001735.4_TAIR10.1_genomic.fna
GCF_000001735.4_TAIR10.1_genomic.gff
GCF_000001735.4_TAIR10.1_protein.faa
GCF_000001735.4_TAIR10.1_rna.fna
NB-ARC.hmm

2.利用HMM文件搜寻pep文件

#利用hmmsearch搜寻pep
$hmmsearch --cut_tc --domtblout NBS-ABC.out NBS-ARC.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa
#过滤筛选得到E-value小于1*10-20,先拿到序列号
$grep -v "#" NBS-ABC.out|awk '($7 + 0) < 1E-20'|cut -f1 -d  " "|sort -u > NBS-ARC_qua_id.txt
#再根据筛选的序列号，从拟南芥pep文件中提取序列
$less Arabidopsis_thaliana.TAIR10.pep.all.fa | /home/spider/project/yuantao/soft/seqkit grep -f NBS-ARC_qua_id.txt > NBS-ARC_qua.fa

hmmsearch结果文件

3.利用clustalw进行多序列比对，构建物种特定的NB-ARC基因家族hmm模型

流程如下：

clustalw -> 1 -> NB-ARC_qua.aln ->2 ->1--> NBS-ARC_qua.aln ->NB-ARC_qua.dnd ->X

最终生成NBS-ARC_qua.aln和NBS-ARC_qua.aln文件，其中NB.aln文件就是用于构建新的马尔科夫模型文件。

$hmmbuild new-NB-ARC.hmm NBS-ARC_qua.aln

在利用新生成的拟南芥特异的马尔科夫模型寻找pep文件中的NB-ARC基因家族序列。

$hmmsearch --cut_tc --domtblout NBS-ARC.second.out new-NB-ARC.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa

对结果文件筛选，E-value值小于0.01，ID取重复重复，最后将ID存为final_E02.NB_ids.list。

grep -v "#" NBS-ARC.second.out|awk '($7 + 0) < 1E-20' | cut -f1 -d " "|sort -u >final.NBS.list
$less Arabidopsis_thaliana.TAIR10.pep.all.fa | /home/spider/project/yuantao/soft/seqkit grep -f final.NBS.list > final_NBS-ARC_qua.fa

这个序列就可以用于后续的进化树，保守结构域，motif等分析。

二：通过blast鉴定基因家族成员（适用于没有该基因家族HMM模型）

利用CDD，Pfam，interProscan确定基因家族成员保守结构域。
利用Plant CARE和PLACE预测成员顺式作用原件。
利用别的物种的基因家族的蛋白序列，实例用玉米和水稻的NBS序列寻找拟南芥的NBS序列。

1.从NCBI下载NBS序列（尽可能下载多个其他物种的NBS序列）

取NCBI得protein条目下搜索NBS序列。

2.blastp比对并筛选目标物种中符合要求的序列。

用query=ref_nbs.plant.fa去拟南芥蛋白序列中搜索

$makeblastdb -in Arabidopsis_thaliana.TAIR10.pep.all.fa -dbtype prot -out TAIR
$blastp -num_threads 16 -db TAIR -query ref_nbs.plant.fa -evalue 1e-20 -outfmt 7 -seg yes > reftoTAIR_blastp.out
#删除#开头的行
$sed -i '/^#/d' reftoTAIR_blastp.out
$awk '{print$2}' reftoTAIR_blastp.out >NBS_second.ids
##最后利用sort和uniq函数对NBS_second.ids去重复，再将去重复的HMM和blast找到的id取交集。
$cat NBS_second.ids|sort|uniq >uniq_NBS_second.ids
$comm -12 uniq_NBS_second.ids final.NBS.list > common.list
#最后根据交集中的ID取拟南芥的pep文件中搜索可信的NB-ARC序列
$less Arabidopsis_thaliana.TAIR10.pep.all.fa | /home/spider/project/yuantao/soft/seqkit grep -f common.list > final_all_NBS-ARC_qua.fa

拿到这些蛋白序列之后，就需要到NCBI，Pfam等网站上搜索是否含有NBS保守的结构域。