Error: unhandled exception: comphet.nim(46, 12)

[Drdreammaerd]

Hi @brentp,

I ran the comphet function after calling the denovo and comp-side, and it gave me fatal error that I don't know how to fix. I have ran the same command line for another vcf, which worked pretty well.

slivar compound-hets -v $vcf --sample-field comphet_side --sample-field denovo -p $ped > $out

Error fatal.nim(49) sysFatal Error: unhandled exception: comphet.nim(46, 12) kid.dad != nil and kid.mom != nil and kid.dad.i != -1 and kid.mom.i != -1 [AssertionDefect]

Thank you for your time and any suggestion is aprreciated,

David

[brentp]

Hi, this looks like a problem with your pedigree file. Can you share that here or send off-line?

[Drdreammaerd]

Hi Brent,

The attached file is the ped file and group file. For calling denovo and compside, I used group file.

Command for denovo and compside

slivar expr --vcf $vcf --alias $group --gnotate $gnomad --js $jsycw
--out-vcf $out \
    --info "variant.call_rate > 0.9" \
    --group-expr "denovo:denovo(kid, mom, dad) && INFO.ExAC_ALL < 0.0004" \
    --group-expr "comphet_side:comphet_side(kid, mom, dad) && (
INFO.ExAC_ALL < 0.001 || INFO.x3b1000g2015aug_all < 0.001 ||
INFO.esp6500siv2_all < 0.001) "

Command for comphet

slivar compound-hets -v $vcf \
    --sample-field comphet_side --sample-field denovo -p $ped > $out

Thanks,

David ---------------------------------------------PED--------------------------------------------------- AUS005 AUS005P AUS005F AUS005M 2 2 AUS007 AUS007P AUS007F AUS007M 2 2 AUS008 AUS008P AUS008F AUS008M 1 2 AUS009 AUS009P AUS009F AUS009M 1 2 AUS012 AUS012P AUS012F AUS012M 2 2 AUS017 AUS017P AUS017F AUS017M 1 2 AUS019 AUS019P AUS019F AUS019M 1 2 ------------------------------------------group----------------------------------------------------- #kid dad mom AUS005P AUS005F AUS005M AUS007P AUS007F AUS007M AUS008P AUS008F AUS008M AUS009P AUS009F AUS009M AUS012P AUS012F AUS012M AUS017P AUS017F AUS017M AUS019P AUS019F AUS019M

[brentp]

ok. that's too much for me to guess. I've added more info to the error message. can you try again with this attached binary (just gunzip and chmod+x) so we can narrow it down.

I'm curious why you chose to use group-expr rather than family-expr?

You can ask other questions here slivar.gz

[Drdreammaerd]

I am just wondering where should I put the file since I am using slivar in docker.

The reason I am using group-expr, just because I thought the group file is more easy to generate for a large size cohort. Can I also use group file for family-expr? Also, if there is singleton, dual, quad. How should I indicate them in group and ped that is able to recognize by family-expr?

Thanks,

David

[brentp]

if your pedigree file is properly speicified, then slivar will automatically detect familys and the --family-expr stuff shown in the wiki will work for families of all sizes, including duos, multi-generational pedigrees, and so on.

I am just wondering where should I put the file since I am using slivar in docker.

which file? just mount the directory containing whichever file you need.

[Drdreammaerd]

I am sorry if my question is not clear. Kind of new to the area. Will do my best to clarify my question. What I mean is the binary file called "slivar.gz" you attached. Do you mind to give me more detailed instruction how to use it? I am using the docker through institution cluster

Thanks,

David

[brentp]

oh. I see. so unzip the file and chmod +x it. Then if you put in in your working directory you can do, for example:

docker run -v$(pwd):/load ...  /load/slivar expr ... 

[Drdreammaerd]

Cool, Thanks. I have chomd +x it

-rwxr-xr-x@ 1 dreammaerd  staff    17M Dec 20 09:22 slivar_brent

How do I execute it?

