<commit_before>import os, sys, re, gc
import optparse
import shutil
import subprocess
import pandas
import h5py
import numpy
import operator
from pyteomics import mzml
#from pyopenms import *
from Bio import SeqIO
from subprocess import Popen
import uniprot as uni
#####################################
#This is a script to combine the outputs of many tools into
#an Rdata image!
#
#VERSION 0.985A
version="0.985A"
#DATE: 7/25/2017
date="7/25/2017"
#####################################
print "-----------------------------------------------------------------------"
print "Welcome to the Rdata merger for Galaxy, Wohlschlegel Lab UCLA"
print "Written by <NAME>"
print "Version: ",version
print "Date: ",date

basedir=os.getcwd()


####################################
#Argument parsing! So much fun!
#We'll use OptParse even though some
#people really rave about argparse...
#
#
# NB: With Optparse, if an option is
#     not specified, it will take a
#     value of None
####################################

parser = optparse.OptionParser()
#parser.add_option("--folder",action="store",type="string",dest="operation_folder",default=".")
parser.add_option("--skyline_chromatograms",action="store",type="string",dest="skyline_chromatograms")
parser.add_option("--skyline_boundaries",action="store",type="string",dest="skyline_boundaries")
parser.add_option("--nsaf_table",action="store",type="string",dest="nsaf_table")
parser.add_option("--spc_table",action="store",type="string",dest="spc_table")
parser.add_option("--peptide_nsaf_table",action="store",type="string",dest="peptide_nsaf_table")
parser.add_option("--peptide_spc_table",action="store",type="string",dest="peptide_spc_table")
parser.add_option("--saint_table",action="store",type="string",dest="saint_table")
parser.add_option("--fido_roc",action="store",type="string",dest="fido_roc")
parser.add_option("--fasta",action="store",type="string",dest="fasta")
parser.add_option("--job_id",action="store",type="string",dest="job_id")
parser.add_option("--mzml",action="store",type="string",dest="mzml")
parser.add_option("--msstats_comparison",action="store",type="string",dest="msstats_comparison")
parser.add_option("--msstats_condition",action="store",type="string",dest="msstats_condition")
parser.add_option("--msstats_quantification",action="store",type="string",dest="msstats_quantification")
parser.add_option("--msstats_skyline_input",action="store",type="string",dest="msstats_skyline_input")
parser.add_option("--RScriptOutput",action="store",type="string",dest="RScriptOutput")
parser.add_option("--RDataOutput",action="store",type="string",dest="RDataOutput")
parser.add_option("--MS2Output",action="store",type="string",dest="MS2Output")
parser.add_option("--MS2Output_id",action="store",type="string",dest="MS2Output_id")
parser.add_option("--experiment_file",action="store",type="string",dest="experiment_file")
parser.add_option("--quant_experiment_file",action="store",type="string",dest="quant_experiment_file")
parser.add_option("--qual_experiment_file",action="store",type="string",dest="qual_experiment_file")
parser.add_option("--psm_table",action="store",type="string",dest="psm_table")
parser.add_option("--filter_decoys",action="store_true",dest="filter_decoys")
parser.add_option("--zlib",action="store_true",dest="zlib")
parser.add_option("--renameProteinType",action="store",type="string",default="accToGene",dest="renameProteinType") #Rename proteins from accToGene, accToProtein, geneToProtein (for TGGT/ToxoDB)
parser.add_option("--analysis_type",action="store",type="string",dest="analysis_type")

(options,args) = parser.parse_args()

print "Now we're going to prepare the R script for Rdata image generation"

#group_information = pandas.read_csv(options.qual_experiment_file,sep='\t')


arguments=['python','/galaxy-central/tools/wohl-proteomics/wohl_skyline/job_history_runner.py','--job_id',str(options.job_id),"--output_file","job_history.json"]
print "Running job:"
print " ".join(arguments)
job_process=subprocess.Popen(arguments,stdout=subprocess.PIPE, stderr=subprocess.PIPE)

