<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Building database
I have 6 SS WPS2 genomes, so I'll concatenate sequences into a single combined fasta file
Version 1.3
Step1:  OK - verified
Running Blast query
Step2:  Optimizing number of threads - No real difference
Step6:  Defining functions
Step7:  Making all_WPS-2 Database [ln3] again more scalable way
Step8:  Reading Raw input file and creating specific input file
Step9:  After iteration0 is done
Step10:  Outputting iteration 1
Step11:  Repeating what we've done
Step12:  adding those to the database
Step13:  Blasting assembled and unasembled contigs
Step14:  Parsing Blast Outputs
Trying pandas instead
Step15:  Adding "Metric", "Seq_size", "Seq_nt" fields
indeces of metagenome sequences and the recruited from blast dataframe are placed in lists of the same size
Step16:  Number of unique scaffolds found
Step17:  Evalue sorting
Step18:  Filtering out by evalue cutoff
Step19:  Efect of evalue on # of unique scaffolds found
Step20:  Doing some stats on the recruited data
Step21:  Retrieving sequences from dataframe list
Step22:  Adding Seq_nt and Size entries
Step24:  Indexing approach
Less memory intensive, close index after it's done
Step25:  Filtering Step
Step26:  Batch export to csv and multiple FASTA
<ASSISTANT_TASK:>
Python Code:
from Bio import SeqIO
import time
import os
import shutil
import pandas
#parameters
version = 'v1.4'
project_name = 'wps2_bl_metagenome'
e_value = 1e-20
iden = 95.0
metric = 50.0
!cat ../ss_genomes/AP_WPS-2_bacterium* > ../ss_genomes/all_AP_WPS-2_bacterium.fna
!makeblastdb -in ../ss_genomes/all_AP_WPS-2_bacterium.fna -dbtype nucl -title "combined_ss_WPS2" -out ../blast_db/combined_ss_WPS2/combined_ss_WPS2 -parse_seqids
Explanation: Building database
I have 6 SS WPS2 genomes, so I'll concatenate sequences into a single combined fasta file
Version 1.3:
Functional programming
Combined assembled and unassembled
End of explanation
import time
start =time.time()
!blastp -db test_pp_pmo -query 76969.assembled.faa -out assembl_contigs_vs_test_pp_pmo2.tab -evalue .00001 -outfmt 6 -num_threads 8
print "command was executed in %d seconds"%(time.time()-start)
Explanation: OK - verified
Running Blast query
End of explanation
ex_time=list()
num_thread=list()
for i in [1,2,4,6,8]:
    start =time.time()
    !blastp -db test_pp_pmo -query 76969.assembled.faa -out assembl_contigs_vs_test_pp_pmo2.tab -evalue .00001 -outfmt 6 -num_threads {i}
    elapsed=time.time()-start
    print 'number of threads: ', i
    print 'exec time: ', elapsed
    ex_time.append(elapsed)
    num_thread.append(i)
Explanation: Optimizing number of threads - No real difference
End of explanation
def parse_contigs_ind(f_name):
    
    Returns sequences index from the input files(s)
    remember to close index object after use
    
    handle = open(f_name, "rU")
    record_dict = SeqIO.index(f_name,"fasta")
    handle.close()
    return record_dict
#returning specific sequences and overal list
def retrive_sequence(contig_lst, rec_dic):
    
    Returns list of sequence elements from dictionary/index of SeqIO objects specific to the contig_lst parameter
    
    contig_seqs = list()
    #record_dict = rec_dic
    #handle.close()
    for contig in contig_lst:
        contig_seqs.append(rec_dic[contig].seq.tostring())
    return contig_seqs
def filter_seq_dict(key_lst, rec_dic):
    
    Returns filtered dictionary element from rec_dic according to sequence names passed in key_lst
    
    return { key: rec_dic[key] for key in key_lst }
def unique_scaffold_topEval(dataframe):
#returns pandas series object
    variables = list(dataframe.columns.values)
    scaffolds=dict()
    rows=list()
    for row in dataframe.itertuples():
        
        #if row[1]=='Ga0073928_10002560':
        if row[1] not in scaffolds:
            scaffolds[row[1]]=row
        else:
            if row[11]<scaffolds[row[1]][11]:
                scaffolds[row[1]]=row
    rows=scaffolds.values()
    #variables=['quid', 'suid', 'iden', 'alen', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits']
    df = pandas.DataFrame([[getattr(i,j) for j in variables] for i in rows], columns = variables)
    return df
