Skip to content

Commit c0042c3

Browse files
theferrit32theferrit32
committed
Include sequence value in output vcf if the length of that sequence is <= MAX_LITERAL_STATE_LENGTH, defaulted to 15
1 parent bce957a commit c0042c3

File tree

3 files changed

+33
-25
lines changed

3 files changed

+33
-25
lines changed

‎src/ga4gh/vrs/extras/annotator/vcf.py‎

Lines changed: 24 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
from ga4gh.vrs import VRS_VERSION, __version__
1717
from ga4gh.vrs.dataproxy import _DataProxy
1818
from ga4gh.vrs.extras.translator import AlleleTranslator
19-
from ga4gh.vrs.models import Allele, Range
19+
from ga4gh.vrs.models import Allele, Range, sequenceString
2020

2121
_logger = logging.getLogger(__name__)
2222

@@ -52,6 +52,19 @@ class FieldName(str, Enum):
5252
}
5353
)
5454

55+
# Short State.sequence will be included in output VCF if <= this value,
56+
# otherwise output will be emitted as the "." character.
57+
MAX_LITERAL_STATE_LENGTH = 15
58+
59+
60+
def _sequence_value_for_info(sequence: str | sequenceString) -> str | None:
61+
"""Return a literal sequence suitable for INFO fields if it is short enough."""
62+
if not sequence:
63+
return None
64+
seq = getattr(sequence, "root", sequence)
65+
seq_str = str(seq)
66+
return seq_str if len(seq_str) <= MAX_LITERAL_STATE_LENGTH else None
67+
5568