#### We'll read in the Skyline chromatograms data, group on the File column, and remove decoys
#decoys are duplicates of the following columns : fileName ModifiedSequence PrecursorCharge FragmentIon ProductCharge IsotopeLabelType
#let's get to it.
#
# NB NB NB NB: This code assumes that 1. Decoys come after targets. AND 2. Decoys are 1:1 with targets.  There are no targets without decoys, and no target transitions without decoy transitions.
if options.analysis_type=="lfq":
    chromatogram_df=pandas.read_csv(options.skyline_chromatograms,sep="\t")
    #chromatogram_df=pandas.read_csv(options.skyline_chromatograms,sep="\t",low_memory=False)

    if options.filter_decoys:
        #groups=chromatogram_df.groupby(by='FileName')
        #groups=chromatogram_df.groupby(by=['FileName','PeptideModifiedSequence','PrecursorCharge','FragmentIon','ProductCharge','IsotopeLabelType'])
        groups=chromatogram_df.groupby(by=['FileName','PeptideModifiedSequence','IsotopeLabelType'])
        filtered_df=[]
        first=True
        for groupKey,group in groups:
            filtered_df.append(group.head(n=(len(group)/2)))
        filtered_dataframe=pandas.concat(filtered_df)
    else:
        filtered_dataframe=chromatogram_df
    filtered_dataframe.to_csv("filtered_chromatograms.tsv",sep="\t",index=False)
    print "Finished outputting the filtered chromatograms..."






with open(options.fasta) as fasta_file:
    sequences = []
    identifiers = []
    descriptions = []
    lengths=[]
    for seq_record in SeqIO.parse(fasta_file, 'fasta'):  # (generator)
        identifiers.append(seq_record.id)
        sequences.append(str(seq_record.seq))
        descriptions.append(seq_record.description.replace("\t",","))
        lengths.append(len(seq_record.seq))
fasta_df=pandas.DataFrame({'id':identifiers,'description':descriptions,'Sequence':sequences})

fasta_df.to_csv("fasta_df.tsv",sep="\t",index=False)

#We're going to read in the MS2 scans via Pyteomics, and then store them in a pandas df.
#Then we'll write them to a TSV file <----- New version implements HDF5 instead to keep it out of memory!
mzml_files=options.mzml.split("___")
mzml_file_scans=[]
ms2fail=False
first=True
h5file=h5py.File(options.MS2Output,'w')

#group_information    
for each_file in mzml_files:
    if first:
        first=False
    else:
        h5file.close()
        gc.collect()
        h5file=h5py.File(options.MS2Output,'a')
    #this_group=h5file.create_group(str(group_information.loc[group_information['Original File Name']==each_file.rsplit(".",1)[0],'Crux File Integer'].iloc[0]))  #This makes groups by crux integer
    this_group=h5file.create_group(each_file)
    try:
        print "Reading in ",each_file
        mzml_reader=mzml.read(each_file,iterative=True)
        for each_scan in mzml_reader:
            if each_scan['ms level']==2:
                this_dataset=this_group.create_dataset(str(each_scan['index']),compression="gzip",compression_opts=9,dtype="float32",data=numpy.column_stack((each_scan['m/z array'],each_scan['intensity array'])).T)
                this_dataset.attrs["scan_index"]=each_scan['index']
        print "done with it........\n"
        del mzml_reader
        del each_scan
        h5file.flush()
    except:
        print "Failed to read in file",each_file
        ms2fail=True
#if not ms2fail:
    #scans_df=pandas.DataFrame.from_dict(mzml_file_scans)
    #scans_df.to_csv("ms2_scans.tsv",sep="\t",index=False)
    #del scans_df

h5file.close()
#job_process.wait()
#out, err = job_process.communicate()
#print out,err
#print "Finished the job processing."


job_process.wait()
out, err = job_process.communicate()
print out,err
print "Finished the job processing."


print "reading PSM table:"
try:
    psm_table=pandas.read_csv(options.psm_table,sep='\t',engine="python")
    #psm_table=pandas.read_csv(options.psm_table,sep='\t',low_memory=True)
except:
    psm_table=pandas.read_csv(options.psm_table,sep='\t',low_memory=False) #Sometimes causes segfaults....
print "finsihed reading psm table"

