xshell-7 本地构建blastp并快速计算POCP值202

*#### 目的： 进行保守蛋白相似性，即 percentage of conserved proteins (POCP) 从而区分属水平的距离。为此，首先需要构建本地blastp,并用现有的Python程序进行计算！

1. 本地blastp的构建：

• (注意：因为后续用的Python路径写的是E，所以必须是在E盘下！）

• 在ncbi下载对应的ncbi-blast-2.10.0+-win64，文件，比如：选择64位的exe, 点击并在E盘构建一个POCP文件夹，其中建立database文件，并安装上述程序！

• 下载地址：

ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

• 原因如下：

cd E:\POCP

图片.png

• 并创建三个文件夹：名称如下：

• Rawdata
• Database
• Result

• 不得错，区分大小写！

图片.png

• 把上面下载的exe安装到Database下：

• 重点来了：安装时生成的bin文件下的所有文件都要拿出来，直接放在POCP文件夹下！

图片.png image.png

---

2. 下载和安装pocp计算程序：

• 下载地址：

https://github.com/2015qyliang/POCP/commit/f2adfb205fb6d91a28d7afe308f209c175f266da#diff-5124dfd2e4a9ae1bf5b439932a2b65892f0ba108d5643c346e62ccc125109c64

• 或者：

https://github.com/2015qyliang/POCP

• 重点来了：改动的代码都是那个615py里的！其中的最终代码为：

print("hello world")
import os
import pandas as pd
from collections import OrderedDict

def Total(filename):
    os.chdir(rawFilePath)
    file = open(filename)
    text = file.read()
    proteinTotal = text.count('>')
    file.close()
    return proteinTotal


def aaTotal(filename):
    os.chdir(rawFilePath)
    file = open(filename)
    text = file.read()
    file.close()
    ltsp = text.split('>')
    aaDict = OrderedDict()
    for line in ltsp:
        if line != '':
            firstBreakIndex = line.index('\n')
            key = line[:firstBreakIndex].split(' ')[0]
            aaDict[key] = len(line[firstBreakIndex:])
        else:
            pass
    return aaDict


def creatdb():
    os.chdir(scriptPath)
    proc1 = r'.\makeblastdb.exe -dbtype prot -parse_seqids -in .\Rawdata\\'
    proc2 = r' -out .\Database\\'
    for fn in fnlist:
        process = proc1 + fn + proc2 + fn.split('.')[0]
        print('='*10,process)
        os.system(process)

# 定义流程化blastp函数
def blast():
    proc1 = r'.\blastp.exe -evalue 1e-5 -max_target_seqs 1 -outfmt 6'
    proc2 = r' -query .\Rawdata\\'
    proc3 = r' -db .\Database\\'
    proc4 = r' -out .\Result\\'
    for fn in fnlist:
        os.chdir(scriptPath)
        index = fnlist.index(fn)
        if index != len(fnlist)-1:
            A = fn
            B = fnlist[index + 1]
            # A2B
            processAB = proc1 + proc2 + A + proc3 + B.split('.')[0] + proc4 + A.split('.')[0] + 'TO' + B.split('.')[0] + '.tab'
            # B2A
            processBA = proc1 + proc2 + B + proc3 + A.split('.')[0] + proc4 + B.split('.')[0] + 'TO' + A.split('.')[0] + '.tab'
            print('='*10,processAB)
            os.system(processAB)
            print('='*10,processBA)
            os.system(processBA)
            os.chdir(resultPath)
            resultFileName1 = A.split('.')[0] + 'TO' + B.split('.')[0] + '.tab'
            resultFileName2 = B.split('.')[0] + 'TO' + A.split('.')[0] + '.tab'
            dframe1 = pd.read_csv(resultFileName1, sep = '\t')
            dframe2 = pd.read_csv(resultFileName2, sep = '\t')
            D1 = aaTotal(A)
            for i in range(dframe1.shape[0]):
                dframe1.iloc[i,11] = (dframe1.iloc[i,3])/(D1[dframe1.iloc[i,0]])
            cdframe1 = dframe1[dframe1.iloc[:,2] > 40]
            cdframe11 = cdframe1[cdframe1.iloc[:,11] > 0.5]
            C1 = cdframe11.shape[0]
            D2 = aaTotal(B)
            for i in range(dframe2.shape[0]):
                dframe2.iloc[i,11] = (dframe2.iloc[i,3])/(D2[dframe2.iloc[i,0]])
            cdframe2 = dframe2[dframe2.iloc[:,2] > 40]
            cdframe22 = cdframe2[cdframe2.iloc[:,11] > 0.5]
            C2 = cdframe22.shape[0]
            T1 = Total(A)
            T2 = Total(B)
            POCP = round((((C1 + C2)/(T1 + T2))*100),2)

            pocpResult.append(A.split('.')[0] + '\t' + B.split('.')[0] + '\t' + str(POCP) + '\n')
        else:
            pass


rawFilePath = r'D:\POCP\Rawdata'
scriptPath = r'D:\POCP'
resultPath = r'D:\POCP\Result'
fnlist = os.listdir(rawFilePath)
pocpResult = []
creatdb()
blast()
os.chdir(scriptPath)




from datetime import datetime

with open('POCP_Result_'+datetime.now().date().isoformat()+'.txt', 'w')as file:
    for line in pocpResult:
        file.write(line)
    file.close()
print('='*50)

• 下载后解压到POCP文件下，打开后20180615.py那个文件里的路径是默认E盘的，所以需要改路径！是r下的路径，改为D开头的，然后再处理！

rawFilePath = r'D:\POCP\Rawdata'
scriptPath = r'D:\POCP'
resultPath = r'D:\POCP\Result'

• 改行代码，只保留两位小数：一行只保留两位小数的代码：用round函数，先乘以100，再取2位小数：

POCP = round((((C1 + C2)/(T1 + T2))*100),2)

• 其中，文件名那行，我是已经写了代码，即输出结果的时候，可以直接输出为POCP-带当天日期！

from datetime import datetime

with open('POCP'+datetime.now().date().isoformat()+'.txt', 'w')as file:
    for line in pocpResult:
        file.write(line)
    file.close()
print('='*50)

图片.png

---

3. 下载细菌基因组对用的蛋白质序列！注意是蛋白质序列！！！！

• 比如：先到NCBI,搜索基因组,点击出来的基因组后的protein，下载后解压，得到faa文件，将其后缀改为fasta，即可！
• 注意：文件名之间不得有空格！所以最简单的方法就是命名如：strain1.fasta
• 同样的方法将两个想要比较的基因组的蛋白质序列进行比较，然后再进行最后一步，进行比较！

https://www.ncbi.nlm.nih.gov/genome/?term=Blautia+faecis

图片.png

4. 进行pocp的计算：

• 在Windows上，打开powershell： 将鼠标点击Windows窗口处，然后输入powershell即可，然后运行如下：

图片.png

• 运行之后进入上面的设置的路径：

1. 进入路径：

 cd E:\POCP\

2. 运行里面放着的程序(注意：事先确保电脑上装过Python，如果没有安装过，可以点击最下面那个下载并安装即可）：

• 计算pocp值：（计算可能需要五分钟左右，计算完了这一组就得删除这一组蛋白质序列，放上新的一对，才能进行下一组的计算！）

python .\pocpPipeline.py

• 运行结束页面：

  
  图片.png

• 修改后的代码给出的结果在如下txt文件下：

• 然后就在pocp文件夹下看到一个TXT结果，请注意：不是在result！

  
  图片.png
• 查看结果：0.4866850321395776，即48.67%

图片.png

---

• 改好的代码输出的结果氏这样的：

image.png