annotate codon_usage.py @ 28:e85a1feaaf38 draft

Uploaded
author fabio
date Wed, 12 Dec 2018 08:12:25 -0500
parents
children 3cb2af2435d3
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
28
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
1 #!/usr/bin/env python
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
2
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
3 import Bio as Bio
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
4 from Bio import SeqIO
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
5 from Bio.Data import CodonTable
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
6 import re
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
7 from pprint import pprint
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
8 import argparse as ap
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
9 import sys
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
10 import os
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
11 import pandas as pd
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
12 #from BCBio import GFF
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
13
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
14
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
15 def read_input(data = "example.fna", type_ = "fasta"):
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
16
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
17 """
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
18 Accepted formats:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
19 - fasta (multifasta)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
20 - gff
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
21 - gbk
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
22
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
23 """
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
24
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
25
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
26 seqs = ""
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
27
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
28 if type_ == "fasta":
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
29 with open(data, "rU") as handle:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
30 for record in SeqIO.parse(handle, type_):
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
31 seqs = seqs + str(record.seq)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
32
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
33
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
34 #elif type_ == "gff":
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
35 # with open(data, "rU") as handle:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
36 # for record in GFF.parse(handle):
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
37 # seqs = seqs + str(record.seq)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
38
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
39 elif type_ == "gbk":
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
40 with open(data, "rU") as input_handle:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
41 for record in SeqIO.parse(input_handle, "genbank"):
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
42 seqs = seqs + str(record.seq)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
43
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
44
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
45 return seqs
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
46
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
47 def codon_usage(seqs, codonTable):
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
48
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
49 codon_usage = {}
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
50 tmp = [x for x in re.split(r'(\w{3})', seqs) if x != ""]
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
51
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
52 b_cod_table = CodonTable.unambiguous_dna_by_name[codonTable].forward_table
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
53
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
54
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
55 for cod in CodonTable.unambiguous_dna_by_name[codonTable].stop_codons:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
56 b_cod_table[cod] = "_Stop"
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
57
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
58 for cod in CodonTable.unambiguous_dna_by_name[codonTable].start_codons:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
59 b_cod_table[cod + " Start"] = b_cod_table[cod]
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
60 b_cod_table.pop(cod)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
61
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
62 aas = set(b_cod_table.values())
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
63
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
64
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
65 for aa in aas:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
66 codon_usage[aa] = {}
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
67 for codon in b_cod_table.keys():
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
68 if b_cod_table[codon] == aa:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
69 codon_usage[aa][codon] = tmp.count(codon.split(" ")[0])
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
70
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
71
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
72 tups = {(outerKey, innerKey): values for outerKey, innerDict in codon_usage.iteritems() for innerKey, values in innerDict.iteritems()}
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
73
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
74 codon_usage_ = pd.DataFrame(pd.Series(tups), columns = ["Count"])
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
75 codon_usage_.index = codon_usage_.index.set_names(["AA", "Codon"])
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
76 codon_usage_['Proportion'] = codon_usage_.groupby(level=0).transform(lambda x: (x / x.sum()).round(2))
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
77
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
78 return {"Dictionary": codon_usage, "Tuples": tups, "Table": codon_usage_}
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
79
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
80 if __name__ == '__main__':
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
81
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
82 parser = ap.ArgumentParser(description=
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
83 'This script takes as input gff, gbk and single or multifasta files and \n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
84 'compute the codon usage for a specified codon table.\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
85 'Usage:\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
86 'python codon_usage.py -i example.gbk -t genebank -o gbk_example -c Bacterial\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
87 'python codon_usage.py -i example.ffn -t fasta -o fasta_example -c Bacterial\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
88 'python codon_usage.py -i example.gff -t gff -o gff_example -c Bacterial\n',
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
89 formatter_class=ap.RawTextHelpFormatter)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
90
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
91 parser.add_argument('-i','--input', help='The path to the input file',required=True)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
92 parser.add_argument('-t','--type', help=
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
93 'The format of the file [genebank, fasta, gff ...]', required=True)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
94 parser.add_argument('-c','--codonTable', help=
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
95 'The codon table to be used [Standard, Bacterial, Archaeal ...]\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
96 'Alternative Flatworm Mitochondrial,\\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
97 'Alternative Yeast Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
98 'Archaeal,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
99 'Ascidian Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
100 'Bacterial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
101 'Blastocrithidia Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
102 'Blepharisma Macronuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
103 'Candidate Division SR1,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
104 'Chlorophycean Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
105 'Ciliate Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
106 'Coelenterate Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
107 'Condylostoma Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
108 'Dasycladacean Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
109 'Echinoderm Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
110 'Euplotid Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
111 'Flatworm Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
112 'Gracilibacteria,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
113 'Hexamita Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
114 'Invertebrate Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
115 'Karyorelict Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
116 'Mesodinium Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
117 'Mold Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
118 'Mycoplasma,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
119 'Pachysolen tannophilus Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
120 'Peritrich Nuclear,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
121 'Plant Plastid,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
122 'Protozoan Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
123 'Pterobranchia Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
124 'SGC0,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
125 'SGC1,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
126 'SGC2,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
127 'SGC3,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
128 'SGC4,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
129 'SGC5,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
130 'SGC8,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
131 'SGC9,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
132 'Scenedesmus obliquus Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
133 'Spiroplasma,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
134 'Standard,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
135 'Thraustochytrium Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
136 'Trematode Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
137 'Vertebrate Mitochondrial,\n'
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
138 'Yeast Mitochondrial\n', required=True)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
139
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
140 parser.add_argument('-o','--output', help='Description for bar argument', required=True)
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
141 args = vars(parser.parse_args())
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
142
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
143 seqs = read_input(data=args['input'], type_=args['type'])
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
144 out = codon_usage(seqs, args['codonTable'])
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
145
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
146 with open(args['output'], "w") as outf:
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
147 out["Table"].to_csv(outf, sep="\t", index_label=["AA", "Codon"])
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
148
e85a1feaaf38 Uploaded
fabio
parents:
diff changeset
149