在上一篇文章中，教大家怎么去鉴定一个基因家族，但是不是所有的基因家族都是只有一个结构域的，有些基因家族就包含很多个结构域，例如：GT family、NBS-LRR family、2OGD family等等。他们包含三个以上的结构域。那我们还在用之前的方法，一个一个去hmmsearch，那不得累死。
那就要学会用工具，用循环的方式去解决问题，话不多说，直接开始搞~
##假如你有很多的结构域模型（.hmm文件），都放在一个文件夹中取名为hmm_model
#进入hmm_model文件夹并打印出hmm_model所在的位置
cd hmm_model && pwd
/mnt/d/tmp/hmm_model
得到/mnt/d/tmp/hmm_model文件夹的信息后，我们就要学习一下如何去编写简单的脚本
#!/usr/bin/bash
#编辑变量
data_dir=/mnt/d/tmp/hmm_model
#蛋白序列所在的位置
pep=
#cds列所在的位置
cds=
#生成文件的前缀（一般为基因家族的名称）
prefix=
#设置参考物种的相应的家族基因蛋白序列位置
Ath=
#生成结构域文件所在路径的文件
ls $data_dir/*.hmm > hmm.txt
#开始循环进行hmmsearch
for i in hmm.txt
do
echo $i
n=$(echo ${i} | sed 's/\.hmm//;s#\/mnt\/d\/tmp\/hmm_model\/##g')#引入一个新的变量
hmmsearch --cut_tc --domtblout $n.domtblout -o $n.hmmout $i $pep
awk '$7<0.01 && $1 !~ /^#/ {print $0 }' $n.domtblout >$n.domtblout.filter
awk '{print $1}' $n.domtblout.filter | sort -u > $n.domtblout.filter.id
done
#blast比对的话还是一样的，毕竟参考的模式物种选一个或者两个就行
blastp -query $pep -db $Ath -evalue 1e-10 -outfmt '6 std qlen slen ' -out $prefix.blastout
#blast结果基于identity30%过滤
awk '$3 > 30 && $4> 155 {print $1}' $prefix.blastout | sort -u > $prefix.blastout.filter.id
###具体情况具体分析，二选一即可
#合并hmm和blast结果,也就是取并集（避免缺漏）
cat *.domtblout.filter.id $prefix.blastout.filter.id | sort -u > ${prefix}_final.geneID
##当然也可以取交集（更加精确）
cat *.domtblout.filter.id $prefix.blastout.filter.id | sort | uniq -C | awk '$1 == 2 {print $2}' > ${prefix}_final.geneID
#提取目标基因的pep和cds序列
seqtk subseq $pep ${prefix}_final.geneID >${prefix}_final.pep.fasta
seqtk subseq $cds ${prefix}_final.geneID > ${prefix}_final.cds.fasta
这样的话，就可以省去很多重复的命令，并且由nohup提交到后台，省事！！！要做其他物种或者其他多结构的家族基因，换一换变量就可以用了~懒人要有懒办法
网友评论