search

openvax / isovar

Assembly of RNA reads to determine the effect of a cancer mutation on protein sequence
Apache License 2.0
25 stars 13 forks source link

Variant sequence ... shorter than min allowed #99

Open gilhornung opened 6 years ago

gilhornung commented 6 years ago

Hi,

I've tried to run isovar-protein-sequences.py and noticed that a few variants with RNA support were missing. Below is an example of two genes, THSD1 and ARHGEF10, in which a variant was skipped in isovar. The stdout has an error Variant sequence ... shorter than min allowed, and Unable to determine reading frame for VariantSequence.... I tried to reduce the --min-transcript-prefix-length, but I still get the error.

All the files to reproduce the error can be found in the link https://owncloud.incpm.weizmann.ac.il/owncloud/index.php/s/L83i1GsOXqO1oOt

python isovar-protein-sequences.py --vcf variants.vcf --genome GRCh38 --bam rna.bam --min-alt-rna-reads 1 --min-transcript-prefix-length 1

2018-06-03 12:22:43,374 - isovar.allele_reads:200 - INFO - Gathering reads for Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38')
2018-06-03 12:22:43,414 - isovar.allele_reads:208 - INFO - Gathering variant reads for variant Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38') (chromosome = 13, gene names = [u'THSD1'])
2018-06-03 12:22:43,438 - isovar.locus_reads:314 - INFO - Found 9 reads overlapping locus 13: 52377690-52377692
2018-06-03 12:22:43,438 - isovar.variant_sequences:422 - INFO - Initial pool of 2 variant sequences (min length=63, max length=63)
2018-06-03 12:22:43,439 - isovar.assembly:140 - INFO - Collapsed 2 -> 2 sequences
2018-06-03 12:22:43,439 - isovar.variant_sequences:440 - INFO - After overlap assembly: 2 variant sequences (min length=63, max length=63)
2018-06-03 12:22:43,440 - isovar.variant_sequences:207 - INFO - Coverage: [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1] (len=63)
2018-06-03 12:22:43,442 - isovar.variant_sequences:207 - INFO - Coverage: [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1] (len=63)
2018-06-03 12:22:43,443 - isovar.variant_sequences:330 - INFO - Kept 1/2 variant sequences after read coverage trimming to >=2x
2018-06-03 12:22:43,443 - isovar.variant_sequences:456 - INFO - After coverage & length filtering: 1 variant sequences (min length=0, max length=0)
2018-06-03 12:22:44,482 - pyensembl.sequence_data:118 - INFO - Loaded sequence dictionary from /users/pmgil/.cache/pyensembl/GRCh38/ensembl87/Homo_sapiens.GRCh38.cdna.all.fa.gz.pickle
2018-06-03 12:22:44,617 - pyensembl.sequence_data:118 - INFO - Loaded sequence dictionary from /users/pmgil/.cache/pyensembl/GRCh38/ensembl87/Homo_sapiens.GRCh38.ncrna.fa.gz.pickle
2018-06-03 12:22:45,031 - pyensembl.sequence_data:118 - INFO - Loaded sequence dictionary from /users/pmgil/.cache/pyensembl/GRCh38/ensembl87/Homo_sapiens.GRCh38.pep.all.fa.gz.pickle
2018-06-03 12:22:45,036 - isovar.effect_prediction:65 - INFO - Predicted total 2 effects for variant Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38')
2018-06-03 12:22:45,036 - isovar.effect_prediction:73 - INFO - Keeping 2/2 effects which affect protein coding sequence for Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38'): <EffectCollection with 2 elements>
  -- Substitution(variant=Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38'), transcript_name=THSD1-001, transcript_id=ENST00000258613, effect_description=p.R760P)
  -- Substitution(variant=Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38'), transcript_name=THSD1-002, transcript_id=ENST00000349258, effect_description=p.R707P)
2018-06-03 12:22:45,036 - isovar.effect_prediction:84 - INFO - Keeping 2 effects with predictable AA sequences for Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38'): [Substitution(variant=Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38'), transcript_name=THSD1-001, transcript_id=ENST00000258613, effect_description=p.R760P), Substitution(variant=Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38'), transcript_name=THSD1-002, transcript_id=ENST00000349258, effect_description=p.R707P)]
2018-06-03 12:22:45,037 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on THSD1-001 for variant Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38') = 2457:2458
2018-06-03 12:22:45,037 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on THSD1-002 for variant Variant(contig='13', start=52377691, ref='C', alt='G', reference_name='GRCh38') = 2664:2665
2018-06-03 12:22:45,037 - isovar.variant_sequence_in_reading_frame:334 - INFO - Variant sequence VariantSequence(prefix='', alt='', suffix='', sequence='', reads=frozenset([AlleleRead(prefix='GGGGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGA', allele='G', suffix='GAGCTCTGTGGGGCTCTGTTCTCTCAATTCCGGCCACTAATCCTG', name='A00376:25:HCKVGDMXX:1:2341:4535:8155', sequence='GGGGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGAGGAGCTCTGTGGGGCTCTGTTCTCTCAATTCCGGCCACTAATCCTG'), AlleleRead(prefix='GGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGA', allele='G', suffix='GAGCTCTGGGGGGCTCTGTTCTCTCAATTCCGGC', name='A00376:25:HCKVGDMXX:1:2426:22553:6073', sequence='GGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGAGGAGCTCTGGGGGGCTCTGTTCTCTCAATTCCGGC')])) shorter than min allowed 1 (iter=1)
2018-06-03 12:22:45,037 - isovar.translation:206 - INFO - Unable to determine reading frame for VariantSequence(prefix='', alt='', suffix='', sequence='', reads=frozenset([AlleleRead(prefix='GGGGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGA', allele='G', suffix='GAGCTCTGTGGGGCTCTGTTCTCTCAATTCCGGCCACTAATCCTG', name='A00376:25:HCKVGDMXX:1:2341:4535:8155', sequence='GGGGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGAGGAGCTCTGTGGGGCTCTGTTCTCTCAATTCCGGCCACTAATCCTG'), AlleleRead(prefix='GGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGA', allele='G', suffix='GAGCTCTGGGGGGCTCTGTTCTCTCAATTCCGGC', name='A00376:25:HCKVGDMXX:1:2426:22553:6073', sequence='GGGATATGGGAGAAGACTGCTTCCTTGAGACACTCCTGTGACTGGGGGACGGTCCCCGAGGAGCTCTGGGGGGCTCTGTTCTCTCAATTCCGGC')]))
2018-06-03 12:22:45,038 - isovar.allele_reads:200 - INFO - Gathering reads for Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38')
2018-06-03 12:22:45,038 - isovar.allele_reads:208 - INFO - Gathering variant reads for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') (chromosome = 8, gene names = [u'ARHGEF10', u'RP11-439C15.6'])
2018-06-03 12:22:45,070 - isovar.locus_reads:314 - INFO - Found 21 reads overlapping locus 8: 1903317-1903319
2018-06-03 12:22:45,070 - isovar.variant_sequences:422 - INFO - Initial pool of 1 variant sequences (min length=38, max length=38)
2018-06-03 12:22:45,070 - isovar.variant_sequences:440 - INFO - After overlap assembly: 1 variant sequences (min length=38, max length=38)
2018-06-03 12:22:45,071 - isovar.variant_sequences:207 - INFO - Coverage: [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1] (len=38)
2018-06-03 12:22:45,072 - isovar.variant_sequences:330 - INFO - Kept 1/1 variant sequences after read coverage trimming to >=2x
2018-06-03 12:22:45,072 - isovar.variant_sequences:456 - INFO - After coverage & length filtering: 1 variant sequences (min length=0, max length=0)
2018-06-03 12:22:45,166 - isovar.effect_prediction:48 - INFO - Skipping transcript ARHGEF10-006 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,166 - isovar.effect_prediction:48 - INFO - Skipping transcript ARHGEF10-004 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,167 - isovar.effect_prediction:48 - INFO - Skipping transcript ARHGEF10-005 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,167 - isovar.effect_prediction:48 - INFO - Skipping transcript RP11-439C15.6-001 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,175 - isovar.effect_prediction:65 - INFO - Predicted total 6 effects for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38')
2018-06-03 12:22:45,175 - isovar.effect_prediction:73 - INFO - Keeping 5/6 effects which affect protein coding sequence for Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'): <EffectCollection with 5 elements>
  -- Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-003, transcript_id=ENST00000349830, effect_description=p.R563T)
  -- Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-201, transcript_id=ENST00000398560, effect_description=p.R319T)
  -- Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-202, transcript_id=ENST00000398564, effect_description=p.R588T)
  -- Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-001, transcript_id=ENST00000518288, effect_description=p.R587T)
  -- Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-002, transcript_id=ENST00000520359, effect_description=p.R525T)
