|
5
|
1
|
|
1
|
2 """
|
|
|
3 Python-code: Merge Scaffold Samples Report files
|
|
|
4 @author = Brent Kuenzi
|
|
|
5 @email = Brent.Kuenzi@moffitt.org
|
|
|
6 """
|
|
|
7 #######################################################################################
|
|
|
8 ## Description: ##
|
|
|
9 # This program will filter a fasta file using a data file
|
|
|
10 # containing an "Accession" or "Accession Number" column
|
|
|
11 ## Required input: ##
|
|
|
12 # 1) fasta file to be filtered
|
|
|
13 # 2) data file containing a "Accession" or "Accession Number" column
|
|
|
14
|
|
|
15
|
|
|
16 import sys
|
|
|
17 import os
|
|
|
18 def readtab(infile): # read in tab-delim text
|
|
|
19 with open(infile,'r') as x:
|
|
|
20 output = []
|
|
|
21 for line in x:
|
|
|
22 line = line.strip()
|
|
|
23 temp = line.split('\t')
|
|
|
24 output.append(temp)
|
|
|
25 return output
|
|
|
26 def getAccessions(infile): # get list of protein accessions from your data
|
|
|
27 data = readtab(infile)
|
|
|
28 cnt = 0
|
|
|
29 header_start = 0
|
|
|
30 prot_start = 0
|
|
|
31 for i in data:
|
|
10
|
32 if "Accession Number" in i: # finds the start of header
|
|
1
|
33 header_start = cnt
|
|
|
34 break
|
|
|
35 cnt += 1
|
|
|
36 header = data[header_start]
|
|
14
|
37 for i in header:
|
|
|
38 if i == "Accession":
|
|
|
39 prot_start = header.index("Accession")
|
|
|
40 if i == "Accession Number":
|
|
|
41 prot_start = header.index("Accession Number")
|
|
|
42 if i == "Main Accession":
|
|
|
43 prot_start = header.index("Main Accession")
|
|
1
|
44 proteins = []
|
|
10
|
45 for protein in data[header_start:]:
|
|
|
46 if len(protein) > prot_start:
|
|
|
47 proteins.append(protein[prot_start])
|
|
1
|
48 return proteins
|
|
|
49 def FilterFastaSeq(infile,accession): # fasta file and UniprotID/SwissprotID
|
|
|
50 input_data = readtab(infile)
|
|
|
51 seq=[]
|
|
|
52 header=[]
|
|
|
53 temp=[]
|
|
|
54 flag = 0
|
|
|
55 cnt = 0
|
|
|
56 for i in input_data:
|
|
|
57 cnt+=1
|
|
|
58 if flag == 1: # once we have a hit, start adding the sequences
|
|
|
59 if ">" not in i[0]: # don't add the headers to the sequence
|
|
|
60 temp.append(i[0])
|
|
|
61 if i[0].startswith(">"): # is it a fasta header?
|
|
|
62 if temp != []: # if it is a continued fasta header, add old sequences to the sequence list
|
|
|
63 # will this cutoff the last on of the file?
|
|
|
64 merged = "\n".join(temp)
|
|
|
65 if merged!="":
|
|
|
66 seq.append(merged)
|
|
|
67 temp=[]
|
|
|
68 lst = i[0].split('|')
|
|
|
69 ID1 = lst[1]
|
|
|
70 lst2 = lst[2].split(' ')
|
|
|
71 ID2 = lst2[0]
|
|
|
72 flag = 0
|
|
|
73 if ID1 in accession: # are the IDs part of your data?
|
|
|
74 flag+=1
|
|
|
75 header.append(i[0]) # if it is then append it
|
|
|
76 if ID2 in accession: # are the IDs part of your data?
|
|
|
77 flag+=1
|
|
|
78 header.append(i[0]) # if it is then append it
|
|
|
79 if cnt == len(input_data): # account for last fasta sequence in file
|
|
|
80 if temp != []:
|
|
|
81 merged = "\n".join(temp)
|
|
|
82 if merged!="":
|
|
|
83 seq.append(merged)
|
|
|
84 cnt=0
|
|
9
|
85 x = open("output.txt","w")
|
|
7
|
86 for i in header:
|
|
|
87 x.write(i+'\n'+seq[cnt]+'\n')
|
|
|
88 cnt+=1
|
|
|
89 x.close()
|
|
11
|
90 fasta = sys.argv[1] # fasta file to filter
|
|
|
91 data = sys.argv[2] # scaffold report #2 -- filename
|
|
1
|
92
|
|
|
93 FilterFastaSeq(fasta,getAccessions(data))
|
|
12
|
94 os.rename("output.txt", sys.argv[3]) |