command not found: slivar_brent

[brentp]

./slivar_brent or use the linked path as shown in the example above.

[Drdreammaerd]

Thanks brent. My colleague helps to find the problem and it works well now. The problem is in original pedigree file I didn't provide the information of parents. Somehow it works when I provide those.

Original PED

AUS005  AUS005P AUS005F AUS005M 2       2
AUS007  AUS007P AUS007F AUS007M 2       2
AUS008  AUS008P AUS008F AUS008M 1       2
AUS009  AUS009P AUS009F AUS009M 1       2
AUS012  AUS012P AUS012F AUS012M 2       2

New PED

AUS005  AUS005P AUS005F AUS005M 2       2
AUS005  AUS005F 0       0       1       1
AUS005  AUS005M 0       0       2       1
AUS007  AUS007P AUS007F AUS007M 2       2
AUS007  AUS007F 0       0       1       1
AUS007  AUS007M 0       0       2       1
AUS008  AUS008P AUS008F AUS008M 1       2
AUS008  AUS008F 0       0       1       1
AUS008  AUS008M 0       0       2       1
AUS009  AUS009P AUS009F AUS009M 1       2
AUS009  AUS009F 0       0       1       1
AUS009  AUS009M 0       0       2       1
AUS012  AUS012P AUS012F AUS012M 2       2
AUS012  AUS012F 0       0       1       1
AUS012  AUS012M 0       0       2       1

P.s. Slivar is really a wonderful tool!!!

Thanks for your kind response and help,

David

[Drdreammaerd]

I have another question for Comp-het function, can I ask here or should I request a new thread?

[brentp]

here is fine, we can move to a new thread if needed. glad to hear that slivar is useful for you! I think it makes sense how you solved. The mix of groups/aliases and ped sort of hid the error that usually would be prevented by using --family-expr.

[Drdreammaerd]

Cool thanks.

My question is what field of vep annotation it will be counted for slivar-compound-hets function? I know in the description, it said group by gene. But somehow when I use vep annotation only choose field gene, the slivar cannot run that vep annotated file by just exiting out without error message.

VEP command

/opt/vep/src/ensembl-vep/vep -i $vcf -o $out --cache --dir $cache --force_overwrite --vcf --fields "Gene"

[brentp]

can you show what VEP adds to the header of $out for that? and can you show a full CSQ field from the INFO for one variant?

[Drdreammaerd]

Sure. The header for VEP is CSQ: Chr and position

chr1_902926_T_C T

Only Gene field

CSQ=ENSG00000272438|

All field

CSQ=C|upstream_gene_variant|MODIFIER|
|ENSG00000272438|Transcript|ENST00000607769|lncRNA|||||||||||1908|1|||

[brentp]

I mean can you show the actual VCF header line?

[Drdreammaerd]

The meta-information line for VEP

##VEP="v105" time="2021-12-18 23:21:53" cache="./.vep/homo_sapiens/105_GRCh38" db="[email protected]" ensembl-variation=105.ac8178e ensembl=105.f357e33 ensembl-funcgen=105.660df8f ensembl-io=105.2a0a40c 1000genomes="phase3" COSMIC="92" ClinVar="202106" ESP="V2-SSA137" HGMD-PUBLIC="20204" assembly="GRCh38.p13" dbSNP="154" gencode="GENCODE 39" genebuild="2014-07" gnomAD="r2.1.1" polyphen="2.2.2" regbuild="1.0" sift="sift5.2.2"
##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Gene|IMAPCT">

And it adds to "INFO" column in vcf.

[brentp]

This can't work:

 ##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format:Gene|IMAPCT">

"IMAPCT" will not be understood and it needs standard VEP annotations in the output. These should be something like:

##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence type as predicted by VEP. Format: Consequence|Codons|Amino_acids|Gene|SYMBOL|Feature|EXON|PolyPhen|SIFT|Protein_position|BIOTYPE">

[Drdreammaerd]

Thanks. Can you share your VEP command please?