2018-06-03 12:22:45,176 - isovar.effect_prediction:84 - INFO - Keeping 5 effects with predictable AA sequences for Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'): [Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-003, transcript_id=ENST00000349830, effect_description=p.R563T), Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-201, transcript_id=ENST00000398560, effect_description=p.R319T), Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-202, transcript_id=ENST00000398564, effect_description=p.R588T), Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-001, transcript_id=ENST00000518288, effect_description=p.R587T), Substitution(variant=Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38'), transcript_name=ARHGEF10-002, transcript_id=ENST00000520359, effect_description=p.R525T)]
2018-06-03 12:22:45,176 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-003 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') = 1872:1873
2018-06-03 12:22:45,176 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-201 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') = 1442:1443
2018-06-03 12:22:45,177 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-202 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') = 1762:1763
2018-06-03 12:22:45,177 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-001 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') = 1922:1923
2018-06-03 12:22:45,178 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-002 for variant Variant(contig='8', start=1903318, ref='G', alt='C', reference_name='GRCh38') = 1746:1747
2018-06-03 12:22:45,178 - isovar.variant_sequence_in_reading_frame:334 - INFO - Variant sequence VariantSequence(prefix='', alt='', suffix='', sequence='', reads=frozenset([AlleleRead(prefix='CCGACA', allele='C', suffix='TCTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG', name='A00376:25:HCKVGDMXX:1:1404:23421:10551', sequence='CCGACACTCTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG')])) shorter than min allowed 1 (iter=1)
2018-06-03 12:22:45,178 - isovar.translation:206 - INFO - Unable to determine reading frame for VariantSequence(prefix='', alt='', suffix='', sequence='', reads=frozenset([AlleleRead(prefix='CCGACA', allele='C', suffix='TCTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG', name='A00376:25:HCKVGDMXX:1:1404:23421:10551', sequence='CCGACACTCTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG')]))
2018-06-03 12:22:45,178 - isovar.allele_reads:200 - INFO - Gathering reads for Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38')
2018-06-03 12:22:45,178 - isovar.allele_reads:208 - INFO - Gathering variant reads for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') (chromosome = 8, gene names = [u'ARHGEF10', u'RP11-439C15.6'])
2018-06-03 12:22:45,180 - isovar.locus_reads:314 - INFO - Found 21 reads overlapping locus 8: 1903318-1903320
2018-06-03 12:22:45,181 - isovar.variant_sequences:422 - INFO - Initial pool of 1 variant sequences (min length=39, max length=39)
2018-06-03 12:22:45,181 - isovar.variant_sequences:440 - INFO - After overlap assembly: 1 variant sequences (min length=39, max length=39)
2018-06-03 12:22:45,181 - isovar.variant_sequences:207 - INFO - Coverage: [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1] (len=39)
2018-06-03 12:22:45,181 - isovar.variant_sequences:330 - INFO - Kept 1/1 variant sequences after read coverage trimming to >=2x
2018-06-03 12:22:45,182 - isovar.variant_sequences:456 - INFO - After coverage & length filtering: 1 variant sequences (min length=0, max length=0)
2018-06-03 12:22:45,185 - isovar.effect_prediction:48 - INFO - Skipping transcript ARHGEF10-006 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,185 - isovar.effect_prediction:48 - INFO - Skipping transcript ARHGEF10-004 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,185 - isovar.effect_prediction:48 - INFO - Skipping transcript ARHGEF10-005 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,185 - isovar.effect_prediction:48 - INFO - Skipping transcript RP11-439C15.6-001 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') because it's incomplete
2018-06-03 12:22:45,185 - isovar.effect_prediction:65 - INFO - Predicted total 6 effects for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38')
2018-06-03 12:22:45,185 - isovar.effect_prediction:73 - INFO - Keeping 5/6 effects which affect protein coding sequence for Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'): <EffectCollection with 5 elements>
  -- Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-003, transcript_id=ENST00000349830, effect_description=p.R563S)
  -- Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-201, transcript_id=ENST00000398560, effect_description=p.R319S)
  -- Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-202, transcript_id=ENST00000398564, effect_description=p.R588S)
  -- Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-001, transcript_id=ENST00000518288, effect_description=p.R587S)
  -- Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-002, transcript_id=ENST00000520359, effect_description=p.R525S)
