Missing Data

[mullerbsf]

Dear Vikas Bansal,

I ran CRISP with the following command:

CRISP --bams ${LIST_OF_BAMS} --ref ${FASTAREF} --regions ${REGIONS} --poolsize 50 --perms 50000 --mmq 10 --minc 3 \ --EM 1 --VCF ${FILE_CODE}_snps_and_indels_pt${SLURM_ARRAY_TASK_ID}.vcf > variantcalls_pt${SLURM_ARRAY_TASK_ID}.log

However, of the total 172696 variants detected, only 29973 does not have 100% of missing data. When I filtered for 40% of missing data, it returned only 6567 variants.

Why CRISP is calling a lot of missing variants, but reads are covering that particular position (example below for SNP: Chla_rei_1_19528)?

Thank you!

Example:

CRISP raw output:

1 19528 Chla_rei_1_19528 G C 6186 LowDepth NP=38;DP=485,301,12;VT=SNV;CT=-0.4;VP=36;VF=EMpass;AC=500;AF=0.26742;EMstats=618.69:-245.27;HWEstats=0.0;MQS=0,0,0,1097;FLANKSEQ=aaccttcact:G:ccatcgtcgc AC:GQ:DP:ADf:ADr:ADb .:0:50:16,5:9,4:0,1 .:0:33:13,3:7,2:0,0 .:0:28:11,4:6,3:0,0 .:0:43:11,6:13,2:0,0 .:0:20:7,3:4,0:0,0 .:0:43:16,8:6,1:0,0 .:0:30:16,1:3,2:1,0 .:0:15:3,2:3,2:0,0 .:0:23:7,1:8,3:0,1

CRISP converted output (using convert_pooled_vcf.py):

1 19528 Chla_rei_1_19528 G C 6186 LowDepth NP=38;DP=485,301,12;VT=SNV;CT=-0.4;VP=36;VF=EMp ass;AC=500;AF=0.26742;EMstats=618.69:-245.27;HWEstats=0.0;MQS=0,0,0,1097;FLANKSEQ=aaccttcact:G:ccatcgtcgc GT:GQ:DP:ADf:A Dr:ADb ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:50:16,5:9,4:0,1 ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:33:13,3:7,2:0,0 ././ ./././././././././././././././././././././././././././././././././././././././././././././././.:0:28:11,4:6,3:0,0 ./././././. /./././././././././././././././././././././././././././././././././././././././././././.:0:43:11,6:13,2:0,0 ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:20:7,3:4,0:0,0 ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:43:16,8:6,1:0,0 ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:30:16,1:3,2:1,0 ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:15:3,2:3,2:0,0 ./././././././././././././././././././././././././././././././././././././././././././././././././.:0:23:7,1:8,3:0,1

[vibansal]

CRISP will not call genotypes for the pools when the coverage for a pool is less than 2 times the pool-size since the allele counts cannot be estimated reliably.

[mullerbsf]

Thanks for the answer.

In my case, I used the --poolsize 50. However, the average depth coverage was only 16. Do I need to use the --poolsize 8 or 10 to reduce the number of missing data?

I appreciate your help!

[vibansal]

Yes, that will reduce the frequency of missing calls.