Mercurial > repos > galaxyp > peptide_genomic_coordinate

comparison peptide_genomic_coordinate.py @ 0:5f49ffce52cb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit be7e9677908b7864ef0b965a1e219a1840eeb2ec
author galaxyp
date Wed, 03 Apr 2019 04:04:18 -0400
parents
children cb0378d2d487
comparison
equal deleted inserted replaced
-1:000000000000 0:5f49ffce52cb
1 #!/usr/bin/env python
2 #
3 # Author: Praveen Kumar
4 # University of Minnesota
5 #
6 # Get peptide's genomic coordinate from the protein's genomic mapping sqlite file (which is derived from the https://toolshed.g2.bx.psu.edu/view/galaxyp/translate_bed/038ecf54cbec)
7 #
8 # python peptideGenomicCoordinate.py <peptide_list> <mz_to_sqlite DB> <genomic mapping file DB> <output.bed>
9 #
10 import sys
11 import sqlite3
12
13
14 def main():
15 conn = sqlite3.connect(sys.argv[2])
16 c = conn.cursor()
17 c.execute("DROP table if exists novel")
18 conn.commit()
19 c.execute("CREATE TABLE novel(peptide text)")
20 pepfile = open(sys.argv[1],"r")
21
22 pep_seq = []
23 for seq in pepfile.readlines():
24 seq = seq.strip()
25 pep_seq.append(tuple([seq]))
26
27 c.executemany("insert into novel(peptide) values(?)", pep_seq)
28 conn.commit()
29
30 c.execute("SELECT distinct psm.sequence, ps.id, ps.sequence from db_sequence ps, psm_entries psm, novel n, proteins_by_peptide pbp where psm.sequence = n.peptide and pbp.peptide_ref = psm.id and pbp.id = ps.id")
31 rows = c.fetchall()
32
33 conn1 = sqlite3.connect(sys.argv[3])
34 c1 = conn1.cursor()
35
36 outfh = open(sys.argv[4], "w")
37
38 master_dict = {}
39 for each in rows:
40 peptide = each[0]
41 acc = each[1]
42 acc_seq = each[2]
43
44 c1.execute("SELECT chrom,start,end,name,strand,cds_start,cds_end FROM feature_cds_map map WHERE map.name = '"+acc+"'")
45 coordinates = c1.fetchall()
46
47 if len(coordinates) != 0:
48 pep_start = 0
49 pep_end = 0
50 flag = 0
51 splice_flag = 0
52 spliced_peptide = []
53 for each_entry in coordinates:
54 chromosome = each_entry[0]
55 start = int(each_entry[1])
56 end = int(each_entry[2])
57 strand = each_entry[4]
58 cds_start = int(each_entry[5])
59 cds_end = int(each_entry[6])
60 pep_pos_start = (acc_seq.find(peptide)*3)
61 pep_pos_end = pep_pos_start + (len(peptide)*3)
62 if pep_pos_start >= cds_start and pep_pos_end <= cds_end:
63 if strand == "+":
64 pep_start = start + pep_pos_start - cds_start
65 pep_end = start + pep_pos_end - cds_start
66 pep_thick_start = 0
67 pep_thick_end = len(peptide)
68 flag == 1
69 else:
70 pep_end = end - pep_pos_start + cds_start
71 pep_start = end - pep_pos_end + cds_start
72 pep_thick_start = 0
73 pep_thick_end = len(peptide)
74 flag == 1
75 spliced_peptide = []
76 splice_flag = 0
77 else:
78 if flag == 0:
79 if strand == "+":
80 if pep_pos_start >= cds_start and pep_pos_start <= cds_end and pep_pos_end > cds_end:
81 pep_start = start + pep_pos_start - cds_start
82 pep_end = end
83 pep_thick_start = 0
84 pep_thick_end = (pep_end-pep_start)
85 spliced_peptide.append([pep_start,pep_end,pep_thick_start,pep_thick_end])
86 splice_flag = splice_flag + 1
87 if splice_flag == 2:
88 flag = 1
89 elif pep_pos_end >= cds_start and pep_pos_end <= cds_end and pep_pos_start < cds_start:
90 pep_start = start
91 pep_end = start + pep_pos_end - cds_start
92 pep_thick_start = (len(peptide)*3)-(pep_end-pep_start)
93 pep_thick_end = (len(peptide)*3)
94 spliced_peptide.append([pep_start,pep_end,pep_thick_start,pep_thick_end])
95 splice_flag = splice_flag + 1
96 if splice_flag == 2:
97 flag = 1
98 else:
99 pass
100 else:
101 if pep_pos_start >= cds_start and pep_pos_start <= cds_end and pep_pos_end >= cds_end:
102 pep_start = start
103 pep_end = end - pep_pos_start - cds_start
104 pep_thick_start = 0
105 pep_thick_end = (pep_end-pep_start)
106 spliced_peptide.append([pep_start,pep_end,pep_thick_start,pep_thick_end])
107 splice_flag = splice_flag + 1
108 if splice_flag == 2:
109 flag = 1
110 elif pep_pos_end >= cds_start and pep_pos_end <= cds_end and pep_pos_start <= cds_start:
111 pep_start = end - pep_pos_end + cds_start
112 pep_end = end
113 pep_thick_start = (len(peptide)*3)-(pep_end-pep_start)
114 pep_thick_end = (len(peptide)*3)
115 spliced_peptide.append([pep_start,pep_end,pep_thick_start,pep_thick_end])
116 splice_flag = splice_flag + 1
117 if splice_flag == 2:
118 flag = 1
119 else:
120 pass
121
122 if len(spliced_peptide) == 0:
123 if strand == "+":
124 bed_line = [chromosome, str(pep_start), str(pep_end), peptide, "255", strand, str(pep_start), str(pep_end), "0", "1", str(pep_end-pep_start), "0"]
125 else:
126 bed_line = [chromosome, str(pep_start), str(pep_end), peptide, "255", strand, str(pep_start), str(pep_end), "0", "1", str(pep_end-pep_start), "0"]
127 outfh.write("\t".join(bed_line)+"\n")
128 else:
129 if strand == "+":
130 pep_entry = spliced_peptide
131 pep_start = min([pep_entry[0][0], pep_entry[1][0]])
132 pep_end = max([pep_entry[0][1], pep_entry[1][1]])
133 blockSize = [str(min([pep_entry[0][3], pep_entry[1][3]])),str(max([pep_entry[0][3], pep_entry[1][3]])-min([pep_entry[0][3], pep_entry[1][3]]))]
134 blockStarts = ["0", str(pep_end-pep_start-(max([pep_entry[0][3], pep_entry[1][3]])-min([pep_entry[0][3], pep_entry[1][3]])))]
135 bed_line = [chromosome, str(pep_start), str(pep_end), peptide, "255", strand, str(pep_start), str(pep_end), "0", "2", ",".join(blockSize), ",".join(blockStarts)]
136 outfh.write("\t".join(bed_line)+"\n")
137 else:
138 pep_entry = spliced_peptide
139 pep_start = min([pep_entry[0][0], pep_entry[1][0]])
140 pep_end = max([pep_entry[0][1], pep_entry[1][1]])
141 blockSize = [str(min([pep_entry[0][3], pep_entry[1][3]])),str(max([pep_entry[0][3], pep_entry[1][3]])-min([pep_entry[0][3], pep_entry[1][3]]))]
142 blockStarts = ["0", str(pep_end-pep_start-(max([pep_entry[0][3], pep_entry[1][3]])-min([pep_entry[0][3], pep_entry[1][3]])))]
143 bed_line = [chromosome, str(pep_start), str(pep_end), peptide, "255", strand, str(pep_start), str(pep_end), "0", "2", ",".join(blockSize), ",".join(blockStarts)]
144 outfh.write("\t".join(bed_line)+"\n")
145 c.execute("DROP table novel")
146 conn.commit()
147 conn.close()
148 conn1.close()
149 outfh.close()
150 pepfile.close()
151
152 return None
153 if __name__ == "__main__":
154 main()