P-GRe/get_seq.py at main · sebc31/P-GRe · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import sys
from re import sub

########################################################################################################################
#### FUNCTIONS

def insertNewlines(string, every=60):
    string=str(string)
    return '\n'.join(string[i:i+every] for i in range(0, len(string), every))

########################################################################################################################
#### MAIN

#ARGV
#1. P-GRe annotation
#2. Miniprot raw results

#### MAIN
######## READING P-GRe OUPUT

#Loading GFF to keep pseudogenes id

pseudogene = set()
with open(sys.argv[1]) as results:
    for line in results:
        id = line.split("=")[1].split(";")[0]
        pseudogene.add(id)

#### MAIN
######## READING RAW MINIPROT GFF

sequence_dic = {}
j = 0

with open(sys.argv[2]) as raw_mp:
    for line in raw_mp:
        if line[0:5] == "##ATA":
            sequence = line
        elif "mRNA" in line and line[0:2] != "##":
            id = line.split("=")[1].split(";")[0]
            if id in pseudogene:
                j += 1
                pseudogene.remove(id)
                sequence_dic[id] = sequence
                sequence = sub(r"[^a-zA-Z*]", "", sequence)
                print(">" + id + "\n" + insertNewlines(sequence[3:]))
            else:
                sequence = ""