*#### 目的： 进行保守蛋白相似性，即 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
- 运行之后进入上面的设置的路径：
- 进入路径：
cd E:\POCP\
- 运行里面放着的程序(注意：事先确保电脑上装过Python，如果没有安装过，可以点击最下面那个下载并安装即可）：
- 计算pocp值：（计算可能需要五分钟左右，计算完了这一组就得删除这一组蛋白质序列，放上新的一对，才能进行下一组的计算！）
python .\pocpPipeline.py
-
运行结束页面：
图片.png
-
修改后的代码给出的结果在如下txt文件下：
-
然后就在pocp文件夹下看到一个TXT结果，请注意：不是在result！
图片.png
- 查看结果：0.4866850321395776，即48.67%
图片.png
-
改好的代码输出的结果氏这样的：
image.png
网友评论