SWISS-PROT比对

第一步: 在uniprot下载UniProt 上植物dat格式的注释文件。

wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_plants.dat.gz
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_trembl_plants.dat.gz

将两个dat合并到成一个文件

zcat uniprot_sprot_plants.dat.gz uniprot_trembl_plants.dat.gz > uniprot_plants.dat

第二步: 从dat中提取fasta序列

dat=uniprot_plants.dat
awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' $1 > ${1%%.dat}_AC.fasta

第三步: 建立DIAMOND或NCBI BLAST+索引

diamond makedb --in uniprot_plants_AC.fasta -d uniprot_plants_AC

第四步: 使用DIAMOND或NCBI BLAST+进行比对

diamond blastp -d /data/database/UniProt-Plant/uniprot_plants_AC.dmnd -q proteins.fasta --evalue 1e-5 > blastp.outfmt6

第五步: 从DIMAMOND或NCBI BLAST+的比对结果中筛选每个query的最佳subject

python -m jcvi.formats.blast best -n 1 blastp.outfmt6

第六步: 使用add_annotation_from_dat.py(代码在GitHub上)根据blastp输出从dat中提取GO/KEGG/同源基因。运行在Python2/3环境中，需要安装BioPython

python ~/myscripts/add_annotation_from_dat.py blastp.outfmt6.best /data/database/UniProt-Plant/uniprot_plants.dat

之后会输出swiss_annotation.tsv， 输出信息包括如下几列

    gene id
    uniprot accession
    identity
    homology species
    EnsemblPlants
    GO ID
    GO component, CC/MF/BP
    evidence

————————————————
版权声明：本文为CSDN博主「徐洲更hoptop」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/u012110870/article/details/115511924