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