annotate snpeff_to_peptides.py @ 1:80eff5812b2a default tip

Check for both classic and sequence_ontology name for EFF: "NON_SYNONYMOUS_CODING" or "missense_variant"
author Jim Johnson <jj@umn.edu>
date Thu, 06 Nov 2014 13:57:31 -0600
parents fcb7188fa0d2
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
1 #!/usr/bin/env python
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
2 """
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
3 #
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
4 #------------------------------------------------------------------------------
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
5 # University of Minnesota
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
6 # Copyright 2013, Regents of the University of Minnesota
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
7 #------------------------------------------------------------------------------
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
8 # Author:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
9 #
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
10 # James E Johnson
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
11 #
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
12 #------------------------------------------------------------------------------
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
13 """
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
14
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
15
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
16 """
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
17 This tool takes a SnpEff VCF file and an Ensembl pep.all.fa file ( e.g. Homo_sapiens.GRCh37.73.pep.all.fa )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
18 It outputs a peptide fasta file with the variant peptide sequence that result from NON_SYNONYMOUS_CODING effects
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
19
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
20 """
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
21
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
22 import sys,re,os.path
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
23 import tempfile
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
24 import optparse
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
25 from optparse import OptionParser
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
26 import logging
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
27
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
28 ## dictionary for Amino Acid Abbreviations
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
29 aa_abbrev_dict = dict()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
30 aa_abbrev_dict['Phe'] = 'F'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
31 aa_abbrev_dict['Leu'] = 'L'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
32 aa_abbrev_dict['Ser'] = 'S'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
33 aa_abbrev_dict['Tyr'] = 'Y'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
34 aa_abbrev_dict['Cys'] = 'C'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
35 aa_abbrev_dict['Trp'] = 'W'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
36 aa_abbrev_dict['Pro'] = 'P'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
37 aa_abbrev_dict['His'] = 'H'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
38 aa_abbrev_dict['Gln'] = 'Q'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
39 aa_abbrev_dict['Arg'] = 'R'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
40 aa_abbrev_dict['Ile'] = 'I'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
41 aa_abbrev_dict['Met'] = 'M'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
42 aa_abbrev_dict['Thr'] = 'T'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
43 aa_abbrev_dict['Asn'] = 'N'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
44 aa_abbrev_dict['Lys'] = 'K'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
45 aa_abbrev_dict['Val'] = 'V'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
46 aa_abbrev_dict['Ala'] = 'A'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
47 aa_abbrev_dict['Asp'] = 'D'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
48 aa_abbrev_dict['Glu'] = 'E'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
49 aa_abbrev_dict['Gly'] = 'G'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
50
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
51 ## Get the peptide ID and sequence a given ID
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
52 def get_sequence(id,seq_file):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
53 fh = open(seq_file, 'r')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
54 try:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
55 for (ln,line) in enumerate(fh):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
56 if line.find(id) >= 0:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
57 fields = line.split('\t')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
58 return ( ' '.join(fields[0:-1]),fields[-1].rstrip() if fields and len(fields) > 0 else None )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
59 except Exception, e:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
60 print >> sys.stderr, "failed: %s" % e
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
61 finally:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
62 fh.close()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
63
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
64 def fasta_to_tabular(fasta_file,tabular_file):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
65 inFile = open(fasta_file,'r')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
66 outFile = open(tabular_file,'w')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
67 for i, line in enumerate( inFile ):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
68 line = line.rstrip( '\r\n' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
69 if not line or line.startswith( '#' ):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
70 continue
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
71 if line.startswith( '>' ):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
72 #Don't want any existing tabs to trigger extra columns:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
73 line = line.replace('\t', ' ')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
74 if i > 0:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
75 outFile.write('\n')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
76 outFile.write(line[1:])
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
77 outFile.write('\t')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
78 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
79 outFile.write(line)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
80 if i > 0:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
81 outFile.write('\n')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
82 if inFile:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
83 inFile.close()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
84 if outFile:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
85 outFile.close()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
86
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
87 def __main__():
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
88 #Parse Command Line
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
89 parser = optparse.OptionParser()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
90 parser.add_option( '-i', '--input', dest='input', help='The input snpeff vcf file with HGVS annotations (else read from stdin)' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
91 parser.add_option( '-o', '--output', dest='output', help='The output fasta (else write to stdout)' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
92 parser.add_option( '-p', '--protein_fasta', dest='protein_fasta', default=None, help='The Esembl protein fasta in tabular format' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
93 parser.add_option( '-l', '--leading_aa_num', dest='leading_aa_num', type='int', default=None, help='leading number of AAs to output' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
94 parser.add_option( '-t', '--trailing_aa_num', dest='trailing_aa_num', type='int', default=None, help='trailing number of AAs to output' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
95 parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stdout' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
96 (options, args) = parser.parse_args()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
97
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
98 # need protein_fasta file
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
99 fastaFile = options.protein_fasta
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
100 if options.protein_fasta == None:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
101 print >> sys.stderr, "Ensembl protein_fasta tabular file required"
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
102 exit(4)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
103 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
104 # determine if fasta is already in tabular format
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
105 is_tabular = False
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
106 standard_aa = '^[AC-IK-WY]+$'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
107 standard_na = '^[ACGTN]+$'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
108 inFile = open(fastaFile,'r')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
109 try:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
110 nseq = 0
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
111 for i, line in enumerate( inFile ):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
112 line = line.rstrip( '\r\n' )
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
113 if not line or line.startswith( '#' ):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
114 continue
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
115 fields = line.split('\t')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
116 if len(fields) < 2:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
117 is_tabular = False
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
118 if line[0] != '>':
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
119 print >> sys.stderr, "failed: %s does not appear to be a fasta file" % fastaFile
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
120 exit(4)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
121 break
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
122 if re.match('^[A-Z]+$',fields[-1].upper()):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
123 is_tabular = True
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
124 nseq += 1
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
125 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
126 if line[0] != '>':
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
127 print >> sys.stderr, "failed: %s does not appear to be a fasta file" % fastaFile
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
128 exit(4)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
129 if nseq > 10:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
130 break
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
131 finally:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
132 if inFile:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
133 inFile.close()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
134 if not is_tabular:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
135 fastaFile = tempfile.NamedTemporaryFile(prefix='pep_fasta_',suffix=".tab",dir=os.getcwd()).name
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
136 fasta_to_tabular(options.protein_fasta,fastaFile)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
137 # vcf input
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
138 if options.input != None:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
139 try:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
140 inputPath = os.path.abspath(options.input)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
141 inputFile = open(inputPath, 'r')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
142 except Exception, e:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
143 print >> sys.stderr, "failed: %s" % e
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
144 exit(2)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
145 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
146 inputFile = sys.stdin
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
147 # output
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
148 if options.output != None:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
149 try:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
150 outputPath = os.path.abspath(options.output)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
151 outputFile = open(outputPath, 'w')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
152 except Exception, e:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
153 print >> sys.stderr, "failed: %s" % e
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
154 exit(3)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
155 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
156 outputFile = sys.stdout
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
157 ## Amino_Acid_Change notations
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
158 # G528R
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
159 # p.Gly528Arg/c.1582G>C
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
160 aa_change_regex = '([A-Z])(\d+)([A-Z])' # G528R
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
161 aa_hgvs_regex = 'p\.([A-Z][a-z][a-z])(\d+)([A-Z][a-z][a-z])(/c\.(\d+)([ACGTN])>([ACGTN]))' # p.Gly528Arg/c.1582G>C
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
162 # Save VCF file header, not currently used
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
163 vcf_header = []
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
164 reading_entries = False
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
165 try:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
166 for linenum,line in enumerate(inputFile):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
167 ## print >> sys.stderr, "%d: %s\n" % (linenum,line)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
168 if line.startswith('##'):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
169 vcf_header.append(line)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
170 # May need to check SnpEff version in the header, the EFF info changed between versions 2 and 3
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
171 ##SnpEffVersion
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
172 elif line.startswith('#CHROM'):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
173 reading_entries = True
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
174 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
175 fields = line.split('\t')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
176 # This is the current format of the EFF entry:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
177 # EFF=missense(MODERATE|MISSENSE|Ggg/Cgg|G528R|802|SCNN1D|protein_coding|CODING|ENST00000379116|12|1);OICR=(ENST00000379116|1808)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
178 # If this becomes variable, will need to dynamically pattern this on the defintion in the vcf header:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
179 ##INFO=<ID=EFF,Number=.,Type=String,Description="Predicted effects for this variant.Format: 'Effect ( Effect_Impact | Functional_Class | Codon_Change | Amino_Acid_Change| Amino_Acid_length | Gene_Name | Transcript_BioType | Gene_Coding | Transcript_ID | Exon_Rank | Genotype_Number [ | ERRORS | WARNINGS ] )' ">
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
180 (chrom,pos,id,ref,alts,qual,filter,info) = fields[0:8]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
181 for info_item in info.split(';'):
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
182 try:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
183 if info_item.find('=') < 0:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
184 continue
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
185 (key,val) = info_item.split('=',1)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
186 if key == 'EFF':
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
187 effects = val.split(',')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
188 for effect in effects:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
189 (eff,effs) = effect.rstrip(')').split('(')
1
80eff5812b2a Check for both classic and sequence_ontology name for EFF: "NON_SYNONYMOUS_CODING" or "missense_variant"
Jim Johnson <jj@umn.edu>
parents: 0
diff changeset
190 if not (eff == 'NON_SYNONYMOUS_CODING' or eff == 'missense_variant'):
0
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
191 continue
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
192 eff_fields = effs.split('|')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
193 (impact,functional_class,codon_change,aa_change,aa_len,gene_name,biotype,coding,transcript,exon) = eff_fields[0:10]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
194 if transcript:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
195 aa_pos = None # 1-based position
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
196 alt_aa = '_'
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
197 # parse aa_change
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
198 # get AA change position and alternate Animo Acid
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
199 sap = aa_change
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
200 m = re.match(aa_change_regex,aa_change)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
201 if m:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
202 aa_pos = int(m.groups()[1])
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
203 alt_aa = m.groups()[2]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
204 else:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
205 m = re.match(aa_hgvs_regex,aa_change)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
206 if m:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
207 aa_pos = int(m.groups()[1])
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
208 ref_aa = aa_abbrev_dict[m.groups()[0]]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
209 alt_aa = aa_abbrev_dict[m.groups()[2]]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
210 sap = "%s%d%s" % (ref_aa,aa_pos,alt_aa)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
211 if not aa_pos:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
212 continue
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
213 # get AA sequence
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
214 aa_offset = aa_pos - 1
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
215 (pep_id,pep_seq) = get_sequence(transcript,fastaFile)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
216 if not pep_seq:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
217 continue
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
218 start_pos = max(aa_offset - options.leading_aa_num, 0) if options.leading_aa_num else 0
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
219 end_pos = min(aa_offset + options.trailing_aa_num + 1, len(pep_seq)) if options.trailing_aa_num else len(pep_seq)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
220 # transform sequence
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
221 alt_seq = pep_seq[start_pos:aa_offset] + alt_aa + pep_seq[aa_offset+1:end_pos]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
222 # >ENSP00000363782 pep:known chromosome:GRCh37:1:22778472:22853855:1 gene:ENSG00000184677 transcript:ENST00000374651 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:1:22778472 codon_change:Gtg/Atg sap:V885M
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
223 pep_id = re.sub('pep:[a-z]*','pep:sap',pep_id)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
224 hdr = ">%s snp_location:%s:%s codon_change:%s sap:%s\n" % (pep_id, chrom, pos, codon_change, sap)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
225 outputFile.write(hdr)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
226 if options.debug:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
227 trimmed_seq = pep_seq[start_pos:end_pos]
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
228 outputFile.write(trimmed_seq)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
229 outputFile.write('\n')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
230 outputFile.write(alt_seq)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
231 outputFile.write('\n')
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
232 except Exception, e:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
233 print >> sys.stderr, "failed: %s" % e
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
234 except Exception, e:
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
235 print >> sys.stderr, "failed: %s" % e
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
236 exit(1)
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
237
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
238 if __name__ == "__main__" : __main__()
fcb7188fa0d2 Uploaded
jjohnson
parents:
diff changeset
239