### The code below is a mixture of things that we need in LFQ and Qual analysis...
### It is meant to query uniprot, and generate name mappings for the MSstats csv table, and the PSM table.
print "Handling the renaming!" #This is done only for both the PSM_table and the MSstats outputs...
if options.analysis_type=="lfq":
    combined_results=pandas.read_csv(options.msstats_comparison,sep=',',low_memory=False)
    #combined_results=pandas.read_csv(options.msstats_comparison,sep=',',index_col=False)
    #print combined_results
    combined_results['backup']=combined_results['Protein']

protein_series=psm_table['protein id']
uni_mapping_dict={}
uniprot_pattern=r'([A-NR-Z][0-9][A-Z][A-Z0-9][A-Z0-9][0-9]|[OPQ][0-9][A-Z0-9][A-Z0-9][A-Z0-9][0-9])'
protein_extracted=protein_series.str.extractall(uniprot_pattern)[0].unique()#for the PSM table...
protein_full_name_series_list=protein_series.str.split(",",expand=False) #This will be used later to properly adjust positioning of AA's...


if "accToGene" in options.renameProteinType or options.renameProteinType is None: #This means we'll take the acc's and replace with GENE NAMES
    if options.analysis_type=="lfq":
        combined_results['uniprot_acc']=combined_results['Protein'].str.extract(uniprot_pattern,expand=False)
    #uni_mapping_dict_combined=uni.map(list(combined_results['uniprot_acc'].unique()),f='ACC',t='GENENAME')
    try:
        uni_mapping_dict_combined=uni.map(list(protein_extracted),f='ACC',t='GENENAME')
    except:
        uni_mapping_dict_combined={} #failure....
    for each_acc in uni_mapping_dict_combined:
        catch=uni_mapping_dict_combined[each_acc]
        if len(catch)>0:
            uni_mapping_dict[str(each_acc)]=str(catch.pop())
    second_query=[]
    for each_acc in protein_extracted:<commit_msg>Pandas C-engine chokes on larger files from MilkyWay.  This is a non-issue when using the python engine.  We have switched this tool to the python engine.<commit_after>import os, sys, re, gc
import optparse
import shutil
import subprocess
import pandas
import h5py
import numpy
import operator
from pyteomics import mzml
#from pyopenms import *
from Bio import SeqIO
from subprocess import Popen
import uniprot as uni
#####################################
#This is a script to combine the outputs of many tools into
#an Rdata image!
#
#VERSION 0.985A
version="0.985A"
#DATE: 7/25/2017
date="7/25/2017"
#####################################
print "-----------------------------------------------------------------------"
print "Welcome to the Rdata merger for Galaxy, Wohlschlegel Lab UCLA"
print "Written by <NAME>"
print "Version: ",version
print "Date: ",date

basedir=os.getcwd()


####################################
#Argument parsing! So much fun!
#We'll use OptParse even though some
#people really rave about argparse...
#
#
# NB: With Optparse, if an option is
#     not specified, it will take a
#     value of None
####################################

