数据库下载与构建

下载

wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz

构建

解压缩
```
gunzip -d uniprot_sprot.fasta.gz
```

构建 blast + 数据库

makeblastdb -in uniprot_sprot.fasta -dbtype prot -out uniprot_sprot -parse_seqids

构建 DIAMOND 数据库

diamond makedb --in uniprot_sprot.fasta -d uniprot_sprot_diamond

比对

blastp 蛋白比对

blastp -query F01.faa -out F01.swissprot -db /new_data/hualin/db/uniprot_sprot -outfmt 6 -num_threads 30 -evalue 1e-5

diamond 蛋白比对

单个基因组对比

diamond blastp -d /new_data/hualin/db/uniprot_sprot_diamond -q F01.faa -e 1e-5 -f 6 -o F01.diamond -k 1 --sensitive -p 30 --query-cover 50

多个个基因组对比

不会 shell 没办法，写 Perl 脚本 (run_diamond.pl) 来完成。

#!/usr/bin/perl
   use strict;
   use warnings;
# Author: Liu hualin
# Date: Sep 28, 2021

   my @faa = glob("*.faa");# 读取所有后缀为“.faa”的文件，可以自己更改
   foreach  (@faa) {
    $_=~/(.+).faa/;
    my $out = $1 . ".diamond";
    # 将/new_data/hualin/db/uniprot_sprot_diamond换成自己的数据库路径; -p表示线程数，在笔记本上用6个即可
    system("diamond blastp -d /new_data/hualin/db/uniprot_sprot_diamond -q $_ -e 1e-5 -f 6 -o $out -k 1 --sensitive -p 30 --query-cover 50");
   }

将上述代码复制到文件中，命名为 “run_diamond.pl”，将其和序列文件放在同一目录下，并在终端中输入如下命令，完成分析：

perl run_diamond.pl

将比对结果 mapping 至其他数据库

打开网址 https://www.uniprot.org/uploadlists/, 上传比对上的 swissprot ID，可以将比对结果转换为诸如 KEGG 等其他数据库的 ID。个人感觉不是很好用。

我们可以把 mapping 文件下载下来，自己写脚本来提取信息，虽然麻烦些，但得到的更多。

下载 mapping 文件

wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/idmapping_selected.tab.gz

解压缩
```
gunzip -d idmapping_selected.tab.gz
```

写脚本提取对应信息

Diamond 比对的结果文件内容如下，第一列是自己的氨基酸序列 ID，第二列是 SwissProt 数据库中序列的 ID，而我们真正需要的是第二列中两个竖线中间的内容，在稍后的脚本中将通过正则表达式把它给揪出来。

F01_00001	sp|Q73G44|MDH_WOLPM	47.2	72	38	0	10	81	243	314	9.55e-16	72.8
   F01_00003	sp|D9PU00|TFRA_METTM	41.3	569	301	7	7	574	4	540	4.89e-131	397
   F01_00004	sp|P9WN88|FRDB_MYCTO	32.7	208	118	6	19	215	23	219	3.84e-28	110
   F01_00005	sp|Q021N6|SUCC_SOLUE	62.8	384	141	2	1	383	1	383	1.45e-155	446

开始写脚本，保存为 “run_mapping.pl”。

#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 28, 2021

my %maps;
my @diaout = glob("*.diamond");# 读取所有的diamond比对后的输出文件
foreach  (@diaout) {
	$_=~/(\S+).diamond/;
	my %hash;
	open IN, "$_" || die;
	while (<IN>) {
		chomp;
		$_=~s/[\r\n]+//g;
		my @lines = split;
		$lines[1]=~/.+\|(.+)\|.+/;
		$hash{$1}++;
	}
	close IN;

	open IN, "idmapping_selected.tab" || die;
	while (<IN>) {
		chomp;
		$_=~s/[\r\n]+//g;
		my @lines = split;
		if (exists $hash{$lines[0]}) {
			$maps{$lines[0]} = $_;
		}
	}
	close IN;
}

my @diaout2 = glob("*.diamond");# 读取所有的diamond比对后的输出文件
foreach  (@diaout2) {
	$_=~/(\S+).diamond/;
	my $out = $1 . ".mapped";
	open IN, "$_" || die;
	open OUT, ">$out" || die;
	print OUT "qseqid\tsseqid\tpident\tlength\tmismatch\tgapopen\tqstart\tqend\tsstart\tsend\tevalue\tbitscore\tUniProtKB-AC	UniProtKB-ID	GeneID (EntrezGene)	RefSeq	GI	PDB	GO	UniRef100	UniRef90	UniRef50	UniParc	PIR	NCBI-taxon	MIM	UniGene	PubMed	EMBL	EMBL-CDS	Ensembl	Ensembl_TRS	Ensembl_PRO	Additional PubMed\n";
	while (<IN>) {
		chomp;
		$_=~s/[\r\n]+//g;
		my @lines = split;
		$lines[1]=~/.+\|(.+)\|.+/;
		if (exists $maps{$1}) {
			print OUT $_ . "\t" . $maps{$1} . "\n";
		}else {
			print OUT $_ . "\n";
		}
	}
	close IN;
	close OUT;
}

将脚本与 diamond 比对的结果文件以及下载的 mapping 文件放在同一目录下，在终端里输入如下命令即可得到 mapping 后的结果：

perl run_mapping.pl

GO 注释

从 map 后的文件中提取基因 ID 和 GO number，各列以制表符分隔，没有 GO 注释的只输出 gene ID。

准备脚本，命名为 “get_GO.pl”，与上一步生成的 “*.mapped” 文件放在同一目录下。

#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 28, 2021

my @mapped = glob("*.mapped");
foreach  (@mapped) {
	$_=~/(.+).mapped/;
	open IN, "$_" || die;
	my $out = $1 . ".GO";
	open OUT, ">$out" || die;
	<IN>;
	while (<IN>) {
		chomp;
		$_=~s/[\r\n]+//g;
		my @lines = split /\t/;
		print OUT $lines[0];
		if ($lines[18]=~/.+\; /) {
			my @terms = split /\; /, $lines[18];# 18代表文件的第19列
			print OUT "\t" . join("\t", @terms) . "\n";
		}elsif ($lines[18]=~/\S+/) {
			print OUT "\t" . $lines[18] . "\n";
		}else {
			print OUT "\n";
		}
	}
	close IN;
	close OUT;
}