# HG changeset patch
# User pravs
# Date 1513935192 18000
# Node ID 032f1ee171edcf5403a624c8490d09aa544df788
planemo upload
diff -r 000000000000 -r 032f1ee171ed pep_pointer.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/pep_pointer.py Fri Dec 22 04:33:12 2017 -0500
@@ -0,0 +1,221 @@
+
+#
+# Author: Praveen Kumar
+# Updated: Nov 8th, 2017
+#
+#
+#
+
+import re
+
+
+def main():
+ import sys
+ if len(sys.argv) == 4:
+ inputFile = sys.argv
+ infh = open(inputFile[1], "r")
+ # infh = open("Mus_musculus.GRCm38.90.chr.gtf", "r")
+
+ gtf = {}
+ gtf_transcript = {}
+ gtf_gene = {}
+ for each in infh.readlines():
+ a = each.split("\t")
+ if re.search("^[^#]", each):
+ if re.search("gene_biotype \"protein_coding\"", a[8]) and int(a[4].strip()) != int(a[3].strip()):
+ type = a[2].strip()
+ if type == "gene" or type == "exon" or type == "CDS" or type == "five_prime_utr" or type == "three_prime_utr":
+ chr = "chr" + a[0].strip()
+ strand = a[6].strip()
+ if strand == "+":
+ start = a[3].strip()
+ end = a[4].strip()
+ elif strand == "-":
+ if int(a[4].strip()) > int(a[3].strip()):
+ start = a[3].strip()
+ end = a[4].strip()
+ elif int(a[4].strip()) < int(a[3].strip()):
+ start = a[4].strip()
+ end = a[3].strip()
+ else:
+ print "Something fishy in start end coordinates"
+ else:
+ print "Something fishy in reading"
+ if not gtf.has_key(strand):
+ gtf[strand] = {}
+ if not gtf[strand].has_key(type):
+ gtf[strand][type] = []
+ b = re.search("gene_id \"(.+?)\";", a[8].strip())
+ gene = b.group(1)
+ if type == "gene":
+ transcript = ""
+ else:
+ b = re.search("transcript_id \"(.+?)\";", a[8].strip())
+ transcript = b.group(1)
+ data = (chr, start, end, gene, transcript, strand, type)
+ gtf[strand][type].append(data)
+
+ if type == "exon":
+ if gtf_transcript.has_key(chr+"#"+strand):
+ if gtf_transcript[chr+"#"+strand].has_key(transcript+"#"+gene):
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene][0].append(int(start))
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene][1].append(int(end))
+ else:
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene] = [[],[]]
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene][0].append(int(start))
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene][1].append(int(end))
+ else:
+ gtf_transcript[chr+"#"+strand] = {}
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene] = [[],[]]
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene][0].append(int(start))
+ gtf_transcript[chr+"#"+strand][transcript+"#"+gene][1].append(int(end))
+
+ if type == "gene":
+ if gtf_gene.has_key(chr+"#"+strand):
+ gtf_gene[chr+"#"+strand][0].append(int(start))
+ gtf_gene[chr+"#"+strand][1].append(int(end))
+ gtf_gene[chr+"#"+strand][2].append(gene)
+ else:
+ gtf_gene[chr+"#"+strand] = [[0],[0],["no_gene"]]
+ gtf_gene[chr+"#"+strand][0].append(int(start))
+ gtf_gene[chr+"#"+strand][1].append(int(end))
+ gtf_gene[chr+"#"+strand][2].append(gene)
+
+
+
+ # "Starting Reading Intron . . ."
+
+ gtf["+"]["intron"] = []
+ gtf["-"]["intron"] = []
+ for chr_strand in gtf_transcript.keys():
+ chr = chr_strand.split("#")[0]
+ strand = chr_strand.split("#")[1]
+
+ for transcript_gene in gtf_transcript[chr_strand].keys():
+ start_list = gtf_transcript[chr_strand][transcript_gene][0]
+ end_list = gtf_transcript[chr_strand][transcript_gene][1]
+ sorted_start_index = [i[0] for i in sorted(enumerate(start_list), key=lambda x:x[1])]
+ sorted_end_index = [i[0] for i in sorted(enumerate(end_list), key=lambda x:x[1])]
+ if sorted_start_index == sorted_end_index:
+ sorted_start = sorted(start_list)
+ sorted_end = [end_list[i] for i in sorted_start_index]
+ for x in range(len(sorted_start))[1:]:
+ intron_start = sorted_end[x-1]+1
+ intron_end = sorted_start[x]-1
+ transcript = transcript_gene.split("#")[0]
+ gene = transcript_gene.split("#")[1]
+ data = (chr, str(intron_start), str(intron_end), gene, transcript, strand, "intron")
+ gtf[strand]["intron"].append(data)
+
+
+ # "Starting Reading Intergenic . . ."
+
+ gtf["+"]["intergenic"] = []
+ gtf["-"]["intergenic"] = []
+ for chr_strand in gtf_gene.keys():
+ chr = chr_strand.split("#")[0]
+ strand = chr_strand.split("#")[1]
+ start_list = gtf_gene[chr_strand][0]
+ end_list = gtf_gene[chr_strand][1]
+ gene_list = gtf_gene[chr_strand][2]
+ sorted_start_index = [i[0] for i in sorted(enumerate(start_list), key=lambda x:x[1])]
+ sorted_end_index = [i[0] for i in sorted(enumerate(end_list), key=lambda x:x[1])]
+
+ sorted_start = sorted(start_list)
+ sorted_end = [end_list[i] for i in sorted_start_index]
+ sorted_gene = [gene_list[i] for i in sorted_start_index]
+ for x in range(len(sorted_start))[1:]:
+ intergene_start = sorted_end[x-1]+1
+ intergene_end = sorted_start[x]-1
+ if intergene_start < intergene_end:
+ intergene_1 = sorted_gene[x-1]
+ intergene_2 = sorted_gene[x]
+ gene = intergene_1 + "-#-" + intergene_2
+ data = (chr, str(intergene_start), str(intergene_end), gene, "", strand, "intergenic")
+ gtf[strand]["intergenic"].append(data)
+
+ import sqlite3
+ # conn = sqlite3.connect('gtf_database.db')
+ conn = sqlite3.connect(":memory:")
+ c = conn.cursor()
+ # c.execute("DROP TABLE IF EXISTS gtf_data;")
+ # c.execute("CREATE TABLE IF NOT EXISTS gtf_data(chr text, start int, end int, gene text, transcript text, strand text, type text)")
+ c.execute("CREATE TABLE gtf_data(chr text, start int, end int, gene text, transcript text, strand text, type text)")
+
+ for strand in gtf.keys():
+ if strand == "+":
+ st = "positive"
+ elif strand == "-":
+ st = "negative"
+ else:
+ print "Something fishy in writing . . ."
+
+ for type in gtf[strand].keys():
+ data = gtf[strand][type]
+ c.executemany('INSERT INTO gtf_data VALUES (?,?,?,?,?,?,?)', data)
+
+ conn.commit()
+
+ infh = open(inputFile[2], "r")
+ # infh = open("Mouse_Data_All_peptides_withNewDBs.txt", "r")
+ data = infh.readlines()
+ # output file
+ outfh = open(inputFile[3], 'w')
+ # outfh = open("classified_1_Mouse_Data_All_peptides_withNewDBs.txt", "w")
+
+ for each in data:
+ a = each.split("\t")
+ chr = a[0].strip()
+ pep_start = a[1].strip()
+ pep_end = a[2].strip()
+ strand = a[5].strip()
+ c.execute("select * from gtf_data where type = 'CDS' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\tCDS\n")
+ else:
+ c.execute("select * from gtf_data where type = 'five_prime_utr' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\tfive_prime_utr\n")
+ else:
+ c.execute("select * from gtf_data where type = 'three_prime_utr' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\tthree_prime_utr\n")
+ else:
+ c.execute("select * from gtf_data where type = 'exon' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\texon\n")
+ else:
+ c.execute("select * from gtf_data where type = 'intron' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\tintron\n")
+ else:
+ c.execute("select * from gtf_data where type = 'gene' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\tgene\n")
+ else:
+ c.execute("select * from gtf_data where type = 'intergenic' and chr = '"+chr+"' and start <= "+pep_start+" and end >= "+pep_end+" and strand = '"+strand+"' ")
+ rows = c.fetchall()
+ if len(rows) > 0:
+ outfh.write(each.strip() + "\tintergene\n")
+ else:
+ outfh.write(each.strip() + "\tOVERLAPPING_ON_TWO_REGIONS: PLEASE_LOOK_MANUALLY (Will be updated in next version)\n")
+
+ conn.close()
+ outfh.close()
+ else:
+ print "USAGE: python pep_pointer.py "
+ return None
+
+if __name__ == "__main__":
+ main()
+
+
+
+
+
diff -r 000000000000 -r 032f1ee171ed pep_pointer.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/pep_pointer.xml Fri Dec 22 04:33:12 2017 -0500
@@ -0,0 +1,36 @@
+
+
+ python
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+@misc{peppointer,
+ author={Kumar, Praveen},
+ year={2017},
+ title={PepPointer}
+}
+
+
+