parser = optparse.OptionParser()
#parser.add_option("--folder",action="store",type="string",dest="operation_folder",default=".")
parser.add_option("--skyline_chromatograms",action="store",type="string",dest="skyline_chromatograms")
parser.add_option("--skyline_boundaries",action="store",type="string",dest="skyline_boundaries")
parser.add_option("--nsaf_table",action="store",type="string",dest="nsaf_table")
parser.add_option("--spc_table",action="store",type="string",dest="spc_table")
parser.add_option("--peptide_nsaf_table",action="store",type="string",dest="peptide_nsaf_table")
parser.add_option("--peptide_spc_table",action="store",type="string",dest="peptide_spc_table")
parser.add_option("--saint_table",action="store",type="string",dest="saint_table")
parser.add_option("--fido_roc",action="store",type="string",dest="fido_roc")
parser.add_option("--fasta",action="store",type="string",dest="fasta")
parser.add_option("--job_id",action="store",type="string",dest="job_id")
parser.add_option("--mzml",action="store",type="string",dest="mzml")
parser.add_option("--msstats_comparison",action="store",type="string",dest="msstats_comparison")
parser.add_option("--msstats_condition",action="store",type="string",dest="msstats_condition")
parser.add_option("--msstats_quantification",action="store",type="string",dest="msstats_quantification")
parser.add_option("--msstats_skyline_input",action="store",type="string",dest="msstats_skyline_input")
parser.add_option("--RScriptOutput",action="store",type="string",dest="RScriptOutput")
parser.add_option("--RDataOutput",action="store",type="string",dest="RDataOutput")
parser.add_option("--MS2Output",action="store",type="string",dest="MS2Output")
parser.add_option("--MS2Output_id",action="store",type="string",dest="MS2Output_id")
parser.add_option("--experiment_file",action="store",type="string",dest="experiment_file")
parser.add_option("--quant_experiment_file",action="store",type="string",dest="quant_experiment_file")
parser.add_option("--qual_experiment_file",action="store",type="string",dest="qual_experiment_file")
parser.add_option("--psm_table",action="store",type="string",dest="psm_table")
parser.add_option("--filter_decoys",action="store_true",dest="filter_decoys")
parser.add_option("--zlib",action="store_true",dest="zlib")
parser.add_option("--renameProteinType",action="store",type="string",default="accToGene",dest="renameProteinType") #Rename proteins from accToGene, accToProtein, geneToProtein (for TGGT/ToxoDB)
parser.add_option("--analysis_type",action="store",type="string",dest="analysis_type")

(options,args) = parser.parse_args()

print "Now we're going to prepare the R script for Rdata image generation"

#group_information = pandas.read_csv(options.qual_experiment_file,sep='\t')


arguments=['python','/galaxy-central/tools/wohl-proteomics/wohl_skyline/job_history_runner.py','--job_id',str(options.job_id),"--output_file","job_history.json"]
print "Running job:"
print " ".join(arguments)
job_process=subprocess.Popen(arguments,stdout=subprocess.PIPE, stderr=subprocess.PIPE)

#### We'll read in the Skyline chromatograms data, group on the File column, and remove decoys
#decoys are duplicates of the following columns : fileName ModifiedSequence PrecursorCharge FragmentIon ProductCharge IsotopeLabelType
#let's get to it.
#
# NB NB NB NB: This code assumes that 1. Decoys come after targets. AND 2. Decoys are 1:1 with targets.  There are no targets without decoys, and no target transitions without decoy transitions.
if options.analysis_type=="lfq":
    chromatogram_df=pandas.read_csv(options.skyline_chromatograms,sep="\t",engine='python')
    #chromatogram_df=pandas.read_csv(options.skyline_chromatograms,sep="\t",low_memory=False)

    if options.filter_decoys:
        #groups=chromatogram_df.groupby(by='FileName')
        #groups=chromatogram_df.groupby(by=['FileName','PeptideModifiedSequence','PrecursorCharge','FragmentIon','ProductCharge','IsotopeLabelType'])
        groups=chromatogram_df.groupby(by=['FileName','PeptideModifiedSequence','IsotopeLabelType'])
        filtered_df=[]
        first=True
        for groupKey,group in groups:
            filtered_df.append(group.head(n=(len(group)/2)))
        filtered_dataframe=pandas.concat(filtered_df)
    else:
        filtered_dataframe=chromatogram_df
    filtered_dataframe.to_csv("filtered_chromatograms.tsv",sep="\t",index=False)
    print "Finished outputting the filtered chromatograms..."






with open(options.fasta) as fasta_file:
    sequences = []
    identifiers = []
    descriptions = []
    lengths=[]
    for seq_record in SeqIO.parse(fasta_file, 'fasta'):  # (generator)
        identifiers.append(seq_record.id)
        sequences.append(str(seq_record.seq))
        descriptions.append(seq_record.description.replace("\t",","))
        lengths.append(len(seq_record.seq))
fasta_df=pandas.DataFrame({'id':identifiers,'description':descriptions,'Sequence':sequences})

fasta_df.to_csv("fasta_df.tsv",sep="\t",index=False)

