# HG changeset patch # User bornea # Date 1458071964 14400 # Node ID 019e60bd3f7f7c7a09ee61a0dd7cd5aa95fb87a4 # Parent 945f600f34cbf14362c1ce89d56c1dd6492176aa Uploaded diff -r 945f600f34cb -r 019e60bd3f7f SAINT_preprocessing.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SAINT_preprocessing.py Tue Mar 15 15:59:24 2016 -0400 @@ -0,0 +1,277 @@ +####################################################################################### +# Python-code: SAINT pre-processing from Scaffold "Samples Report" output +# Author: Brent Kuenzi +####################################################################################### +# This program reads in a raw Scaffold "Samples Report" output and a user generated +# bait file and autoformats it into prey and interaction files for SAINTexpress +# analysis +####################################################################################### +# Copyright (C) Brent Kuenzi. +# Permission is granted to copy, distribute and/or modify this document +# under the terms of the GNU Free Documentation License, Version 1.3 +# or any later version published by the Free Software Foundation; +# with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. +# A copy of the license is included in the section entitled "GNU +# Free Documentation License". +####################################################################################### +## REQUIRED INPUT ## + +# 1) infile: Scaffold "Samples Report" output +# 2) baitfile: SAINT formatted bait file generated in Galaxy +# 3) fasta_db: fasta database for use (defaults to SwissProt_HUMAN_2014_08.fasta) +# 4) prey: Y or N for generating a prey file +# 5) make_bait: String of bait names, assignment, and test or control boolean +####################################################################################### + +import sys +import os.path + + +infile = sys.argv[1] +#Scaffold "Samples Report" output. +prey = sys.argv[2] +# Y or N boolean from Galaxy. +fasta_db = sys.argv[3] +tool_path = sys.argv[8] +if fasta_db == "None": + fasta_db = str(tool_path) + "/SwissProt_HUMAN_2014_08.fasta" +make_bait = sys.argv[6] +bait_bool = sys.argv[9] + + +def bait_create(baits, infile): + # Verifies the Baits are valid in the Scaffold file and writes the Bait.txt. + baits = make_bait.split() + i = 0 + bait_file_tmp = open("bait.txt", "w") + order = [] + bait_cache = [] + while i < len(baits): + if baits[i+2] == "true": + T_C = "C" + else: + T_C = "T" + bait_line = baits[i] + "\t" + baits[i+1] + "\t" + T_C + "\n" + read_infile = open(infile, "r") + for input_line in read_infile: + input_line = input_line.strip() + temp = input_line.split('\t') + if "Quantitative Variance" in str(temp): + if baits[i] in temp: + number_bait = temp.index(str(baits[i])) + number_bait = number_bait - 9 + bait_cache.append((number_bait, str(bait_line))) + # Locates the Bait names in the column names and then sets the Baits in the + # correct order in the cache thus the - 9 because the baits start at the 9th + # column. + else: + print "Error: bad bait " + str(baits[i]) + sys.exit() + else: + pass + i = i + 3 + + bait_cache.sort() + for cache_line in bait_cache: + bait_file_tmp.write(cache_line[1]) + + bait_file_tmp.close() + +if bait_bool == 'false': + bait_create(make_bait, infile) + baitfile = "bait.txt" +else: + bait_temp_file = open(sys.argv[10], 'r') + bait_cache = bait_temp_file.readlines() + bait_file_tmp = open("bait.txt", "wr") + for cache_line in bait_cache: + bait_file_tmp.write(cache_line) + bait_file_tmp.close() + baitfile = "bait.txt" + + +class ReturnValue1(object): + def __init__(self, sequence, gene): + self.seqlength = sequence + self.genename = gene + + +class ReturnValue2(object): + def __init__(self, getdata, getproteins, getheader): + self.data = getdata + self.proteins = getproteins + self.header = getheader + + +def main(Scaffold_input, baits): + bait_check(baitfile, Scaffold_input) + make_inter(Scaffold_input) + if prey == 'true': + make_prey(Scaffold_input) + no_error_inter(Scaffold_input) + os.rename('prey.txt', sys.argv[5]) + elif prey == 'false': + if os.path.isfile('error proteins.txt') == True: + no_error_inter(Scaffold_input) + pass + elif prey != 'true' or 'false': + sys.exit("Invalid Prey Argument: Y or N") + + +def get_info(uniprot_accession_in): + # Get aminoacid lengths and gene name. + error = open('error proteins.txt', 'a+') + data = open(fasta_db, 'r') + data_lines = data.readlines() + db_len = len(data_lines) + seqlength = 0 + count = 0 + for data_line in data_lines: + if ">sp" in data_line: + namer = data_line.split("|")[2] + if uniprot_accession_in == data_line.split("|")[1]: + match = count + 1 + if 'GN=' in data_line: + lst = data_linedata_line.split('GN=') + lst2 = lst[1].split(' ') + genename = lst2[0] + if 'GN=' not in data_line: + genename = 'NA' + while ">sp" not in data_lines[match]: + if match <= db_len: + seqlength = seqlength + len(data_lines[match].strip()) + match = match + 1 + else: + break + return ReturnValue1(seqlength, genename) + elif uniprot_accession_in == namer.split(" ")[0]: + match = count + 1 + # Ensures consistent spacing throughout. + if 'GN=' in data_line: + lst = data_line.split('GN=') + lst2 = lst[1].split(' ') + genename = lst2[0] + if 'GN=' not in data_line: + genename = 'NA' + while ">sp" not in data_lines[match]: + if match <= db_len: + seqlength = seqlength + len(data_lines[match].strip()) + match = match + 1 + else: + break + return ReturnValue1(seqlength, genename) + count = count + 1 + if seqlength == 0: + error.write(uniprot_accession_in + '\t' + "Uniprot not in Fasta" + '\n') + error.close + seqlength = 'NA' + genename = 'NA' + return ReturnValue1(seqlength, genename) + + +def readtab(infile): + with open(infile, 'r') as input_file: + # read in tab-delim text + output = [] + for input_line in input_file: + input_line = input_line.strip() + temp = input_line.split('\t') + output.append(temp) + return output + + +def read_Scaffold(Scaffold_input): + # Get data, proteins and header from Scaffold output + dupes = readtab(Scaffold_input) + cnt = 0 + for Scaffold_line in dupes: + cnt += 1 + if Scaffold_line[0] == '#': + # Finds the start of second header. + header_start = cnt-1 + header = dupes[header_start] + prot_start = header.index("Accession Number") + data = dupes[header_start+1:len(dupes)-2] + # Cut off blank line and END OF FILE. + proteins = [] + for Scaffold_line in data: + Scaffold_line[4] = Scaffold_line[4].split()[0] + # Removes the (+##) that sometimes is attached. + for protein in data: + proteins.append(protein[prot_start]) + return ReturnValue2(data, proteins, header) + + +def make_inter(Scaffold_input): + bait = readtab(baitfile) + data = read_Scaffold(Scaffold_input).data + header = read_Scaffold(Scaffold_input).header + proteins = read_Scaffold(Scaffold_input).proteins + bait_index = [] + for bait_line in bait: + bait_index.append(header.index(bait_line[0])) + # Find just the baits defined in bait file. + with open('inter.txt', 'w') as inter_file: + a = 0; l = 0 + for bb in bait: + for lst in data: + inter_file.write(header[bait_index[l]] + '\t' + bb[1] + '\t' + proteins[a] + '\t' + + lst[bait_index[l]] + '\n') + a += 1 + if a == len(proteins): + a = 0; l += 1 + + +def make_prey(Scaffold_input): + proteins = read_Scaffold(Scaffold_input).proteins + output_file = open("prey.txt", 'w') + for protein in proteins: + protein = protein.replace("\n", "") + # Remove \n for input into function. + protein = protein.replace("\r", "") + # Ditto for \r. + seq = get_info(protein).seqlength + GN = get_info(protein).genename + if seq != 'NA': + output_file.write(protein + "\t" + str(seq) + "\t" + str(GN) + "\n") + output_file.close() + + +def no_error_inter(Scaffold_input): + # Remake inter file without protein errors from Uniprot. + err = readtab("error proteins.txt") + bait = readtab(baitfile) + data = read_Scaffold(Scaffold_input).data + header = read_Scaffold(Scaffold_input).header + bait_index = [] + for bait_line in bait: + bait_index.append(header.index(bait_line[0])) + proteins = read_Scaffold(Scaffold_input).proteins + errors = [] + for e in err: + errors.append(e[0]) + with open('inter.txt', 'w') as y: + l = 0; a = 0 + for bb in bait: + for lst in data: + if proteins[a] not in errors: + y.write(header[bait_index[l]] + '\t' + bb[1] + '\t' + proteins[a] + '\t' + + lst[bait_index[l]] + '\n') + a += 1 + if a == len(proteins): + l += 1; a = 0 + + +def bait_check(bait, Scaffold_input): + # Check that bait names share Scaffold header titles. + bait_in = readtab(bait) + header = read_Scaffold(Scaffold_input).header + for i in bait_in: + if i[0] not in header: + sys.exit("Bait must share header titles with Scaffold output") + +if __name__ == '__main__': + main(infile, baitfile) + +os.rename("inter.txt", sys.argv[4]) +os.rename("bait.txt", sys.argv[7])