my_src/TagGenRefSeq_control.py at master · geparada/my_src · 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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
import sys
import csv


RNA = []
Genoma = []
SeqName = []
Resultados=[]


def SJtag(AliSeq, len_tag, GenomicTags=""):
    '''Generates tags of the splice-juntions
    *Input file = USCS BLAT alignments table (22 columns ) joined with fasta converted to tab (3 columns). This file has 25 columns.
    *Usage = python Tag-Fulinfo inputfile tag length
    *Output = seq ID , intron name, strand, tag sequence, tag length '''


    reader1 = csv.reader(open(AliSeq), dialect='excel-tab' )

         # row[0]  matches
         # row[1]  misMatches
         # row[2]  repMaches
         # row[3]  nCount
         # row[4]  qNumInsert
         # row[5]  qBaseInsert
         # row[6]  tNumInsert
         # row[7]  tBaseInsert
         # row[8]  strand
         # row[9] qName
         # row[10] qSize
         # row[11] qStart
         # row[12] qEnd
         # row[13] tName
         # row[14] tSize
         # row[15] tStart
         # row[16] tEnd
         # row[17] blockCount
         # row[18] blockSizes
         # row[19] qStarts
         # row[20] tStarts

         # row[21] qName
         # row[22] 1
         # row[23] sequence


    for row in reader1:

        qstarts = map (int, row[19].split(",")[1:-1])                      # obteniendo lista int de las distintas variables
        tstarts = map(int, row[20].split(",")[0:-1])
        blocksizes = map(int, row[18].split(",")[0:-1])

        n = len_tag/2                                                      #preparando argumento para su uso


        for t, b, t2, b2, q in zip(tstarts, blocksizes, tstarts[1:], blocksizes[1:]  , qstarts):

            if b>=0 and b2>=0:                                   #filtro por tamano de exones
                tbsum = t+b                                      # generando primera coordenada del intron al cual se le va a hacer el tag
                start = q-n                                     #generando intervalos para hacer el slicing del la secuencia
                end = q+n


                if row[8] == '+' :                          #Para los que aliniean en la hebra +


                    if start<0:                             #Precausiones generar buenos tag del primer y ultimo tag
                        start = 0
                    if end>len(row[23]):
                        end=len((row[23]))
                    tag = row[23][q:end] + row[23][start:q]


                if row[8] == '-' :                                #Para los que alinian en la hebra - es todo al inverso

                    if end>len(row[23]):
                        end=len(row[23])
                    tag = row[23][-q:-start] + row[23][-end:-q]
                    if start<=0:
                        tag = row[23][-end:]


                if GenomicTags!="tags":
                    name_seq_tag = row[9], row[13] + str(tbsum) + row[8] + str(t2), tag
                    SeqName.append(name_seq_tag)
                if GenomicTags=="tags":
                    #print row[13] + str(tbsum) + row[8] + str(t2) + row[9], row[13] + str(tbsum) + row[8] + str(t2) + tag
		    #print row[13] + str(tbsum) + row[8] + str(t2), tag, len(tag)
		    print ">" + row[9] + row[13] + str(tbsum) + row[8] + str(t2) + "\n" + tag

def switch (file2=""):
    if file2=="tags":
        pass
    else:
        CompareTags(sys.argv[3], sys.argv[4])


def CompareTags(GalaxyGenomico, output="F"):

    reader2 = csv.reader(open(GalaxyGenomico), dialect='excel-tab' )
        # row [0] intronName
        # row [9] seq

    for intro in reader2:
        B = intro[0], intro[9]
        Genoma.append(B)

    set_Genoma = set(Genoma)
    result = [x for x in SeqName if (x[1], x[2]) in set_Genoma]
    for c, d, e in result:
        if output=="F":
            print c, d, e
        if output=="D":
            filtrados = c, d, e
            Resultados.append(filtrados)

    if output=="D":
        sacados = set(SeqName)-set(Resultados)
        for c, d, e in sacados:
            print c, d, e


if __name__ == '__main__':                                      #Permite tomar argumentos mediante el terminal
    SJtag(sys.argv[1], int(sys.argv[2]), sys.argv[3])
    switch(sys.argv[3])

#SJtag("BLAToincDNAmm9",  30, "BLATIntronescDNAmm9")
#switch("BLATIntronescDNAmm9")
#CompareTags("BLATIntronescDNAmm9")