annotate codon_usage.py @ 6:5e0f80e08344 draft

Uploaded
author fabio
date Tue, 11 Dec 2018 16:10:00 -0500
parents 5b61f1b564b3
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
1 #!/home/gianmarco/galaxy-python/python
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
2
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
3 import Bio
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
4 from Bio import SeqIO
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
5 from Bio.Data import CodonTable
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
6 import re
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
7 import sys
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
8 import os
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
9 import pandas as pd
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
10
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
11 def read_input(data = "example.fna"):
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
12
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
13 seqs = ""
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
14 with open(data, "rU") as handle:
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
15 for record in SeqIO.parse(handle, "fasta"):
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
16 seqs = seqs + str(record.seq)
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
17
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
18 return seqs
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
19
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
20 def codon_usage(seqs, codonTable):
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
21
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
22 codon_usage = {}
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
23 tmp = [x for x in re.split(r'(\w{3})', seqs) if x != ""]
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
24
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
25 b_cod_table = CodonTable.unambiguous_dna_by_name[codonTable].forward_table
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
26
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
27
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
28 for cod in CodonTable.unambiguous_dna_by_name[codonTable].stop_codons:
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
29 b_cod_table[cod] = "_Stop"
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
30
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
31 for cod in CodonTable.unambiguous_dna_by_name[codonTable].start_codons:
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
32 b_cod_table[cod + " Start"] = b_cod_table[cod]
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
33 b_cod_table.pop(cod)
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
34
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
35 aas = set(b_cod_table.values())
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
36
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
37
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
38 for aa in aas:
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
39 codon_usage[aa] = {}
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
40 for codon in b_cod_table.keys():
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
41 if b_cod_table[codon] == aa:
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
42 codon_usage[aa][codon] = tmp.count(codon.split(" ")[0])
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
43
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
44
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
45 tups = {(outerKey, innerKey): values for outerKey, innerDict in codon_usage.iteritems() for innerKey, values in innerDict.iteritems()}
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
46
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
47 #aas_ = set(tups.keys())
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
48
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
49 #stops_ = {el for el in aas_ if el[0] == "Stop"}
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
50 #aas_ = list(aas_.difference(stops_))
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
51 #stops_ = list(stops_)
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
52 #aas_.sort()
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
53 #stops_.sort()
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
54
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
55 codon_usage_ = pd.DataFrame(pd.Series(tups), columns = ["Count"])
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
56 codon_usage_.index = codon_usage_.index.set_names(["AA", "Codon"])
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
57 #codon_usage_.index.reindex(pd.MultiIndex.from_tuples([aas_, stops_], names=('AA', 'Codon')), level=[0,1])
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
58
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
59
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
60 codon_usage_['Proportion'] = codon_usage_.groupby(level=0).transform(lambda x: (x / x.sum()).round(2))
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
61
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
62 return {"Dictionary": codon_usage, "Tuples": tups, "Table": codon_usage_}
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
63
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
64
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
65
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
66 if __name__ == '__main__':
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
67
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
68
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
69 seqs = read_input(data=sys.argv[1])
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
70 out = codon_usage(seqs,"Bacterial")
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
71
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
72
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
73 with open(sys.argv[2], "w") as outf:
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
74 out["Table"].to_csv(outf, sep="\t")
5b61f1b564b3 Uploaded
fabio
parents:
diff changeset
75 #sys.stdout.write(out['Table'])