Explanation: Defining functions
End of explanation
iteration = 0
blast_db_homeDir = "../blast_db/"
#!mkdir blast_db_homeDir
Explanation: Making all_WPS-2 Database [ln3] again more scalable way
End of explanation
#command_line file of reference sequences
reference_0 = "../ss_genomes/all_AP_WPS-2_bacterium.fna"
input_dir = "input_files"
if os.path.exists(input_dir):
    shutil.rmtree(input_dir)
try:
    os.mkdir(input_dir)
except OSError:
    raise
#returns 
records_0 = parse_contigs_ind(reference_0)
ref_out_0 = "input_files/reference0.fna"
with open(ref_out_0, "w") as handle:
    SeqIO.write(records_0.values(), handle, "fasta")
    #NO NEED TO CLOSE with statement will automatically close the file
records_0.close()
Explanation: Reading Raw input file and creating specific input file
End of explanation
output_1 = "../ss_genomes/iteration1.fna"
iteration1_all =cutoff_contigs['1e-20_cutoff_contigs']
print len(iteration1_all) 
len((iteration1_all[iteration1_all['iden']>=95.0]))
print len((iteration1_all[iteration1_all['iden']>=95.0]))
print len((iteration1_all[iteration1_all['Metric']>=50]))
iteration1_filtered = iteration1_all[(iteration1_all['iden']>=95.0)&(iteration1_all['Metric']>=50)]
len(iteration1_filtered)
names = iteration1_filtered['quid'].tolist()
records_0 = dict(parse_contigs_ind(input_0))
#all_records = parse_contigs_ind(assembled_contigs)
records_1 = filter_seq_dict(names, all_records)
print type(records_1)
records_1.update(records_0)
len(records_0.items())
len(records_1.items())
Explanation: After iteration0 is done:
filter assembled e-20: >=95 iden, >=50 metric
add records to output_1
End of explanation
#records_0 = parse_contigs_ind(input_0)
output_1_0 = "../ss_genomes/iteration1_0.fna"
output_1 = "../ss_genomes/iteration1_0.fna"
with open(output_1_0, "w") as handle:
  SeqIO.write(records_0.values(), handle, "fasta")
with open(output_1, "w") as handle:
  SeqIO.write(records_1.values(), handle, "fasta")
