PSI-BLAST(Position-Specific Iterative Basic Local Alignment Search Tool)是blastp的一种，用于发现与查询序列更远距离相关的序列。
PSSM（position-specific scoring matrix），PSI-BLAST中要用到的中间结果，也被称为profile（配置文件）。PSSM捕获保持模式的对齐并将其存储为对齐中每个位置的分数矩阵 - 高度保守的位置接收高分数，弱保守位置接收分数接近零。该配置文件用于代替原始替换矩阵，以进一步搜索数据库以检测与PSSM指定的保护模式匹配的序列。
步骤
PSI-BLAST一般分为以下几个步骤：
1.用提交的序列在数据库中进行搜索，（与blastp相同）；
2.对初始的结果进行多序列比对，构建初始profile，即PSSM；
3.用PSSM作为query，搜索数据库，发现较远的相似性序列，并给出E-value；
4.构建新的PSSM，重复步骤3。这个过程不断迭代，直到没有符合要求(E-value threshold,比如0.005,有时可放宽至0.01)的新序列或者到达迭代次数要求。
运行
../ncbi-blast-2.2.25+/bin/psiblast -query ../protein/1A01_HUMAN.fasta -evalue .001 -inclusion_ethresh .002 -db ../databases/nr90-2012/nr90 -num_iterations 3 -seg yes -outfmt '7 std qseq sseq stitle' -out 1A01_HUMAN_100.output -max_target_seqs 100
../ncbi-blast-2.2.25+/bin/psiblast -query ../protein/1A01_HUMAN.fasta -evalue .001 -inclusion_ethresh .002 -db ../databases/nr90-2012/nr90 -num_iterations 3 -seg yes -outfmt '7 std qseq sseq stitle' -out 1A01_HUMAN_100.output -max_target_seqs 500
-query:搜索的序列；
-num_iterations：迭代次数；
-seg:是否SEG过滤；
-max_target_seqs：最大匹配数量，用100、500做了两次
-outfmt:输出格式设置，改了一下方便后续处理。
结果处理
outfmt7是带有注释行的tab文件，对于R比较友好，同时使用的格式转换的软件MView也支持outfmt7。
psiblast输出的output文件包含了所有迭代过程所得到的结果，我们需要的只有最后一次的搜索结果，因此需要过滤一下。
library(tidyverse)
text <- readLines("./1A01_HUMAN_100.output")
firstStr <- grep("#",str_sub(text,1,1))
start <- tail(firstStr,2)[1]+1
end <- tail(firstStr,2)[2]-1
cat(text[c(1:6,start:(end+1))],sep = "\n",file = "./1A01_HUMAN_last.output")
输出的1A01_HUMAN_last.output使用MView转换为aln文件。
MView
MView是一个命令行程序，用于提取和重新格式化序列数据库搜索或多重对齐的结果，可选择为Web页面布局添加HTML标记。 它还可以用作过滤器，以提取和转换搜索或对齐到常见格式。
安装见：https://desmid.github.io/mview/index.html
./mview/bin/mview -in blast 1A01_HUMAN_100_last.output -out aln > 1A01_HUMAN_100.aln
具体的使用没有细看，-in -out 两个参数设置一下就可以做一些基本的使用了。
fasta的提取需要用到reutils包
接上面：
id <- text[start:end] %>%
str_split("\t",simplify = T) %>%
.[,2] %>%
str_split("\\|",simplify = T) %>%
.[,2]
library(reutils)
reutils::epost(id,'protein') %>%
efetch(retmode = "text" , rettype = "fasta",outfile = "./1A01_HUMAN_100.fasta")
efetch()函数在id数量大于500时要添加outfile参数，否则会有报错，数量不多的时候可以如下：
reutils::epost(id,'protein') %>%
efetch(retmode = "text" , rettype = "fasta") %>%
content() %>%
cat(file = "./1A01_HUMAN_100.fasta")
E-utilities不使用R包，通过URL也可以得到结果，具体参考https://www.ncbi.nlm.nih.gov/books/NBK25500/
，或者使用web-server也可以。
参考：
http://yangli.name/2015/11/11/20151111psiblast/
https://www.ncbi.nlm.nih.gov/books/NBK2590/
https://www.ncbi.nlm.nih.gov/books/NBK25500/
