Python notebook to obtain the HMM-based probability distance matrix¶
In [1]:
# load libraries
import pandas as pd
import hhsuitedb as hhdb
import glob, os, subprocess, shlex
from Bio import SeqIO
import csb
from csb.bio.io.hhpred import HHOutputParser
from matplotlib import pyplot as plt
import numpy as np
import networkx as nx
import pickle
import tempfile
from ete3 import PhyloTree
import seaborn as sns
In [3]:
# Path to the UniClust database
uniclust = '/home/cactuskid13/mntpt/HHBLITsdb/uniclust30_2018_08/uniclust30_2018_08'
qdir = 'orco/'
runName ='ORCO'
ncores = 8
# path to all query files
queries = glob.glob(qdir + 'queries/*.fasta')
In [4]:
# Functions to run programs from the HHsuite software
def runHHblits( aln , name, path , outdir, db , iterations , ncores , runName='' , SS= False , ohhm = False , verbose = True , Z = 2000 , B = 2000 , xargs = ''):
if verbose == True:
print( [aln , name, path , outdir, db , iterations , ncores , runName] )
outhhr= outdir+name+runName+".hhr"
args = path + ' -cpu '+ str(ncores) +' -d ' + db + ' -i ' + aln +' -o '+ outhhr + ' -n ' + str(iterations) + ' -B '+ str(B) + ' -Z ' + str(Z) +' '+ xargs
if SS == True:
args += ' -ssm 2 -ssw .5 '
if ohhm == True:
outa3m = outdir+name+runName+'.hhm'
args += ' -ohhm ' + outa3m
else:
outa3m = None
if verbose == True:
print(args)
args = shlex.split( args)
p = subprocess.run( args )
return p , [outhhr,outa3m]
def runHHmake( aln , name, path = 'hhmake' , outdir='./', verbose = False, SS = False):
if verbose == True:
print( [aln , name, path , outdir] )
outhhm= outdir+name+".hhm"
args = path + ' -i '+ aln +' -o '+ outhhm + ' -M 50'
if SS == True:
#todo : make ss prediction here
pass
args = shlex.split(args)
print(args)
p = subprocess.Popen(args )
return p , [outhhm]
def hhrparse(hhr , coverage , proba ):
profile = HHOutputParser(alignments=False).parse_file(hhr)
qname = profile.query_name
for hit in profile:
proba = hit.probability
i = hit.id
In [8]:
# create output directories, as well as query files
queries = glob.glob(qdir + '*.fasta')
if not os.path.exists(qdir +'queries'):
os.mkdir(qdir + 'queries')
if not os.path.exists(qdir +'HHfiles'):
os.mkdir(qdir + 'HHfiles')
if not os.path.exists(qdir +'ALLVSALL'):
os.mkdir(qdir + 'ALLVSALL')
for qfile in queries:
print(qfile)
fasta = SeqIO.parse(qfile, 'fasta')
for s in fasta:
SeqIO.write( [s] , qdir + 'queries/' + s.id + '.fasta' , 'fasta' )
In [10]:
# run HHblits on each query
print(queries)
for q in queries:
p,output = runHHblits(q , name = q.split('.')[0].split('/')[-1] + 'Profphylo' , path= 'hhblits ' , outdir = qdir+'HHfiles/' , db = uniclust , iterations= 3 , ncores = ncores , ohhm = True, verbose = True , runName=runName, xargs = ' -mact .5')
In [18]:
# We can assign an HHblits score to each hit
newhits = {}
hhfiles = glob.glob(qdir+'HHfiles/*.hhr')
print(hhfiles)
for file in hhfiles:
results = HHOutputParser(alignments=False).parse_file(file)
for hit in results:
newhits.update( { hit.id.split('|')[1].strip() : { 'prob':hit.probability, 'len' : hit.length , 'query':results.query_name }} )
df = pd.DataFrame.from_dict(newhits, orient='index')
#adhoc score = len * proba
df['score'] = df['prob'] * df['len']
#sort the data frame and print
print(df.sort_values( 'score' , ascending = False) )
df.to_csv( './orco/orco_hhblits_results.csv' )
In [29]:
#plot the scores and lengths of hits
newhitsdf = pd.DataFrame.from_dict(newhits, orient= 'index')
newhitsdf.hist()
plt.show()
#filter out poor matches
filterdf = newhitsdf[newhitsdf.prob >.9]
filterdf = filterdf[newhitsdf.len >200]
filterdf.hist()
plt.show()
print(filterdf)
In [7]:
#make a db with HMMS
if os.path.exists(qdir+runName+'profilephylo_v2_hhm_db*'):
os.rm(qdir+runName+'profilephylo_v2_hhm_db*')
import hhsuitedb
hhsuitedb.add_new_files( qdir+'HHfiles/*.hhm' , "hhm" , qdir+runName+'profilephylo_v2_hhm_db')
In [8]:
#run all v all comparison
hhms = glob.glob( qdir+'HHfiles/*.hhm' )
print(hhms)
for q in hhms:
p,output = runHHblits(q , name = q.split('.')[0].split('/')[-1] + 'allvall' , path= 'hhsearch ' , outdir = qdir+'ALLVSALL/' , db = qdir+runName+'profilephylo_v2_hhm_db' , iterations= 1 , ncores = ncores , ohhm = False, verbose = True , runName= 'test' , xargs = ' -mact .3')
In [9]:
# parse hhr files and make dist kernel
def cleanID(ID):
if '|PDBID' in ID:
ID = ID.split('|PDBID')[0]
elif '|' in ID:
ID = ID.split('|')[0].split()[0]
elif '.' in ID:
ID = ID.split('.')[0]
elif '/' in ID:
ID = ID.split('/')[0]
elif 'hmmercut' in ID:
ID = ID.split('hmmercut')[0]
ID = ID.split()[0]
ID = ID.strip()
return ID
def HHSearch_parseTo_DMandNX(hhrs):
clusternames = []
for i,hhr in enumerate(hhrs):
print(hhr)
profile = HHOutputParser(alignments=False).parse_file(hhr)
if cleanID(profile.query_name) not in clusternames:
print(cleanID(profile.query_name))
clusternames.append(cleanID(profile.query_name))
print(clusternames)
evalDM = np.ones( (len(clusternames),len(clusternames) ))
pvalDM = np.ones( (len(clusternames),len(clusternames) ))
scoreDM = np.zeros( (len(clusternames),len(clusternames) ))
SSDM = np.zeros( (len(clusternames),len(clusternames) ))
probaDM = np.zeros( (len(clusternames),len(clusternames) ))
lenDM = np.ones( (len(clusternames),len(clusternames) ))
NX = nx.Graph()
for i,hhr in enumerate(hhrs):
protlist = []
profile = HHOutputParser(alignments=False).parse_file(hhr)
for hit in profile:
DMscore = float(hit.evalue)
proba = hit.probability
if 'anchor' not in hit.id and 'anchor' not in profile.query_name:
i = clusternames.index(cleanID(hit.id))
j = clusternames.index(cleanID(profile.query_name))
if hit.evalue < evalDM[i,j]:
evalDM[i,j] = hit.evalue
evalDM[j,i] = evalDM[i,j]
if hit.pvalue < pvalDM[i,j]:
pvalDM[i,j] = hit.pvalue
pvalDM[j,i] = pvalDM[i,j]
if scoreDM[i,j] < hit.score:
scoreDM[i,j] = hit.score
scoreDM[j,i] = scoreDM[i,j]
if SSDM[i,j] < hit.ss_score:
SSDM[i,j] = hit.ss_score
SSDM[j,i] = SSDM[i,j]
if probaDM[i,j] < hit.probability:
probaDM[i,j] = hit.probability
probaDM[j,i] = probaDM[i,j]
#use smallest of the two prots
if lenDM[i,j] == 1 or lenDM[i,j] > hit.qlength:
lenDM[i,j] = hit.qlength
lenDM[j,i] = lenDM[i,j]
if hit.id != profile.query_name :
NX.add_edge( hit.id , profile.query_name )
NX[hit.id][profile.query_name]['score']= hit.score
return probaDM, evalDM ,pvalDM, lenDM , scoreDM, SSDM, NX , clusternames
In [14]:
print(glob.glob(qdir+'ALLVSALL/*.hhr'))
probaDM, evalDM ,pvalDM, lenDM , scoreDM, SSDM, NX , clusternames = HHSearch_parseTo_DMandNX(glob.glob(qdir+ 'ALLVSALL/*.hhr'))
In [21]:
g = sns.clustermap( probaDM , xticklabels=clusternames , yticklabels=clusternames , figsize = (20,20) )
