PHI注释步骤（含批量）

https://www.sciencedirect.com/science/article/pii/S1878614620300015?via%3Dihub

BLASTP was run using NCBI-BLAST-2.7.1+ (Camacho et al., 2009) against fungal sequences in the UniProt database, against the PHI database version 45 (Urban et al., 2017), and against the MEROPS release 12.0 database (Rawlings et al., 2018).

17.PHI病原与宿主互作数据库

将PHI数据库下载下来本地blast，需要先填资料，填完之后即可下载至本地，使用blast进行比对。

PHI-base :: Pathogen Host Interactions

# 下载好的PHI数据库位置如下
/media/aa/DATA/SZQ2/PHI/PHI.fasta
# 在(pfam_scan)下
conda activate pfam_scan
# 新建文件夹
mkdir 17.PHI && cd 17.PHI

根据下载好的PHI数据库，将需要查找的蛋白序列，与数据库进行BLASTP比对，命令如下：

1） diamond blastp

 evalue为1e-10

diamond blastp --db /media/aa/DATA/SZQ2/PHI/PHI.fasta --query /media/aa/DATA/SZQ2/bj/functional_annotation/pep70/$i.fasta --out $i.phi10.xml --outfmt 5 --sensitive --max-target-seqs 20 --evalue 1e-10 --id 20 --tmpdir /dev/shm --index-chunks 1

批量操作

批量 pep70

for i in `cat /media/aa/DATA/SZQ2/bj/functional_annotation/70list.txt`
do
    echo "diamond blastp --db /media/aa/DATA/SZQ2/PHI/PHI.fasta --query /media/aa/DATA/SZQ2/bj/functional_annotation/pep70/$i.fasta --out $i.phi10.xml --outfmt 5 --sensitive --max-target-seqs 20 --evalue 1e-10 --id 20 --tmpdir /dev/shm --index-chunks 1"
done > command.phi.list
ParaFly -c command.phi.list -CPU 48

批量 pepmy

for i in `cat /media/aa/DATA/SZQ2/bj/functional_annotation/pepmylist.txt`
do
    echo "diamond blastp --db /media/aa/DATA/SZQ2/PHI/PHI.fasta --query /media/aa/DATA/SZQ2/bj/functional_annotation/pepmy/$i.fasta --out $i.phi10.xml --outfmt 5 --sensitive --max-target-seqs 20 --evalue 1e-10 --id 20 --tmpdir /dev/shm --index-chunks 1"
done > command.phi.list
ParaFly -c command.phi.list -CPU 48

2）parsing_blast_result.pl

新建文件夹

mkdir phi10.tab && cd phi10.tab

/media/aa/DATA2/bin/parsing_blast_result.pl --evalue 1e-10 --HSP-num 1 --out-hit-confidence --suject-annotation ../$i.phi10.xml > $i.phi10.tab

批量操作

批量 pep70

for i in `cat /media/aa/DATA/SZQ2/bj/functional_annotation/70list.txt`
do
    echo "/media/aa/DATA2/bin/parsing_blast_result.pl --evalue 1e-10 --HSP-num 1 --out-hit-confidence --suject-annotation ../$i.phi10.xml > $i.phi10.tab"
done > command.phi.list
ParaFly -c command.phi.list -CPU 48

 批量 pepmy

for i in `cat /media/aa/DATA/SZQ2/bj/functional_annotation/pepmylist.txt`
do
    echo "/media/aa/DATA2/bin/parsing_blast_result.pl --evalue 1e-10 --HSP-num 1 --out-hit-confidence --suject-annotation ../$i.phi10.xml > $i.phi10.tab"
done > command.phi.list
ParaFly -c command.phi.list -CPU 48

Blast分析后的注释结果中一条序列含有多种结果，后续的分析需要提取得到基因最可能的注释结果。

如“/media/aa/DATA2/bin/gene_annotation_from_SwissProt.pl”的操作，如下：

3）比对结果中筛选每个query的最佳subject

# 在(jcvi)下
conda activate jcvi

python -m jcvi.formats.blast best -n 1 $i.phi10.tab

批量操作

批量 pep70

for i in `cat /media/aa/DATA/SZQ2/bj/functional_annotation/70list.txt`
do
    echo "python -m jcvi.formats.blast best -n 1 $i.phi10.tab"
done > command.jcvi.list
ParaFly -c command.jcvi.list -CPU 48

 批量 pepmy

for i in `cat /media/aa/DATA/SZQ2/bj/functional_annotation/pepmylist.txt`
do
    echo "python -m jcvi.formats.blast best -n 1 $i.phi10.tab"
done > command.jcvi.list
ParaFly -c command.jcvi.list -CPU 48

4）复制并重命名

# 新建文件夹
mkdir best && cd best
# 复制
cp ../*.tcdb10.tab.best ./

首先将denovo组装出来的基因组序列预测得到蛋白文件*pep，与数据库中的宿主库和病原库分别比对。如果某些序列的病原库比对得分更高，就将这些序列与PHI数据库进行比对，找出病原相关基因。