Explanation: Outputting iteration 1
End of explanation
type(records_1)
# output_0 = "../ss_genomes/iteration0.fna"
# with open(output_0, "w") as handle:
#   SeqIO.write(records_0.values(), handle, "fasta")
len(records_0)
Explanation: Repeating what we've done
End of explanation
(iteration1_all[iteration1_all['iden']>=95.0])&(iteration1_all[iteration1_all['alen']>=50])
dct1
in_db = '../ss_genomes/all_AP_WPS-2_bacterium.fna'
title = 'combined_ss_WPS2'
outfile = "../blast_db/combined_ss_WPS2/combined_ss_WPS2"
!makeblastdb -in {infile} -dbtype nucl -title "{title}" -out {outfile} -parse_seqids
Explanation: adding those to the database
End of explanation
import time
q_fname1 = './../IMG Data/76969.assembled.fna'
q_fname2 = './../IMG Data/76969.unassembled_illumina.fna'
out_name1= 'all_wps2_q_assembled76969.tab'
out_name2= 'all_wps2_q_unasembled76969.tab'
database = '../blast_db/combined_ss_WPS2/combined_ss_WPS2'
# start =time.time()
# !blastn -db {database} -query "{q_fname1}" -out {out_name1} -evalue .00001 -outfmt 6 -num_threads 8
# elapsed=time.time()-start
# print 'exec time: ', elapsed
# start =time.time()
# !blastn -db {database} -query "{q_fname2}" -out {out_name2} -evalue .00001 -outfmt 6 -num_threads 8
# elapsed=time.time()-start
# print 'exec time: ', elapsed
Explanation: Blasting assembled and unasembled contigs
End of explanation
del recruited_mg
recruited_mg = []
blast_cols = ['quid', 'suid', 'iden', 'alen', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits']
recruited_mg_1 = pandas.read_csv(out_name1 ,sep="\t", header=None)
recruited_mg_1.columns=['quid', 'suid', 'iden', 'alen', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits']
recruited_mg_2 = pandas.read_csv(out_name2 ,sep="\t", header=None)
recruited_mg_2.columns=['quid', 'suid', 'iden', 'alen', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits']
recruited_mg = [recruited_mg_1, recruited_mg_2]
len(recruited_mg[0])
len(recruited_mg[1])
Explanation: Parsing Blast Outputs
Trying pandas instead
End of explanation
#For linux use symlink
metagenome_1 ="../../pp_metagenome1/IMG Data/76969.assembled.fna"
metagenome_2 ="../../pp_metagenome1/IMG Data/76969.unassembled_illumina.fna"
#!Beware of closing all records enries eventually
all_records_mg_1=parse_contigs_ind(metagenome_1)
all_records_mg_2=parse_contigs_ind(metagenome_2)
all_records=[all_records_mg_1, all_records_mg_2]
start = time.time()
#wps2_eval_cutoff['Seq_nt']=retrive_sequence(contig_list, all_records)
for i in range(len(recruited_mg)):
    #cutoff_contigs[dataframe]=evalue_filter(cutoff_contigs[dataframe])
    recruited_mg[i]=unique_scaffold_topEval(recruited_mg[i])
    contig_list = recruited_mg[i]['quid'].tolist()
    recruited_mg[i]['Seq_nt']=retrive_sequence(contig_list, all_records[i])
    recruited_mg[i]['Seq_size']=recruited_mg[i]['Seq_nt'].apply(lambda x: len(x))
    recruited_mg[i]['Coverage']=recruited_mg[i]['alen'].apply(lambda x: float(x))/recruited_mg[i]['Seq_size']
    recruited_mg[i]['Metric']=recruited_mg[i]['Coverage']*recruited_mg[i]['iden']
    recruited_mg[i] = recruited_mg[i][['quid', 'suid', 'iden', 'alen','Coverage','Metric', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits','Seq_size', 'Seq_nt']]
    #all_records[i].close()# keep open if multiple iterations
    
elapsed=time.time()-start
print elapsed
# don't close just yet
# do after it is not needed
for i in range(len(recruited_mg)):
    all_records[i].close()
print  len(recruited_mg[0])
print len(recruited_mg[1])
Explanation: Adding "Metric", "Seq_size", "Seq_nt" fields
indeces of metagenome sequences and the recruited from blast dataframe are placed in lists of the same size: No of -m (metagenome) sequences
End of explanation
len(set(recruited_mg[0]['quid']))
len(set(recruited_mg[1]['quid']))
Explanation: Number of unique scaffolds found
End of explanation
#result = df.sort(['A', 'B'], ascending=[1, 0])\
# wps2_assembled_sort_eval = wps2_assembled.sort_values(by=['eval'], ascending=False)
# wps2_assembled_sort_eval
Explanation: Evalue sorting
End of explanation
# eval_cutoff = 1e-20
# wps2_assembled_filtered = wps2_assembled_sort_eval[wps2_assembled_sort_eval['eval']<eval_cutoff]
# wps2_assembled_filtered
Explanation: Filtering out by evalue cutoff
End of explanation
# #Doing a 
# unique_scaffolds = list()
# eval_cutoffs = [1e-10, 1e-20, 1e-30, 1e-40, 1e-50, 1e-60, 1e-70, 1e-80, 1e-90, 1e-100, 1e-110, 1e-120, 1e-130, 1e-140, 1e-150]
# #for e_val in eval_cutoffs:
# for e_val in eval_cutoffs:
#     wps2_assembled_filtered = wps2_assembled_sort_eval[wps2_assembled_sort_eval['eval']<e_val]
#     unique_scaffolds.append(len(set(wps2_assembled_filtered['quid'])))
# print eval_cutoffs
# print unique_scaffolds
# import matplotlib
# import matplotlib.pyplot as plt
# %matplotlib inline  
# #plt.semilogy(t, np.exp(-t/5.0))
# plt.plot()
# plt.semilogx(eval_cutoffs, unique_scaffolds, marker = 'o', basex=10)
# plt.title('Number of unique scaffolds as a function of BLAST e-Value')
# plt.xticks(eval_cutoffs)
# plt.ylabel('No. of unique scaffolds')
# plt.xlabel('BLAST e-Value')
# plt.grid(True)
# plt.ylim((0,12000))
# plt.show()
Explanation: Efect of evalue on # of unique scaffolds found
End of explanation
#assembled
ranges = range(5, 210, 15)
unique_scaffolds = list()
eval_ranges = [10**(-x) for x in ranges]
assembled_filtered_lst = list()
#for e_val in eval_cutoffs:
for i in range(len(eval_ranges)):
    assembled_filtered = recruited_mg[0][recruited_mg[0]['eval']<eval_ranges[i]]
    assembled_filtered_lst.append(assembled_filtered)
    #unique_scaffolds.append(len(set(wps2_assembled_filtered['quid'])))
    unique_scaffolds.append(len(set(assembled_filtered_lst[i]['quid'])))
print ranges
print eval_ranges
print unique_scaffolds
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline  
import matplotlib.pylab as pylab
params = {'legend.fontsize': 'x-large',
          'figure.figsize': (15, 5),
         'axes.labelsize': 'x-large',
         'axes.titlesize':'x-large',
         'xtick.labelsize':'x-large',
         'ytick.labelsize':'x-large'}
pylab.rcParams.update(params)
#plt.semilogy(t, np.exp(-t/5.0))
plt.plot()
plt.semilogx(eval_ranges, unique_scaffolds, marker = 'o', basex=10)
plt.title('Number of unique assembled scaffolds as a function of BLAST e-Value', fontsize=25)
plt.xticks(eval_ranges, fontsize=18)
plt.yticks(fontsize=18)
plt.ylabel('No. of unique scaffolds', fontsize=20)
plt.xlabel('BLAST e-Value', fontsize=20)
plt.grid(True)
plt.ylim((0,12000))
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(18.5, 10.5)
fig.savefig('assembled_fn_evalue.pdf', dpi=300)
plt.show()
Explanation: Doing some stats on the recruited data
End of explanation
# dataframe_names= [str(x)+"_cutoff_contigs" for x in eval_ranges]
# dataframe_names
# # cutoff = eval_ranges[3]
# # cutoff_contigs = wps2_assembled_filtered_lst[3]
# cutoff_contigs={}
# # #cutoff_contigs
# for i in range(len(eval_ranges)):
#     cutoff_contigs[dataframe_names[i]]=wps2_assembled_filtered_lst[i]
len(cutoff_contigs['1e-50_cutoff_contigs'])
print list(cutoff_contigs['1e-50_cutoff_contigs'].columns.values)
Explanation: Retrieving sequences from dataframe list
End of explanation
# from Bio import SeqIO
# assembled_contigs = "../IMG Data/76969.assembled.fna"
# handle = open(assembled_contigs, "rU")
# record_dict = SeqIO.to_dict(SeqIO.parse(handle,"fasta"))
# handle.close()
# rec = record_dict["Ga0073928_11111377"]
# rec
Explanation: Adding Seq_nt and Size entries
End of explanation
# start = time.time()
# all_records = parse_contigs_ind(assembled_contigs)
# print time.time()-start
# all_records.close()
# all_records['Ga0073928_10022750']
# all_records=parse_contigs_ind(assembled_contigs)
# start = time.time()
# #wps2_eval_cutoff['Seq_nt']=retrive_sequence(contig_list, all_records)
# for dataframe in cutoff_contigs:
#     #cutoff_contigs[dataframe]=evalue_filter(cutoff_contigs[dataframe])
#     contig_list = cutoff_contigs[dataframe]['quid'].tolist()
#     cutoff_contigs[dataframe]['Seq_nt']=retrive_sequence(contig_list, all_records)
#     cutoff_contigs[dataframe]['Seq_size']=cutoff_contigs[dataframe]['Seq_nt'].apply(lambda x: len(x))
#     cutoff_contigs[dataframe]['Coverage']=cutoff_contigs[dataframe]['alen'].apply(lambda x: float(x))/cutoff_contigs[dataframe]['Seq_size']
#     cutoff_contigs[dataframe]['Metric']=cutoff_contigs[dataframe]['Coverage']*cutoff_contigs[dataframe]['iden']
#     cutoff_contigs[dataframe] = cutoff_contigs[dataframe][['quid', 'suid', 'iden', 'alen','Coverage','Metric', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits','Seq_size', 'Seq_nt']]
# elapsed=time.time()-start
# print elapsed
# all_records.close()
# cutoff_contigs['1e-50_cutoff_contigs']
# cutoff_contigs['1e-50_cutoff_contigs'].loc[7444]['Seq_nt']
# !cutoff_contigs['1e-50_cutoff_contigs'].loc[7444]['Seq_nt'][:30] "../IMG Data/76969.assembled.fna"
# #verification by manual searching and extracting the sequence
# len(''.join(CCTGTTAGTGCTTCGGCGCCCGATACAATCTCGAGCAGCTCCTTGGTTATCGCGGCTTGA
# CGCGCGTTGTTCATAGCGATCGTCAGTTCATCGATCAATTTGCCGGCGTTTTCGGTGGCG
# TTGTTCATTGCAAGCAGCTGAGCTGCGTAGAATGACGCATCGGTCTCCAGCATCGCCGAA
# AATAGTGTGAACTCCAGATATTTCGGCAAGAGTTTTGAAAGTACGAACTCCGGCGAAGGC
# ACGATTTCGACCGCGCCGCGCGGACCCGTGGATGGAGCCGCCGCGGTTTCATGCTCGATC
# GGGACAAGCTGCCGTAGCTCGGGTCGCTGCACCATCATCGAGACATGTTTGGACGAAACG
# AGGATGATGTCGCCGATCTCTCCGGCGGTGAAATCGGCGGTGACGCTCTGCGCGAGTTCG
# TGCGCGGTTTGCAGTTTGGACGGCGCACTCAGCGGCCAGCCCGGCCGATCTCCTAGGCCC
# CAGCGGCGCACAGCGTTGCGAGCTTTTATGCCGACGGTGTAAAAACGCGCGTCCGCGTGC
# TGACGCCAGAAGCGCTCGGCCTCGCGAATGACGTTGGAGTTGAATGCGCCGGCCAG.split()))
# len(cutoff_contigs['1e-50_cutoff_contigs'].loc[7444]['Seq_nt'])
recruited_mg[0]
# wps2_eval_cutoff['Seq_size']=wps2_eval_cutoff['Seq_nt'].apply(lambda x: len(x))
# wps2_eval_cutoff = wps2_eval_cutoff[['quid', 'suid', 'iden', 'alen', 'mism', 'gapo', 'qsta', 'qend', 'ssta', 'send', 'eval', 'bits','Seq_size', 'Seq_nt']]
# wps2_eval_cutoff
Explanation: Indexing approach
Less memory intensive, close index after it's done: record_dict.close()
End of explanation
for i in range(len(recruited_mg)):
    recruited_mg[i]=recruited_mg[i][(recruited_mg[i]['iden']>=iden)&(recruited_mg[i]['Metric']>=metric)&(recruited_mg[i]['eval']<=e_value)]
print  len(recruited_mg[0])
print len(recruited_mg[1])
Explanation: Filtering Step
End of explanation
#wps2_eval_cutoff.to_csv(path_or_buf="wps2_eval_cutoff.csv", sep='\t')
prefix="iteration1_"+str(version)+"_"+str(e_value)+"_recruited_mg"
for i in range(len(recruited_mg)):
    records= []
    outfile1 = prefix+str(i)+".csv"
#     try:
#         os.remove(outfile1)     
#     except OSError:
#         pass
    recruited_mg[i].to_csv(outfile1, sep='\t')
    ids = recruited_mg[i]['quid'].tolist()
    #if len(ids)==len(sequences):
    for j in range(len(ids)):
        records.append(all_records[i][ids[j]])
    outfile2 = prefix+str(i)+".fasta" 
#     try:
#         os.remove(outfile2)        
#     except OSError:
#         pass
    with open(outfile2, "w") as output_handle:
        SeqIO.write(records, output_handle, "fasta")
        
    #print prefix+dataframe_name    
Explanation: Batch export to csv and multiple FASTA
End of explanation
<END_TASK>