annotate getalleleseq.py @ 8:698ede7baba9 draft

Uploaded
author boris
date Tue, 18 Mar 2014 12:25:24 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8
698ede7baba9 Uploaded
boris
parents:
diff changeset
1 #!/usr/bin/env python
698ede7baba9 Uploaded
boris
parents:
diff changeset
2 # Boris Rebolledo-Jaramillo (boris-at-bx.psu.edu)
698ede7baba9 Uploaded
boris
parents:
diff changeset
3 #
698ede7baba9 Uploaded
boris
parents:
diff changeset
4 #usage: getalleleseq.py [-h] [-l INT] [-j FILE] [-d DIR] alleles
698ede7baba9 Uploaded
boris
parents:
diff changeset
5 #
698ede7baba9 Uploaded
boris
parents:
diff changeset
6 #Given a table with minor and major alleles per position, it generates the
698ede7baba9 Uploaded
boris
parents:
diff changeset
7 #minor and major allele sequences in FASTA format
698ede7baba9 Uploaded
boris
parents:
diff changeset
8 #
698ede7baba9 Uploaded
boris
parents:
diff changeset
9 #positional arguments:
698ede7baba9 Uploaded
boris
parents:
diff changeset
10 # alleles Table containing minor and major allele base per
698ede7baba9 Uploaded
boris
parents:
diff changeset
11 # position. cols: [id, chr, pos, A, C, G, T, cvrg,
698ede7baba9 Uploaded
boris
parents:
diff changeset
12 # plody, major, minor, freq_minor]
698ede7baba9 Uploaded
boris
parents:
diff changeset
13 #
698ede7baba9 Uploaded
boris
parents:
diff changeset
14 #optional arguments:
698ede7baba9 Uploaded
boris
parents:
diff changeset
15 # -h, --help show this help message and exit
698ede7baba9 Uploaded
boris
parents:
diff changeset
16 # -l INT, --seq-length INT
698ede7baba9 Uploaded
boris
parents:
diff changeset
17 # Background sequence length. Bases in an artifical
698ede7baba9 Uploaded
boris
parents:
diff changeset
18 # all-N-sequence of length INT will be replaced by
698ede7baba9 Uploaded
boris
parents:
diff changeset
19 # either the major or minor allele base accordingly
698ede7baba9 Uploaded
boris
parents:
diff changeset
20 # -j FILE, --major-seq FILE
698ede7baba9 Uploaded
boris
parents:
diff changeset
21 # File to write major allele sequences in FASTA multiple
698ede7baba9 Uploaded
boris
parents:
diff changeset
22 # alignment format.
698ede7baba9 Uploaded
boris
parents:
diff changeset
23 # -d DIR, --minor-dir DIR
698ede7baba9 Uploaded
boris
parents:
diff changeset
24 # Per sample minor allele sequences will be written to
698ede7baba9 Uploaded
boris
parents:
diff changeset
25 # this directory
698ede7baba9 Uploaded
boris
parents:
diff changeset
26 #
698ede7baba9 Uploaded
boris
parents:
diff changeset
27 # The expected columns in the alleles table follow Nicholas Stoler's
698ede7baba9 Uploaded
boris
parents:
diff changeset
28 # Variant Annotator tool format. See Variant Annotator in Galaxy's tool shed
698ede7baba9 Uploaded
boris
parents:
diff changeset
29 # http://testtoolshed.g2.bx.psu.edu/repos/nick/allele_counts_1 for more details
698ede7baba9 Uploaded
boris
parents:
diff changeset
30 #
698ede7baba9 Uploaded
boris
parents:
diff changeset
31 # Expected columns:
698ede7baba9 Uploaded
boris
parents:
diff changeset
32 # 1. sample_id
698ede7baba9 Uploaded
boris
parents:
diff changeset
33 # 2. chr
698ede7baba9 Uploaded
boris
parents:
diff changeset
34 # 3. position
698ede7baba9 Uploaded
boris
parents:
diff changeset
35 # 4 counts for A's
698ede7baba9 Uploaded
boris
parents:
diff changeset
36 # 5. counts for C's
698ede7baba9 Uploaded
boris
parents:
diff changeset
37 # 6. counts for G's
698ede7baba9 Uploaded
boris
parents:
diff changeset
38 # 7. counts for T's
698ede7baba9 Uploaded
boris
parents:
diff changeset
39 # (8. counts for a's)
698ede7baba9 Uploaded
boris
parents:
diff changeset
40 # (9. counts for c's)
698ede7baba9 Uploaded
boris
parents:
diff changeset
41 # (10. counts for g's)
698ede7baba9 Uploaded
boris
parents:
diff changeset
42 # (11. counts for t's)
698ede7baba9 Uploaded
boris
parents:
diff changeset
43 # 8. (12.) Coverage
698ede7baba9 Uploaded
boris
parents:
diff changeset
44 # 9. (13.) Number of alleles passing a given criteria
698ede7baba9 Uploaded
boris
parents:
diff changeset
45 # 10. (14.) Major allele
698ede7baba9 Uploaded
boris
parents:
diff changeset
46 # 11. (15.) Minor allele
698ede7baba9 Uploaded
boris
parents:
diff changeset
47 # 12. (16.) Minor allele frequency in position
698ede7baba9 Uploaded
boris
parents:
diff changeset
48
698ede7baba9 Uploaded
boris
parents:
diff changeset
49 import sys
698ede7baba9 Uploaded
boris
parents:
diff changeset
50 import os
698ede7baba9 Uploaded
boris
parents:
diff changeset
51 import argparse
698ede7baba9 Uploaded
boris
parents:
diff changeset
52
698ede7baba9 Uploaded
boris
parents:
diff changeset
53 def createseq(sample, allele, seq_size, table):
698ede7baba9 Uploaded
boris
parents:
diff changeset
54 """Generate major or minor allele sequence"""
698ede7baba9 Uploaded
boris
parents:
diff changeset
55 out_sequence = ['N' for i in range(seq_size)]
698ede7baba9 Uploaded
boris
parents:
diff changeset
56 sample_data = [line for line in table if line[0] == sample]
698ede7baba9 Uploaded
boris
parents:
diff changeset
57
698ede7baba9 Uploaded
boris
parents:
diff changeset
58 for entry in sample_data:
698ede7baba9 Uploaded
boris
parents:
diff changeset
59 position = int(entry[2])
698ede7baba9 Uploaded
boris
parents:
diff changeset
60 if len(entry)==12:
698ede7baba9 Uploaded
boris
parents:
diff changeset
61 number_of_alleles = int(entry[8])
698ede7baba9 Uploaded
boris
parents:
diff changeset
62 major_allele = entry[9].strip()
698ede7baba9 Uploaded
boris
parents:
diff changeset
63 minor_allele = entry[10].strip()
698ede7baba9 Uploaded
boris
parents:
diff changeset
64 else:
698ede7baba9 Uploaded
boris
parents:
diff changeset
65 number_of_alleles = int(entry[12])
698ede7baba9 Uploaded
boris
parents:
diff changeset
66 major_allele = entry[13].strip()
698ede7baba9 Uploaded
boris
parents:
diff changeset
67 minor_allele = entry[14].strip()
698ede7baba9 Uploaded
boris
parents:
diff changeset
68
698ede7baba9 Uploaded
boris
parents:
diff changeset
69 if allele == 'major':
698ede7baba9 Uploaded
boris
parents:
diff changeset
70 out_sequence[position-1] = major_allele
698ede7baba9 Uploaded
boris
parents:
diff changeset
71 elif allele == 'minor':
698ede7baba9 Uploaded
boris
parents:
diff changeset
72 if number_of_alleles >= 2:
698ede7baba9 Uploaded
boris
parents:
diff changeset
73 out_sequence[position-1] = minor_allele
698ede7baba9 Uploaded
boris
parents:
diff changeset
74 else:
698ede7baba9 Uploaded
boris
parents:
diff changeset
75 out_sequence[position-1] = major_allele
698ede7baba9 Uploaded
boris
parents:
diff changeset
76 return out_sequence
698ede7baba9 Uploaded
boris
parents:
diff changeset
77
698ede7baba9 Uploaded
boris
parents:
diff changeset
78 def printseq(sample,allele,seq,output):
698ede7baba9 Uploaded
boris
parents:
diff changeset
79 """Print out sequence"""
698ede7baba9 Uploaded
boris
parents:
diff changeset
80 #print >> output, '>{0}_{1}'.format(sample,allele)
698ede7baba9 Uploaded
boris
parents:
diff changeset
81 print >> output, '>{0}{1}'.format(sample,allele)
698ede7baba9 Uploaded
boris
parents:
diff changeset
82 for i in range(0,len(seq),70):
698ede7baba9 Uploaded
boris
parents:
diff changeset
83 print >> output, ''.join(seq[i:i+70])
698ede7baba9 Uploaded
boris
parents:
diff changeset
84
698ede7baba9 Uploaded
boris
parents:
diff changeset
85 def main():
698ede7baba9 Uploaded
boris
parents:
diff changeset
86 parser = argparse.ArgumentParser(description='Given a table with minor and major alleles per position, it generates the minor and major allele sequences in FASTA format', epilog='Boris Rebolledo-Jaramillo (boris-at-bx.psu.edu)')
698ede7baba9 Uploaded
boris
parents:
diff changeset
87 parser.add_argument('alleles', type=str, help='Table containing minor and major allele base per position. cols: [id, chr, pos, A, C, G, T, cvrg, plody, major, minor, freq_minor] ')
698ede7baba9 Uploaded
boris
parents:
diff changeset
88 parser.add_argument('-l','--seq-length', type=int, metavar='INT', help='Background sequence length. Bases in an artifical all-N-sequence of length INT will be replaced by either the major or minor allele base accordingly')
698ede7baba9 Uploaded
boris
parents:
diff changeset
89 parser.add_argument('-j','--major-seq', type=str, metavar='FILE', help='File to write major allele sequences in FASTA multiple alignment format.')
698ede7baba9 Uploaded
boris
parents:
diff changeset
90 parser.add_argument('-d', '--minor-dir', type=str, metavar='DIR', default='.', help="Per sample minor allele sequences will be written to this directory (Default: current directory)")
698ede7baba9 Uploaded
boris
parents:
diff changeset
91 parser.add_argument('-p', '--minor-prefix', type=str, metavar='STR', nargs='?', const='', default='', help=argparse.SUPPRESS) #Galaxy compatibility
698ede7baba9 Uploaded
boris
parents:
diff changeset
92 args = parser.parse_args()
698ede7baba9 Uploaded
boris
parents:
diff changeset
93
698ede7baba9 Uploaded
boris
parents:
diff changeset
94
698ede7baba9 Uploaded
boris
parents:
diff changeset
95 try:
698ede7baba9 Uploaded
boris
parents:
diff changeset
96 table = [line.strip().split('\t') for line in list(open(args.alleles)) if "#" not in line]
698ede7baba9 Uploaded
boris
parents:
diff changeset
97 samples = sorted(list(set([ line[0] for line in table ])))
698ede7baba9 Uploaded
boris
parents:
diff changeset
98 except:
698ede7baba9 Uploaded
boris
parents:
diff changeset
99 sys.exit('\nERROR: Could not open %s\n' % args.alleles)
698ede7baba9 Uploaded
boris
parents:
diff changeset
100 try:
698ede7baba9 Uploaded
boris
parents:
diff changeset
101 major_out = open(args.major_seq, 'w+')
698ede7baba9 Uploaded
boris
parents:
diff changeset
102 except:
698ede7baba9 Uploaded
boris
parents:
diff changeset
103 sys.exit('\nCould not create %s\n' % args.major_seq)
698ede7baba9 Uploaded
boris
parents:
diff changeset
104
698ede7baba9 Uploaded
boris
parents:
diff changeset
105 # Single file for all major allele sequences in FASTA multiple alignment
698ede7baba9 Uploaded
boris
parents:
diff changeset
106 for sample in samples:
698ede7baba9 Uploaded
boris
parents:
diff changeset
107 sequence = createseq(sample,'major',args.seq_length,table)
698ede7baba9 Uploaded
boris
parents:
diff changeset
108 #printseq(sample,'major',sequence,major_out)
698ede7baba9 Uploaded
boris
parents:
diff changeset
109 printseq(sample,'',sequence,major_out)
698ede7baba9 Uploaded
boris
parents:
diff changeset
110 major_out.close()
698ede7baba9 Uploaded
boris
parents:
diff changeset
111
698ede7baba9 Uploaded
boris
parents:
diff changeset
112 # Sample specific minor allele sequence in FASTA format
698ede7baba9 Uploaded
boris
parents:
diff changeset
113 try:
698ede7baba9 Uploaded
boris
parents:
diff changeset
114 os.makedirs(args.minor_dir)
698ede7baba9 Uploaded
boris
parents:
diff changeset
115 except:
698ede7baba9 Uploaded
boris
parents:
diff changeset
116 pass
698ede7baba9 Uploaded
boris
parents:
diff changeset
117
698ede7baba9 Uploaded
boris
parents:
diff changeset
118 for sample in samples:
698ede7baba9 Uploaded
boris
parents:
diff changeset
119 if args.minor_prefix: # to fit Galaxy requirements
698ede7baba9 Uploaded
boris
parents:
diff changeset
120 name = sample.replace('_','')
698ede7baba9 Uploaded
boris
parents:
diff changeset
121 minor_name = "%s_%s_%s" % ('primary',args.minor_prefix,name+'-minor_visible_fasta')
698ede7baba9 Uploaded
boris
parents:
diff changeset
122 else: # for non-Galaxy
698ede7baba9 Uploaded
boris
parents:
diff changeset
123 minor_name = sample+'-minor.fa'
698ede7baba9 Uploaded
boris
parents:
diff changeset
124 minor_out = open(os.path.join(args.minor_dir, minor_name), 'w+')
698ede7baba9 Uploaded
boris
parents:
diff changeset
125 sequence = createseq(sample,'minor',args.seq_length,table)
698ede7baba9 Uploaded
boris
parents:
diff changeset
126 #printseq(sample,'minor',sequence,minor_out)
698ede7baba9 Uploaded
boris
parents:
diff changeset
127 printseq(sample,'_minor',sequence,minor_out)
698ede7baba9 Uploaded
boris
parents:
diff changeset
128 minor_out.close()
698ede7baba9 Uploaded
boris
parents:
diff changeset
129
698ede7baba9 Uploaded
boris
parents:
diff changeset
130 if __name__ == "__main__": main()