5669
def dump_alleles_to_pkl(alleles: list[Allele], output_pkl_path: Path) -> None:
5770
"""Create pkl file of dictionary mapping VRS IDs to ingested alleles.
@@ -394,24 +407,18 @@ def _get_vrs_object(
394407
# Common fields for all state types
395408
start = vrs_obj.location.start
396409
end = vrs_obj.location.end
410+
state = vrs_obj.state
397411

398412
# State-specific fields
399-
state_type = vrs_obj.state.type
400-
401-
if state_type == "LiteralSequenceExpression":
402-
# For LSE, populate sequence in STATES
403-
alt = (
404-
str(vrs_obj.state.sequence.root)
405-
if vrs_obj.state.sequence
406-
else None
407-
)
408-
# TODO TODO Leave in the `sequence` when the length of it is less than 15 characters
409-
elif state_type in (
413+
state_type = state.type
414+
alt = _sequence_value_for_info(getattr(state, "sequence", None))
415+
416+
if state_type in (
410417
"ReferenceLengthExpression",
411418
"LengthExpression",
412419
):
413420
# For RLE and LE, populate length
414-
length = vrs_obj.state.length
421+
length = state.length
415422
if length is None:
416423
err_msg = f"{state_type} requires a non-empty length: {vcf_coords}"
417424
raise VcfAnnotatorError(err_msg)
@@ -420,10 +427,11 @@ def _get_vrs_object(
420427
raise VcfAnnotatorError(err_msg)
421428
# For RLE only, also populate repeatSubunitLength
422429
if state_type == "ReferenceLengthExpression":
423-
repeat_subunit_length = vrs_obj.state.repeatSubunitLength
430+
repeat_subunit_length = state.repeatSubunitLength
424431
else:
425-
err_msg = f"Unsupported state type '{state_type}' for VCF annotation: {vcf_coords}"
426-
raise VcfAnnotatorError(err_msg)
432+
if state_type != "LiteralSequenceExpression":
433+
err_msg = f"Unsupported state type '{state_type}' for VCF annotation: {vcf_coords}"
434+
raise VcfAnnotatorError(err_msg)
427435

428436
vrs_field_data[FieldName.STARTS_FIELD].append(start)
429437
vrs_field_data[FieldName.ENDS_FIELD].append(end)

‎tests/extras/data/test_vcf_expected_altsonly_output.vcf‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -237,8 +237,8 @@
237237
##INFO=<ID=VRS_RepeatSubunitLengths,Number=A,Type=Integer,Description="The repeat subunit length values from ReferenceLengthExpression states for the GA4GH VRS Alleles corresponding to the GT indexes of the ALT alleles">
238238
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG001
239239
chr19 82664 . C T 50 PASS platforms=2;platformnames=10X,PacBio;datasets=2;datasetnames=10XChromiumLR,CCS15kb_20kb;callsets=2;callsetnames=10XLRGATK,CCS15kb_20kbDV;datasetsmissingcall=HiSeqPE300x,CCS15kb_20kb,CGnormal,IonExome,SolidSE75bp;callable=CS_CCS15kb_20kbDV_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable;filt=CS_CCS15kb_20kbDV_filt;arbitrated=TRUE;difficultregion=hg38.segdups_sorted_merged,lowmappabilityall;VRS_Allele_IDs=ga4gh:VA.b9AIUegs1SOYO7nKY8l01yllkyUns-u4;VRS_Starts=82663;VRS_Ends=82664;VRS_States=T;VRS_Lengths=.;VRS_RepeatSubunitLengths=. GT:PS:DP:ADALL:AD:GQ 0/1:.:154:0,0:0,0:120
240-
chr19 284350 . CA C 50 PASS platforms=4;platformnames=Illumina,10X,PacBio,CG;datasets=4;datasetnames=HiSeqPE300x,10XChromiumLR,CCS15kb_20kb,CGnormal;callsets=5;callsetnames=HiSeqPE300xGATK,10XLRGATK,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes;datasetsmissingcall=CCS15kb_20kb,IonExome,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;difficultregion=GRCh38_AllHomopolymers_gt6bp_imperfectgt10bp_slop5,GRCh38_SimpleRepeat_imperfecthomopolgt10_slop5;VRS_Allele_IDs=ga4gh:VA.a04jFsNg0bS0RMIWjKWSbwJS4_vp7S6x;VRS_Starts=284350;VRS_Ends=284366;VRS_States;VRS_Lengths=15;VRS_RepeatSubunitLengths=1 GT:PS:DP:ADALL:AD:GQ 0/1:.:422:117,101:81,75:356
241-
chr19 289464 . T TCACGCCTGTAATCC 50 PASS platforms=4;platformnames=Illumina,PacBio,CG,10X;datasets=4;datasetnames=HiSeqPE300x,CCS15kb_20kb,CGnormal,10XChromiumLR;callsets=6;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes,CCS15kb_20kbDV,10XLRGATK;datasetsmissingcall=IonExome,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;VRS_Allele_IDs=ga4gh:VA.ySvDptXfHB_9WEfu78v32DzBXJfwGgO7;VRS_Starts=289464;VRS_Ends=289466;VRS_States=CACGCCTGTAATCCCA;VRS_Lengths=.;VRS_RepeatSubunitLengths=. GT:PS:DP:ADALL:AD:GQ 0/1:.:518:94,98:116,137:785
240+
chr19 284350 . CA C 50 PASS platforms=4;platformnames=Illumina,10X,PacBio,CG;datasets=4;datasetnames=HiSeqPE300x,10XChromiumLR,CCS15kb_20kb,CGnormal;callsets=5;callsetnames=HiSeqPE300xGATK,10XLRGATK,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes;datasetsmissingcall=CCS15kb_20kb,IonExome,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;difficultregion=GRCh38_AllHomopolymers_gt6bp_imperfectgt10bp_slop5,GRCh38_SimpleRepeat_imperfecthomopolgt10_slop5;VRS_Allele_IDs=ga4gh:VA.a04jFsNg0bS0RMIWjKWSbwJS4_vp7S6x;VRS_Starts=284350;VRS_Ends=284366;VRS_States=AAAAAAAAAAAAAAA;VRS_Lengths=15;VRS_RepeatSubunitLengths=1 GT:PS:DP:ADALL:AD:GQ 0/1:.:422:117,101:81,75:356
241+
chr19 289464 . T TCACGCCTGTAATCC 50 PASS platforms=4;platformnames=Illumina,PacBio,CG,10X;datasets=4;datasetnames=HiSeqPE300x,CCS15kb_20kb,CGnormal,10XChromiumLR;callsets=6;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes,CCS15kb_20kbDV,10XLRGATK;datasetsmissingcall=IonExome,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;VRS_Allele_IDs=ga4gh:VA.ySvDptXfHB_9WEfu78v32DzBXJfwGgO7;VRS_Starts=289464;VRS_Ends=289466;VRS_States;VRS_Lengths=.;VRS_RepeatSubunitLengths=. GT:PS:DP:ADALL:AD:GQ 0/1:.:518:94,98:116,137:785
242242
chr19 28946400 . T C 50 PASS platforms=5;platformnames=Illumina,PacBio,CG,10X,Solid;datasets=5;datasetnames=HiSeqPE300x,CCS15kb_20kb,CGnormal,10XChromiumLR,SolidSE75bp;callsets=7;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbDV,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes,10XLRGATK,SolidSE75GATKHC;datasetsmissingcall=IonExome;callable=CS_HiSeqPE300xGATK_callable,CS_CCS15kb_20kbDV_callable,CS_10XLRGATK_callable,CS_CCS15kb_20kbGATK4_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_CCS15kb_20kbDV_filt,CS_CCS15kb_20kbGATK4_filt;VRS_Allele_IDs=ga4gh:VA.uV5O4M9zpiwk6sftOd-EDvtw_pkSAvdf;VRS_Starts=28946399;VRS_Ends=28946400;VRS_States=C;VRS_Lengths=.;VRS_RepeatSubunitLengths=. GT:PS:DP:ADALL:AD:GQ 1/1:.:874:0,275:115,378:502
243243
chr19 490414 . ACT A 50 PASS platforms=5;platformnames=Illumina,PacBio,CG,10X,Solid;datasets=5;datasetnames=HiSeqPE300x,CCS15kb_20kb,CGnormal,10XChromiumLR,SolidSE75bp;callsets=7;callsetnames=HiSeqPE300xGATK,CCS15kb_20kbDV,CCS15kb_20kbGATK4,CGnormal,HiSeqPE300xfreebayes,10XLRGATK,SolidSE75GATKHC;datasetsmissingcall=IonExome;callable=CS_HiSeqPE300xGATK_callable,CS_CCS15kb_20kbDV_callable,CS_CCS15kb_20kbGATK4_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_10XLRGATK_filt;VRS_Allele_IDs=ga4gh:VA.lok7a3lot_cvUyw626otpJi4yxk0X07v;VRS_Starts=490414;VRS_Ends=490416;VRS_States;VRS_Lengths=0;VRS_RepeatSubunitLengths=2 GT:PS:DP:ADALL:AD:GQ 0/1:.:821:163,158:239,220:1004
244244
chr19 54220024 . G *,A 50 PASS platforms=1;platformnames=PacBio;datasets=1;datasetnames=CCS15kb_20kb;callsets=1;callsetnames=CCS15kb_20kbGATK4;datasetsmissingcall=HiSeqPE300x,CCS15kb_20kb,10XChromiumLR,CGnormal,IonExome,SolidSE75bp;callable=CS_CCS15kb_20kbGATK4_callable;filt=CS_CCS15kb_20kbDV_filt,CS_10XLRGATK_filt,CS_HiSeqPE300xfreebayes_filt;difficultregion=HG001.hg38.300x.bam.bilkentuniv.010920.dups,hg38.segdups_sorted_merged;VRS_Allele_IDs=,ga4gh:VA.I7J3i1B36BACEUINcTwEh7uMv3I-PXT1;VRS_Starts=.,54220023;VRS_Ends=.,54220024;VRS_States=,A;VRS_Lengths=.,.;VRS_RepeatSubunitLengths=.,. GT:PS:DP:ADALL:AD:GQ 1/2:.:45:0,20,25:0,20,25:99

0 commit comments

Comments
 (0)