Question Regarding MSMC Input File Generation Without Individual Masks Using GATK

[Luker121]

Hi,
thanks for the comprehensive tutorial you have provided. I have a question regarding the elif statement in the msmc_2_generateInput_multiInd.shfile. Particularly about the section that mentions: "Creating MSMC input file WITHOUT individual mask (gatk)". From this description, I assume that individual mask files are not required when the VCF files are generated using GATK. Could you please clarify if my understanding is correct? If the use of individual mask files is indeed not necessary in the context of GATK, could you please explain the reasons of this decision?

Thank you for your help in advance.

[jessicarick]

Thanks for the question, and I'm glad to hear that you're finding the tutorial useful. The short answer is that using GATK does not necessarily preclude the need for an individual mask. The elif statement that you mention is mostly a relic from testing out different variant callers and using the msmc-tools scripts from the MSMC2 authors (https://github.com/stschiff/msmc-tools/tree/master). The msmc-tools script bamCaller.py uses samtools, and will generate both a VCF and an individual mask file from your bamfiles, and this is the method that I have mostly used (and thus the relic elif statement that has hung around). If you want instead to create your VCF file another way (e.g., using GATK), then you can use the vcfAllSiteParser.py script to generate individual mask files. You can generate input without using an individual mask, but my impression is that it is not recommended.

Thanks for pointing this out to me, and I'll get the script updated so that it includes the step to generate individual mask files the alternative way as well!

Jessi

[plubbe]

Hi Jessi -

First, thank you very much for your code and the helpful tutorial files you have created.
I have been having quite a lot of issues getting the scripts to run for an analysis on multiple individuals based on GATK outputs - specifically, I have not been able to use msmc_2_generateInput_multiInd.sh.
To begin, I am not sure how to modify the early for loop in that script to work with VCFs generated by a collaborator:
for ind in $IND do INDMASK='ls ${OUTDIR}/mask/ind_mask.${ind}.${SCAFFOLD}.bed.gz' echo "--mask=$INDMASK " >> ${OUTDIR}/mask/${POP}.mask_file.$SCAFFOLD INDVCF='ls ${OUTDIR}/vcf/${ind}.${SCAFFOLD}.msmc.vcf.gz' echo $INDVCF >> ${OUTDIR}/vcf/${POP}.vcf_file.${SCAFFOLD} done
as the files ${OUTDIR}/mask/ind_mask.${ind}.${SCAFFOLD}.bed.gz do not exist; they are seemingly created in msmc_1 during variant calling, but I do not wish to (not would I be able to) repeat the variant calling process just to generate the files. This looks like an unnecessary step as the gatk method section seems to contain code for creating individual masks using vcfAllSiteParser.py, in any case. However, this results in almost exactly the same knock-on problem as the files called in that step don't exist, either, nor can I find any code that generates them elsewhere:
MASK_INDIV=${OUTDIR}/mask/ind_mask.${IND}.${SCAFFOLD}.${METHOD}.bed.gz
Regardless, even if I hash out the INDMASK portion of the above loop (i.e., run only the INDVCF/echo $INDVCF portion) and continue to run the next command:
generate_multihetsep.py --mask=$MASK_GENOME 'cat ${OUTDIR}/vcf/${POP}.vcf_file.${SCAFFOLD}' > $MSMC_INPUT
by pointing to my (scaffold and individual separated) GATK VCFs, I get an error I don't understand:
Traceback (most recent call last):
File "msmc_tools/generate_multihetsep.py", line 225, in
for chrom, snp_pos, alleles in joinedVcfIterator:
File "msmc_tools/generate_multihetsep.py", line 153, in next
self.current_lines[i] = next(self.vcfIterators[i])
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "msmc_tools/generate_multihetsep.py", line 81, in next
return (chrom, pos, tuple(alleles), (int(geno[0]), int(geno[2])), phased)
^^^^^^^^^^^^
ValueError: invalid literal for int() with base 10: '.'

I have really been struggling to make this pipeline work with my files and would appreciate any insight you might have so I can resolve this pernicious problem that has plagued me for weeks. Thanks again for your really excellent documentation, it's wonderfully written really.