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