#We're going to read in the MS2 scans via Pyteomics, and then store them in a pandas df.
#Then we'll write them to a TSV file <----- New version implements HDF5 instead to keep it out of memory!
mzml_files=options.mzml.split("___")
mzml_file_scans=[]
ms2fail=False
first=True
h5file=h5py.File(options.MS2Output,'w')

#group_information    
for each_file in mzml_files:
    if first:
        first=False
    else:
        h5file.close()
        gc.collect()
        h5file=h5py.File(options.MS2Output,'a')
    #this_group=h5file.create_group(str(group_information.loc[group_information['Original File Name']==each_file.rsplit(".",1)[0],'Crux File Integer'].iloc[0]))  #This makes groups by crux integer
    this_group=h5file.create_group(each_file)
    try:
        print "Reading in ",each_file
        mzml_reader=mzml.read(each_file,iterative=True)
        for each_scan in mzml_reader:
            if each_scan['ms level']==2:
                this_dataset=this_group.create_dataset(str(each_scan['index']),compression="gzip",compression_opts=9,dtype="float32",data=numpy.column_stack((each_scan['m/z array'],each_scan['intensity array'])).T)
                this_dataset.attrs["scan_index"]=each_scan['index']
        print "done with it........\n"
        del mzml_reader
        del each_scan
        h5file.flush()
    except:
        print "Failed to read in file",each_file
        ms2fail=True
#if not ms2fail:
    #scans_df=pandas.DataFrame.from_dict(mzml_file_scans)
    #scans_df.to_csv("ms2_scans.tsv",sep="\t",index=False)
    #del scans_df

h5file.close()
#job_process.wait()
#out, err = job_process.communicate()
#print out,err
#print "Finished the job processing."


job_process.wait()
out, err = job_process.communicate()
print out,err
print "Finished the job processing."


print "Handling the renaming!" #This is done only for both the PSM_table and the MSstats outputs...

psm_table=pandas.read_csv(options.psm_table,sep='\t',engine='python')#low_memory=False


### The code below is a mixture of things that we need in LFQ and Qual analysis...
### It is meant to query uniprot, and generate name mappings for the MSstats csv table, and the PSM table.
if options.analysis_type=="lfq":
    combined_results=pandas.read_csv(options.msstats_comparison,sep=',',engine='python')#low_memory=False
    #combined_results=pandas.read_csv(options.msstats_comparison,sep=',',index_col=False)
    #print combined_results
    combined_results['backup']=combined_results['Protein']

protein_series=psm_table['protein id']
uni_mapping_dict={}
uniprot_pattern=r'([A-NR-Z][0-9][A-Z][A-Z0-9][A-Z0-9][0-9]|[OPQ][0-9][A-Z0-9][A-Z0-9][A-Z0-9][0-9])'
protein_extracted=protein_series.str.extractall(uniprot_pattern)[0].unique()#for the PSM table...
protein_full_name_series_list=protein_series.str.split(",",expand=False) #This will be used later to properly adjust positioning of AA's...


if "accToGene" in options.renameProteinType or options.renameProteinType is None: #This means we'll take the acc's and replace with GENE NAMES
    if options.analysis_type=="lfq":
        combined_results['uniprot_acc']=combined_results['Protein'].str.extract(uniprot_pattern,expand=False)
    #uni_mapping_dict_combined=uni.map(list(combined_results['uniprot_acc'].unique()),f='ACC',t='GENENAME')
    try:
        uni_mapping_dict_combined=uni.map(list(protein_extracted),f='ACC',t='GENENAME')
    except:
        uni_mapping_dict_combined={} #failure....
    for each_acc in uni_mapping_dict_combined:
        catch=uni_mapping_dict_combined[each_acc]
        if len(catch)>0:
            uni_mapping_dict[str(each_acc)]=str(catch.pop())
    second_query=[]
    for each_acc in protein_extracted:
{'fraction_non_alphanumeric': 0.16534438324321593, 'fraction_numerical': 0.007681034918389004, 'mean_word_length': 5.975678533662319, 'pattern_counts': {'":': 6, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1920959', 'n_tokens_mistral': 6683, 'n_tokens_neox': 6406, 'n_words': 1116}