2018-06-03 12:22:45,186 - isovar.effect_prediction:84 - INFO - Keeping 5 effects with predictable AA sequences for Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'): [Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-003, transcript_id=ENST00000349830, effect_description=p.R563S), Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-201, transcript_id=ENST00000398560, effect_description=p.R319S), Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-202, transcript_id=ENST00000398564, effect_description=p.R588S), Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-001, transcript_id=ENST00000518288, effect_description=p.R587S), Substitution(variant=Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38'), transcript_name=ARHGEF10-002, transcript_id=ENST00000520359, effect_description=p.R525S)]
2018-06-03 12:22:45,186 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-003 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') = 1873:1874
2018-06-03 12:22:45,187 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-201 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') = 1443:1444
2018-06-03 12:22:45,187 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-202 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') = 1763:1764
2018-06-03 12:22:45,187 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-001 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') = 1923:1924
2018-06-03 12:22:45,188 - isovar.reference_sequence_key:113 - INFO - Interbase offset range on ARHGEF10-002 for variant Variant(contig='8', start=1903319, ref='G', alt='T', reference_name='GRCh38') = 1747:1748
2018-06-03 12:22:45,188 - isovar.variant_sequence_in_reading_frame:334 - INFO - Variant sequence VariantSequence(prefix='', alt='', suffix='', sequence='', reads=frozenset([AlleleRead(prefix='CCGACAC', allele='T', suffix='CTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG', name='A00376:25:HCKVGDMXX:1:1404:23421:10551', sequence='CCGACACTCTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG')])) shorter than min allowed 1 (iter=1)
2018-06-03 12:22:45,188 - isovar.translation:206 - INFO - Unable to determine reading frame for VariantSequence(prefix='', alt='', suffix='', sequence='', reads=frozenset([AlleleRead(prefix='CCGACAC', allele='T', suffix='CTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG', name='A00376:25:HCKVGDMXX:1:1404:23421:10551', sequence='CCGACACTCTGCCTCTTCAGATGGCCCTGACAGAGCTCGAAACACTAGCAGAGAAGTTAAATGAAAGAAAGAGAGATGCTGATCAACGCTGTGAAGTGAAG')]))