Mercurial > repos > jjohnson > snpeff_to_peptides
changeset 0:41a666a3d8a5 draft default tip
Uploaded
| author | jjohnson |
|---|---|
| date | Tue, 17 Dec 2013 18:45:13 -0500 |
| parents | |
| children | |
| files | snpeff_to_peptides.py snpeff_to_peptides.xml test-data/all_pep.fa test-data/all_pep.tabular test-data/peptides_10_10.fa test-data/snpeff.vcf |
| diffstat | 6 files changed, 413 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/snpeff_to_peptides.py Tue Dec 17 18:45:13 2013 -0500 @@ -0,0 +1,239 @@ +#!/usr/bin/env python +""" +# +#------------------------------------------------------------------------------ +# University of Minnesota +# Copyright 2013, Regents of the University of Minnesota +#------------------------------------------------------------------------------ +# Author: +# +# James E Johnson +# +#------------------------------------------------------------------------------ +""" + + +""" +This tool takes a SnpEff VCF file and an Ensembl pep.all.fa file ( e.g. Homo_sapiens.GRCh37.73.pep.all.fa ) +It outputs a peptide fasta file with the variant peptide sequence that result from NON_SYNONYMOUS_CODING effects + +""" + +import sys,re,os.path +import tempfile +import optparse +from optparse import OptionParser +import logging + +## dictionary for Amino Acid Abbreviations +aa_abbrev_dict = dict() +aa_abbrev_dict['Phe'] = 'F' +aa_abbrev_dict['Leu'] = 'L' +aa_abbrev_dict['Ser'] = 'S' +aa_abbrev_dict['Tyr'] = 'Y' +aa_abbrev_dict['Cys'] = 'C' +aa_abbrev_dict['Trp'] = 'W' +aa_abbrev_dict['Pro'] = 'P' +aa_abbrev_dict['His'] = 'H' +aa_abbrev_dict['Gln'] = 'Q' +aa_abbrev_dict['Arg'] = 'R' +aa_abbrev_dict['Ile'] = 'I' +aa_abbrev_dict['Met'] = 'M' +aa_abbrev_dict['Thr'] = 'T' +aa_abbrev_dict['Asn'] = 'N' +aa_abbrev_dict['Lys'] = 'K' +aa_abbrev_dict['Val'] = 'V' +aa_abbrev_dict['Ala'] = 'A' +aa_abbrev_dict['Asp'] = 'D' +aa_abbrev_dict['Glu'] = 'E' +aa_abbrev_dict['Gly'] = 'G' + +## Get the peptide ID and sequence a given ID +def get_sequence(id,seq_file): + fh = open(seq_file, 'r') + try: + for (ln,line) in enumerate(fh): + if line.find(id) >= 0: + fields = line.split('\t') + return ( ' '.join(fields[0:-1]),fields[-1].rstrip() if fields and len(fields) > 0 else None ) + except Exception, e: + print >> sys.stderr, "failed: %s" % e + finally: + fh.close() + +def fasta_to_tabular(fasta_file,tabular_file): + inFile = open(fasta_file,'r') + outFile = open(tabular_file,'w') + for i, line in enumerate( inFile ): + line = line.rstrip( '\r\n' ) + if not line or line.startswith( '#' ): + continue + if line.startswith( '>' ): + #Don't want any existing tabs to trigger extra columns: + line = line.replace('\t', ' ') + if i > 0: + outFile.write('\n') + outFile.write(line[1:]) + outFile.write('\t') + else: + outFile.write(line) + if i > 0: + outFile.write('\n') + if inFile: + inFile.close() + if outFile: + outFile.close() + +def __main__(): + #Parse Command Line + parser = optparse.OptionParser() + parser.add_option( '-i', '--input', dest='input', help='The input snpeff vcf file with HGVS annotations (else read from stdin)' ) + parser.add_option( '-o', '--output', dest='output', help='The output fasta (else write to stdout)' ) + parser.add_option( '-p', '--protein_fasta', dest='protein_fasta', default=None, help='The Esembl protein fasta in tabular format' ) + parser.add_option( '-l', '--leading_aa_num', dest='leading_aa_num', type='int', default=None, help='leading number of AAs to output' ) + parser.add_option( '-t', '--trailing_aa_num', dest='trailing_aa_num', type='int', default=None, help='trailing number of AAs to output' ) + parser.add_option( '-d', '--debug', dest='debug', action='store_true', default=False, help='Turn on wrapper debugging to stdout' ) + (options, args) = parser.parse_args() + + # need protein_fasta file + fastaFile = options.protein_fasta + if options.protein_fasta == None: + print >> sys.stderr, "Ensembl protein_fasta tabular file required" + exit(4) + else: + # determine if fasta is already in tabular format + is_tabular = False + standard_aa = '^[AC-IK-WY]+$' + standard_na = '^[ACGTN]+$' + inFile = open(fastaFile,'r') + try: + nseq = 0 + for i, line in enumerate( inFile ): + line = line.rstrip( '\r\n' ) + if not line or line.startswith( '#' ): + continue + fields = line.split('\t') + if len(fields) < 2: + is_tabular = False + if line[0] != '>': + print >> sys.stderr, "failed: %s does not appear to be a fasta file" % fastaFile + exit(4) + break + if re.match('^[A-Z]+$',fields[-1].upper()): + is_tabular = True + nseq += 1 + else: + if line[0] != '>': + print >> sys.stderr, "failed: %s does not appear to be a fasta file" % fastaFile + exit(4) + if nseq > 10: + break + finally: + if inFile: + inFile.close() + if not is_tabular: + fastaFile = tempfile.NamedTemporaryFile(prefix='pep_fasta_',suffix=".tab",dir=os.getcwd()).name + fasta_to_tabular(options.protein_fasta,fastaFile) + # vcf input + if options.input != None: + try: + inputPath = os.path.abspath(options.input) + inputFile = open(inputPath, 'r') + except Exception, e: + print >> sys.stderr, "failed: %s" % e + exit(2) + else: + inputFile = sys.stdin + # output + if options.output != None: + try: + outputPath = os.path.abspath(options.output) + outputFile = open(outputPath, 'w') + except Exception, e: + print >> sys.stderr, "failed: %s" % e + exit(3) + else: + outputFile = sys.stdout + ## Amino_Acid_Change notations + # G528R + # p.Gly528Arg/c.1582G>C + aa_change_regex = '([A-Z])(\d+)([A-Z])' # G528R + 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 + # Save VCF file header, not currently used + vcf_header = [] + reading_entries = False + try: + for linenum,line in enumerate(inputFile): + ## print >> sys.stderr, "%d: %s\n" % (linenum,line) + if line.startswith('##'): + vcf_header.append(line) + # May need to check SnpEff version in the header, the EFF info changed between versions 2 and 3 + ##SnpEffVersion + elif line.startswith('#CHROM'): + reading_entries = True + else: + fields = line.split('\t') + # This is the current format of the EFF entry: + # EFF=missense(MODERATE|MISSENSE|Ggg/Cgg|G528R|802|SCNN1D|protein_coding|CODING|ENST00000379116|12|1);OICR=(ENST00000379116|1808) + # If this becomes variable, will need to dynamically pattern this on the defintion in the vcf header: + ##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 ] )' "> + (chrom,pos,id,ref,alts,qual,filter,info) = fields[0:8] + for info_item in info.split(';'): + try: + if info_item.find('=') < 0: + continue + (key,val) = info_item.split('=',1) + if key == 'EFF': + effects = val.split(',') + for effect in effects: + (eff,effs) = effect.rstrip(')').split('(') + if not eff == 'NON_SYNONYMOUS_CODING': + continue + eff_fields = effs.split('|') + (impact,functional_class,codon_change,aa_change,aa_len,gene_name,biotype,coding,transcript,exon) = eff_fields[0:10] + if transcript: + aa_pos = None # 1-based position + alt_aa = '_' + # parse aa_change + # get AA change position and alternate Animo Acid + sap = aa_change + m = re.match(aa_change_regex,aa_change) + if m: + aa_pos = int(m.groups()[1]) + alt_aa = m.groups()[2] + else: + m = re.match(aa_hgvs_regex,aa_change) + if m: + aa_pos = int(m.groups()[1]) + ref_aa = aa_abbrev_dict[m.groups()[0]] + alt_aa = aa_abbrev_dict[m.groups()[2]] + sap = "%s%d%s" % (ref_aa,aa_pos,alt_aa) + if not aa_pos: + continue + # get AA sequence + aa_offset = aa_pos - 1 + (pep_id,pep_seq) = get_sequence(transcript,fastaFile) + if not pep_seq: + continue + start_pos = max(aa_offset - options.leading_aa_num, 0) if options.leading_aa_num else 0 + end_pos = min(aa_offset + options.trailing_aa_num + 1, len(pep_seq)) if options.trailing_aa_num else len(pep_seq) + # transform sequence + alt_seq = pep_seq[start_pos:aa_offset] + alt_aa + pep_seq[aa_offset+1:end_pos] + # >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 + pep_id = re.sub('pep:[a-z]*','pep:sap',pep_id) + hdr = ">%s snp_location:%s:%s codon_change:%s sap:%s\n" % (pep_id, chrom, pos, codon_change, sap) + outputFile.write(hdr) + if options.debug: + trimmed_seq = pep_seq[start_pos:end_pos] + outputFile.write(trimmed_seq) + outputFile.write('\n') + outputFile.write(alt_seq) + outputFile.write('\n') + except Exception, e: + print >> sys.stderr, "failed: %s" % e + except Exception, e: + print >> sys.stderr, "failed: %s" % e + exit(1) + +if __name__ == "__main__" : __main__() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/snpeff_to_peptides.xml Tue Dec 17 18:45:13 2013 -0500 @@ -0,0 +1,76 @@ +<?xml version="1.0"?> +<tool id="snpeff_to_peptides" name="SnpEff to Peptide fasta" version="0.0.1"> + <description> to create a Search DB fasta for variant SAP peptides</description> + <command interpreter="python">snpeff_to_peptides.py --input "$snpeff_vcf" --protein_fasta "$all_pep_fasta" --output "$peptide_variant_fasta" + #if $leading_aa_num: + --leading_aa_num $leading_aa_num + #end if + #if $trailing_aa_num: + --trailing_aa_num $trailing_aa_num + #end if + </command> + <inputs> + <param name="snpeff_vcf" type="data" format="vcf" label="SnpEff generated VCF file with NON_SYNONYMOUS_CODING annotations"/> + <param name="all_pep_fasta" type="data" format="fasta,tabular" label="Ensembl all_pep.fa" + help="An Ensembl all_pep.fa file corresponding to the genome build used for SnpEff (May be converted to tabular fasta format)"/> + <param name="leading_aa_num" type="integer" value="30" min="0" optional="true" label="Preceeding AAs" + help="The number of Amino Acids to include before the variant position (leave blank to include all)"/> + <param name="trailing_aa_num" type="integer" value="30" min="0" optional="true" label="Following AAs" + help="The number of Amino Acids to include after the variant position (leave blank to include all)"/> + </inputs> + <stdio> + <exit_code range="1:" level="fatal" description="Error" /> + </stdio> + <outputs> + <data name="peptide_variant_fasta" metadata_source="all_pep_fasta" format="fasta"/> + </outputs> + <tests> + <test> + <param name="snpeff_vcf" value="snpeff.vcf" ftype="vcf" dbkey="hg19"/> + <param name="all_pep_fasta" value="all_pep.fa" ftype="fasta" dbkey="hg19"/> + <param name="leading_aa_num" value="10"/> + <param name="trailing_aa_num" value="10"/> + <output name="peptide_variant_fasta" file="peptides_10_10.fa"/> + </test> + <test> + <param name="snpeff_vcf" value="snpeff.vcf" ftype="vcf" dbkey="hg19"/> + <param name="all_pep_fasta" value="all_pep.tabular" ftype="tabular" dbkey="hg19"/> + <param name="leading_aa_num" value="10"/> + <param name="trailing_aa_num" value="10"/> + <output name="peptide_variant_fasta" file="peptides_10_10.fa"/> + </test> + </tests> + <help> +**SnpEff to Peptide Fasta** + +This generates a fasta file of peptide sequences with SAPs ( Single Amino acid Polymorphisms ) +from the NON_SYNONYMOUS_CODING EFF annnotations from the SnpEff_ application. +The SnpEff VCF may be filtered or annotated using SnpSift. + +The following is appended to the fasta ID line: snp_location:chr:position codon_change:nnn/nnn sap:AposA + +For VCF entry:: + + chr1 22846709 . G A 9.31 . DP=2;VDB=0.0174;AF1=1;AC1=2;DP4=0,0,1,1;MQ=20;FQ=-33;EFF=NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|Gtg/Atg|V885M|1127|ZBTB40|protein_coding|CODING|ENST00000374651|12|1) PL 40,6,0 + +The peptide fasta entry that matches transcript ID: ENST00000374651 would be:: + + >ENSP00000363782 pep:known chromosome:GRCh37:1:22778472:22853855:1 gene:ENSG00000184677 transcript:ENST00000374651 gene_biotype:protein_coding transcript_biotype:protein_coding + +The ID of the output peptide fasta ID would be:: + + >ENSP00000363782 pep:sap chromosome:GRCh37:1:22778472:22853855:1 gene:ENSG00000184677 transcript:ENST00000374651 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:22846709 codon_change:Gtg/Atg sap:V885M + + +.. _SnpEff: http://snpeff.sourceforge.net/index.html + +**Citation** + +SnpEff citation: +"A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3.", Cingolani P, Platts A, Wang le L, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM. Fly (Austin). 2012 Apr-Jun;6(2):80-92. PMID: 22728672 [PubMed - in process] + +SnpSift citation: +"Using Drosophila melanogaster as a model for genotoxic chemical mutational studies with a new program, SnpSift", Cingolani, P., et. al., Frontiers in Genetics, 3, 2012. + + </help> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/all_pep.fa Tue Dec 17 18:45:13 2013 -0500 @@ -0,0 +1,47 @@ +>ENSP00000363782 pep:known chromosome:GRCh37:1:22778472:22853855:1 gene:ENSG00000184677 transcript:ENST00000374651 gene_biotype:protein_coding transcript_biotype:protein_coding +MELPNYSRQLLQQLYTLCKEQQFCDCTISIGTIYFRAHKLVLAAASLLFKTLLDNTDTISIDASVVSPEE +FALLLEMMYTGKLPVGKHNFSKIISLADSLQMFDVAVSCKNLLTSLVNCSVQGQVVRDVSAPSSETFRKE +PEKPQVEILSSEGAGEPHSSPELAATPGGPVKAETEEAAHSVSQEMSVNSPTAQESQRNAETPAETPTTA +EACSPSPAVQTFSEAKKTSTEPGCERKHYQLNFLLENEGVFSDALMVTQDVLKKLEMCSEIKGPQKEVIL +NCCEGRTPKETIENLLHRMTEEKTLTAEGLVKLLQAVKTTFPNLGLLLEKLQKSATLPSTTVQPSPDDYG +TELLRRYHENLSEIFTDNQILLKMISHMTSLAPGEREVMEKLVKRDSGSGGFNSLISAVLEKQTLSATAI +WQLLLVVQETKTCPLDLLMEEIRREPGADAFFRAVTTPEHATLETILRHNQLILEAIQQKIEYKLFTSEE +EHLAETVKEILSIPSETASPEASLRAVLSRAMEKSVPAIEICHLLCSVHKSFPGLQPVMQELAYIGVLTK +EDGEKETWKVSNKFHLEANNKEDEKAAKEDSQPGEQNDQGETGSLPGQQEKEASASPDPAKKSFICKACD +KSFHFYCRLKVHMKRCRVAKSKQVQCKECSETKDSKKELDKHQLEAHGAGGEPDAPKKKKKRLPVTCDLC +GREFAHASGMQYHKLTEHFDEKPFSCEECGAKFAANSTLKNHLRLHTGDRPFMCKHCLMTFTQASALAYH +TKKKHSEGKMYACQYCDAVFAQSIELSRHVRTHTGDKPYVCRDCGKGFRQANGLSIHLHTFHNIEDPYDC +KKCRMSFPTLQDHRKHIHEVHSKEYHPCPTCGKIFSAPSMLERHVVTHVGGKPFSCGICNKAYQQLSGLW +YHNRTHHPDVFAAQNHRSSKFSSLQCSSCDKTFPNTIEHKKHIKAEHADMKFHECDQCKELFPTPALLQV +HVKCQHSGSQPFRCLYCAATFRFPGALQHHVTTEHFKQSETTFPCELCGELFTSQAQLDSHLESEHPKVM +STETQAAASQMAQVIQTPEPVAPTEQVITLEETQLAGSQVFVTLPDSQASQASSELVAVTVEDLLDGTVT +LICGEAK +>ENSP00000346634 pep:known chromosome:GRCh37:1:36690017:36770958:1 gene:ENSG00000054118 transcript:ENST00000354618 gene_biotype:protein_coding transcript_biotype:protein_coding +MSKTNKSKSGSRSSRSRSASRSRSRSFSKSRSRSRSLSRSRKRRLSSRSRSRSYSPAHNRERNHPRVYQN +RDFRGHNRGYRRPYYFRGRNRGFYPWGQYNRGGYGNYRSNWQNYRQAYSPRRGRSRSRSPKRRSPSPRSR +SHSRNSDKSSSDRSRRSSSSRSSSNHSRVESSKRKSAKEKKSSSKDSRPSQAAGDNQGDEAKEQTFSGGT +SQDTKASESSKPWPDATYGTGSASRASAVSELSPRERSPALKSPLQSVVVRRRSPRPSPVPKPSPPLSST +SQMGSTLPSGAGYQSGTHQGQFDHGSGSLSPSKKSPVGKSPPSTGSTYGSSQKEESAASGGAAYTKRYLE +EQKTENGKDKEQKQTNTDKEKIKEKGSFSDTGLGDGKMKSDSFAPKTDSEKPFRGSQSPKRYKLRDDFEK +KMADFHKEEMDDQDKDKAKGRKESEFDDEPKFMSKVIGANKNQEEEKSGKWEGLVYAPPGKEKQRKTEEL +EEESFPERSKKEDRGKRSEGGHRGFVPEKNFRVTAYKAVQEKSSSPPPRKTSESRDKLGAKGDFPTGKSS +FSITREAQVNVRMDSFDEDLARPSGLLAQERKLCRDLVHSNKKEQEFRSIFQHIQSAQSQRSPSELFAQH +IVTIVHHVKEHHFGSSGMTLHERFTKYLKRGTEQEAAKNKKSPEIHRRIDISPSTFRKHGLAHDEMKSPR +EPGYKAEGKYKDDPVDLRLDIERRKKHKERDLKRGKSRESVDSRDSSHSRERSAEKTEKTHKGSKKQKKH +RRARDRSRSSSSSSQSSHSYKAEEYTEETEEREESTTGFDKSRLGTKDFVGPSERGGGRARGTFQFRARG +RGWGRGNYSGNNNNNSNNDFQKRNREEEWDPEYTPKSKKYYLHDDREGEGSDKWVSRGRGRGAFPRGRGR +FMFRKSSTSPKWAHDKFSGEEGEIEDDESGTENREEKDNIQPTTE +>ENSP00000318415 pep:known chromosome:GRCh37:1:89445139:89458455:-1 gene:ENSG00000213516 transcript:ENST00000321792 gene_biotype:protein_coding transcript_biotype:protein_coding +MVEADRPGKLFIGGLNTETNEKALETVFGKYGRIVEVLLIKDRETNKSRGFAFVTFESPADAKDAARDMN +GKSLDGKAIKVEQATKPSFERGRHGPPPPPRSRGPPRGFGAGRGGSGGTRGPPSRGGHMDDGGYSMNFNM +SSSRGPLPVKRGPPPRSGGPSPKRSAPSGLVRSSSGMGGRAPLSRGRDSYGGPPRREPLPSRRDVYLSPR +DDGYSTKDSYSSRDYPSSRDTRDYAPPPRDYTYRDYGHSSSRDDYPSRGYGDRDGYGRDRDYSDHPSGGS +YRDSYESYGNSRSAPLTRGPPPSYGGSSRYDDYSSSRDGYGGSRDSYSSSRSDLYSSCDRVGRQERGLPP +SVERGYPSSRDSYSSSSRGAPRGAGPGGSRSDRGGGRSRY +>ENSP00000446099 pep:known chromosome:GRCh37:1:89445142:89458643:-1 gene:ENSG00000213516 transcript:ENST00000399794 gene_biotype:protein_coding transcript_biotype:protein_coding +MVEADRPGKLFIGGLNTETNEKALETVFGKYGRIVEVLLIKDRETNKSRGFAFVTFESPADAKDAARDMN +GKSLDGKAIKVEQATKPSFERGRHGPPPPPRSRGPPRGFGAGRGGSGGTRGPPSRGGHMDDGGYSMNFNM +SSSRGPLPVKRGPPPRSGGPSPKRSAPSGLVRSSSGMGGRAPLSRGRDSYGGPPRREPLPSRRDVYLSPR +DDGYSTKDSYSSRDYPSSRDTRDYAPPPRDYTYRDYGHSSSRDDYPSRGYGDRDGYGRDRDYSDHPSGGS +YRDSYESYGNSRSAPLTRGPPPSYGGSSRYDDYSSSRDGYGGSRDSYSSSRSDLYSSCDRVGRQERGLPP +SVERGYPSSRDSYSSSSRGAPRGAGPGGSRSDRGGGRSRY
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/all_pep.tabular Tue Dec 17 18:45:13 2013 -0500 @@ -0,0 +1,4 @@ +ENSP00000363782 pep:known chromosome:GRCh37:1:22778472:22853855:1 gene:ENSG00000184677 transcript:ENST00000374651 gene_biotype:protein_coding transcript_biotype:protein_coding MELPNYSRQLLQQLYTLCKEQQFCDCTISIGTIYFRAHKLVLAAASLLFKTLLDNTDTISIDASVVSPEEFALLLEMMYTGKLPVGKHNFSKIISLADSLQMFDVAVSCKNLLTSLVNCSVQGQVVRDVSAPSSETFRKEPEKPQVEILSSEGAGEPHSSPELAATPGGPVKAETEEAAHSVSQEMSVNSPTAQESQRNAETPAETPTTAEACSPSPAVQTFSEAKKTSTEPGCERKHYQLNFLLENEGVFSDALMVTQDVLKKLEMCSEIKGPQKEVILNCCEGRTPKETIENLLHRMTEEKTLTAEGLVKLLQAVKTTFPNLGLLLEKLQKSATLPSTTVQPSPDDYGTELLRRYHENLSEIFTDNQILLKMISHMTSLAPGEREVMEKLVKRDSGSGGFNSLISAVLEKQTLSATAIWQLLLVVQETKTCPLDLLMEEIRREPGADAFFRAVTTPEHATLETILRHNQLILEAIQQKIEYKLFTSEEEHLAETVKEILSIPSETASPEASLRAVLSRAMEKSVPAIEICHLLCSVHKSFPGLQPVMQELAYIGVLTKEDGEKETWKVSNKFHLEANNKEDEKAAKEDSQPGEQNDQGETGSLPGQQEKEASASPDPAKKSFICKACDKSFHFYCRLKVHMKRCRVAKSKQVQCKECSETKDSKKELDKHQLEAHGAGGEPDAPKKKKKRLPVTCDLCGREFAHASGMQYHKLTEHFDEKPFSCEECGAKFAANSTLKNHLRLHTGDRPFMCKHCLMTFTQASALAYHTKKKHSEGKMYACQYCDAVFAQSIELSRHVRTHTGDKPYVCRDCGKGFRQANGLSIHLHTFHNIEDPYDCKKCRMSFPTLQDHRKHIHEVHSKEYHPCPTCGKIFSAPSMLERHVVTHVGGKPFSCGICNKAYQQLSGLWYHNRTHHPDVFAAQNHRSSKFSSLQCSSCDKTFPNTIEHKKHIKAEHADMKFHECDQCKELFPTPALLQVHVKCQHSGSQPFRCLYCAATFRFPGALQHHVTTEHFKQSETTFPCELCGELFTSQAQLDSHLESEHPKVMSTETQAAASQMAQVIQTPEPVAPTEQVITLEETQLAGSQVFVTLPDSQASQASSELVAVTVEDLLDGTVTLICGEAK +ENSP00000346634 pep:known chromosome:GRCh37:1:36690017:36770958:1 gene:ENSG00000054118 transcript:ENST00000354618 gene_biotype:protein_coding transcript_biotype:protein_coding MSKTNKSKSGSRSSRSRSASRSRSRSFSKSRSRSRSLSRSRKRRLSSRSRSRSYSPAHNRERNHPRVYQNRDFRGHNRGYRRPYYFRGRNRGFYPWGQYNRGGYGNYRSNWQNYRQAYSPRRGRSRSRSPKRRSPSPRSRSHSRNSDKSSSDRSRRSSSSRSSSNHSRVESSKRKSAKEKKSSSKDSRPSQAAGDNQGDEAKEQTFSGGTSQDTKASESSKPWPDATYGTGSASRASAVSELSPRERSPALKSPLQSVVVRRRSPRPSPVPKPSPPLSSTSQMGSTLPSGAGYQSGTHQGQFDHGSGSLSPSKKSPVGKSPPSTGSTYGSSQKEESAASGGAAYTKRYLEEQKTENGKDKEQKQTNTDKEKIKEKGSFSDTGLGDGKMKSDSFAPKTDSEKPFRGSQSPKRYKLRDDFEKKMADFHKEEMDDQDKDKAKGRKESEFDDEPKFMSKVIGANKNQEEEKSGKWEGLVYAPPGKEKQRKTEELEEESFPERSKKEDRGKRSEGGHRGFVPEKNFRVTAYKAVQEKSSSPPPRKTSESRDKLGAKGDFPTGKSSFSITREAQVNVRMDSFDEDLARPSGLLAQERKLCRDLVHSNKKEQEFRSIFQHIQSAQSQRSPSELFAQHIVTIVHHVKEHHFGSSGMTLHERFTKYLKRGTEQEAAKNKKSPEIHRRIDISPSTFRKHGLAHDEMKSPREPGYKAEGKYKDDPVDLRLDIERRKKHKERDLKRGKSRESVDSRDSSHSRERSAEKTEKTHKGSKKQKKHRRARDRSRSSSSSSQSSHSYKAEEYTEETEEREESTTGFDKSRLGTKDFVGPSERGGGRARGTFQFRARGRGWGRGNYSGNNNNNSNNDFQKRNREEEWDPEYTPKSKKYYLHDDREGEGSDKWVSRGRGRGAFPRGRGRFMFRKSSTSPKWAHDKFSGEEGEIEDDESGTENREEKDNIQPTTE +ENSP00000318415 pep:known chromosome:GRCh37:1:89445139:89458455:-1 gene:ENSG00000213516 transcript:ENST00000321792 gene_biotype:protein_coding transcript_biotype:protein_coding MVEADRPGKLFIGGLNTETNEKALETVFGKYGRIVEVLLIKDRETNKSRGFAFVTFESPADAKDAARDMNGKSLDGKAIKVEQATKPSFERGRHGPPPPPRSRGPPRGFGAGRGGSGGTRGPPSRGGHMDDGGYSMNFNMSSSRGPLPVKRGPPPRSGGPSPKRSAPSGLVRSSSGMGGRAPLSRGRDSYGGPPRREPLPSRRDVYLSPRDDGYSTKDSYSSRDYPSSRDTRDYAPPPRDYTYRDYGHSSSRDDYPSRGYGDRDGYGRDRDYSDHPSGGSYRDSYESYGNSRSAPLTRGPPPSYGGSSRYDDYSSSRDGYGGSRDSYSSSRSDLYSSCDRVGRQERGLPPSVERGYPSSRDSYSSSSRGAPRGAGPGGSRSDRGGGRSRY +ENSP00000446099 pep:known chromosome:GRCh37:1:89445142:89458643:-1 gene:ENSG00000213516 transcript:ENST00000399794 gene_biotype:protein_coding transcript_biotype:protein_coding MVEADRPGKLFIGGLNTETNEKALETVFGKYGRIVEVLLIKDRETNKSRGFAFVTFESPADAKDAARDMNGKSLDGKAIKVEQATKPSFERGRHGPPPPPRSRGPPRGFGAGRGGSGGTRGPPSRGGHMDDGGYSMNFNMSSSRGPLPVKRGPPPRSGGPSPKRSAPSGLVRSSSGMGGRAPLSRGRDSYGGPPRREPLPSRRDVYLSPRDDGYSTKDSYSSRDYPSSRDTRDYAPPPRDYTYRDYGHSSSRDDYPSRGYGDRDGYGRDRDYSDHPSGGSYRDSYESYGNSRSAPLTRGPPPSYGGSSRYDDYSSSRDGYGGSRDSYSSSRSDLYSSCDRVGRQERGLPPSVERGYPSSRDSYSSSSRGAPRGAGPGGSRSDRGGGRSRY
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/peptides_10_10.fa Tue Dec 17 18:45:13 2013 -0500 @@ -0,0 +1,12 @@ +>ENSP00000363782 pep:sap chromosome:GRCh37:1:22778472:22853855:1 gene:ENSG00000184677 transcript:ENST00000374651 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:22846709 codon_change:Gtg/Atg sap:V885M +FSAPSMLERHMVTHVGGKPFS +>ENSP00000346634 pep:sap chromosome:GRCh37:1:36690017:36770958:1 gene:ENSG00000054118 transcript:ENST00000354618 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:36752433 codon_change:gCc/gTc sap:A201V +QAAGDNQGDEVKEQTFSGGTS +>ENSP00000318415 pep:sap chromosome:GRCh37:1:89445139:89458455:-1 gene:ENSG00000213516 transcript:ENST00000321792 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:89449390 codon_change:atA/atG sap:I40M +KYGRIVEVLLMKDRETNKSRG +>ENSP00000446099 pep:sap chromosome:GRCh37:1:89445142:89458643:-1 gene:ENSG00000213516 transcript:ENST00000399794 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:89449390 codon_change:atA/atG sap:I40M +KYGRIVEVLLMKDRETNKSRG +>ENSP00000318415 pep:sap chromosome:GRCh37:1:89445139:89458455:-1 gene:ENSG00000213516 transcript:ENST00000321792 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:89449434 codon_change:Aca/Gca sap:T26A +NTETNEKALEAVFGKYGRIVE +>ENSP00000446099 pep:sap chromosome:GRCh37:1:89445142:89458643:-1 gene:ENSG00000213516 transcript:ENST00000399794 gene_biotype:protein_coding transcript_biotype:protein_coding snp_location:chr1:89449434 codon_change:Aca/Gca sap:T26A +NTETNEKALEAVFGKYGRIVE
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/snpeff.vcf Tue Dec 17 18:45:13 2013 -0500 @@ -0,0 +1,35 @@ +##fileformat=VCFv4.1 +##samtoolsVersion=0.1.18 (r982:295) +##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth"> +##INFO=<ID=DP4,Number=4,Type=Integer,Description="# high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases"> +##INFO=<ID=MQ,Number=1,Type=Integer,Description="Root-mean-square mapping quality of covering reads"> +##INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability of all samples being the same"> +##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)"> +##INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)"> +##INFO=<ID=G3,Number=3,Type=Float,Description="ML estimate of genotype frequencies"> +##INFO=<ID=HWE,Number=1,Type=Float,Description="Chi^2 based HWE test P-value based on G3"> +##INFO=<ID=CLR,Number=1,Type=Integer,Description="Log ratio of genotype likelihoods with and without the constraint"> +##INFO=<ID=UGT,Number=1,Type=String,Description="The most probable unconstrained genotype configuration in the trio"> +##INFO=<ID=CGT,Number=1,Type=String,Description="The most probable constrained genotype configuration in the trio"> +##INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for strand bias, baseQ bias, mapQ bias and tail distance bias"> +##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL."> +##INFO=<ID=PC2,Number=2,Type=Integer,Description="Phred probability of the nonRef allele frequency in group1 samples being larger (,smaller) than in group2."> +##INFO=<ID=PCHI2,Number=1,Type=Float,Description="Posterior weighted chi^2 P-value for testing the association between group1 and group2 samples."> +##INFO=<ID=QCHI2,Number=1,Type=Integer,Description="Phred scaled PCHI2."> +##INFO=<ID=PR,Number=1,Type=Integer,Description="# permutations yielding a smaller PCHI2."> +##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality"> +##FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)"> +##FORMAT=<ID=DP,Number=1,Type=Integer,Description="# high-quality bases"> +##FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value"> +##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods"> +##SnpEffVersion="3.4 (build 2013-11-27), by Pablo Cingolani" +##SnpEffCmd="SnpEff -i vcf -o vcf -upDownStreamLen 5000 -spliceSiteSize 1 -snp -no-downstream -no-intergenic -no-intron -no-upstream -no-utr -stats /galaxy/PRODUCTION/database/files/000/376/dataset_376197.dat GRCh37.73 /galaxy/PRODUCTION/database/files/000/376/dataset_376194.dat " +##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 ] )' "> +##SnpEffCmd="SnpEff -i vcf -o vcf -upDownStreamLen 5000 -spliceSiteSize 1 -geneId -no-downstream -no-intergenic -no-intron -no-upstream -no-utr -stats /galaxy/PRODUCTION/database/files/000/376/dataset_376199.dat GRCh37.73 /galaxy/PRODUCTION/database/files/000/376/dataset_376196.dat " +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT /galaxy/PRODUCTION/database/tmp/tmp-SAMTOOLS-BmXdVn/bam_input_0.bam +chr1 22846709 . G A 9.31 . DP=2;VDB=0.0174;AF1=1;AC1=2;DP4=0,0,1,1;MQ=20;FQ=-33;EFF=NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|Gtg/Atg|V885M|1127|ZBTB40|protein_coding|CODING|ENST00000374651|12|1) PL 40,6,0 +chr1 36752433 . C T 6.02 . DP=2;VDB=0.0099;AF1=1;AC1=2;DP4=0,0,2,0;MQ=20;FQ=-33;EFF=NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|gCc/gTc|A201V|955|THRAP3|protein_coding|CODING|ENST00000354618|4|1) PL 36,6,0 +chr1 89449390 . T C 15.1 . DP=3;VDB=0.0214;AF1=1;AC1=2;DP4=0,0,1,2;MQ=20;FQ=-36;EFF=EXON(MODIFIER|||||RBMXL1|processed_transcript|CODING|ENST00000413769|3|1),NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|atA/atG|I40M|390|RBMXL1|protein_coding|CODING|ENST00000321792|2|1),NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|atA/atG|I40M|390|RBMXL1|protein_coding|CODING|ENST00000399794|3|1);EFF=EXON(MODIFIER|||||ENSG00000213516|processed_transcript|CODING|ENST00000413769|3|1) PL 47,9,0 +chr1 89449434 . T C 9.31 . DP=2;VDB=0.0120;AF1=1;AC1=2;DP4=0,0,1,1;MQ=20;FQ=-33;EFF=EXON(MODIFIER|||||RBMXL1|processed_transcript|CODING|ENST00000413769|3|1),NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|Aca/Gca|T26A|390|RBMXL1|protein_coding|CODING|ENST00000321792|2|1),NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|Aca/Gca|T26A|390|RBMXL1|protein_coding|CODING|ENST00000399794|3|1);EFF=EXON(MODIFIER|||||ENSG00000213516|processed_transcript|CODING|ENST00000413769|3|1) PL 40,6,0