使用HMMER从转录组数据中寻找基因家族

1.软件准备

1.1HMMER

conda create -n HMMER
conda activate HMMER
conda install -c bioconda hmmer
#创建环境并安装

1.2 Trinity

conda create -n trinity
conda activate trinity
conda install -c bioconda trinity=2.15.1
#创建环境并安装指定版本的trinity

1.3 transdecoder

conda create -n transdecoder
conda activate transdecoder
conda install -c bioconda transdecoder=5.7.1
#创建环境并安装指定版本的transdecoder

2.蛋白序列获得

2.1 转录本拼接

Trinity --seqType fq --max_memory 50G --left reads_1.fq  --right reads_2.fq --CPU 6

--seqType fq指定测序文件的类型
--max_memory指定最大使用内存
--left/right双端测序文件名称
--CPU 最大CPU使用数目

2.2 提取最长转录本

get_longest_isoform_seq_per_trinity_gene.pl Trinity.fasta > longest.fasta

Trinity.fasta输入需要查找最长转录本的转录本文件，longest.fasta输出文件名称。

2.3转录本翻译

TransDecoder.LongOrfs -t longest.fasta --output_dir /path/to/your/PEP

longest.fasta上一步获得的最长转录本，/path/to/your/PEP输出文件夹路径。

3.查找基因

3.1下载结构域

wget https://www.ebi.ac.uk/interpro/wwwapi//entry/pfam/PF02441?annotation=hmm
gunzip PF02441.hmm.gz

从Pfam网站（Browse - InterPro (ebi.ac.uk)）上下载需要的结构域并解压缩。

3.2构建索引与运行

hmmpress PF02441.hmm
hmmscan --tblout output.txt PF02441.hmm protein.fasta 

output.txt是输出文件的名称，protein.fasta是蛋白序列文件。
output.txt的基本内容如下：

#                                                                      --- full sequence ---- --- best 1 domain ---- --- domain number estimation ----
# target name        accession  query name                  accession    E-value  score  bias   E-value  score  bias   exp reg clu  ov env dom rep inc description of target
#------------------- ----------        -------------------- ---------- --------- ------ ----- --------- ------ -----   --- --- --- --- --- --- --- --- ---------------------
Flavoprotein         PF02441.23 TRINITY_DN25061_c0_g5_i2.p1 -            2.9e-34  104.2   0.0   4.3e-34  103.6   0.0   1.3   1   1   0   1   1   1   1 Flavoprotein
#
# Program:         hmmscan
# Version:         3.4 (Aug 2023)
# Pipeline mode:   SCAN
# Query file:      protein.fasta
# Target file:     PF02441.hmm
# Option settings: hmmscan --tblout CHS.txt PF02441.hmm protein.fasta 
# Current dir:     /home/oldhome/monkeyflower/xiaodeng/bin/HMMERDB/HTHC
# Date:            Tue Jan  2 09:22:11 2024
# [ok]