Issue with Generating VCF for PanGenie Using the Pipeline

[HuangYihang1222]

Dear Jana,

I encountered an issue while using the pipeline from this repository to generate the VCF files required for PanGenie.

Upon inspection, I found that the VCF files in the calls folder were generated correctly. The following files are present and seem normal:

all-haplotypes-callable.log
all-haplotypes.vcf
all-haplotypes-callable.vcf
sample1-hap0.vcf.gz
sample1-hap0.vcf.gz.tbi
sample1-hap1.vcf.gz
sample1-hap1.vcf.gz.tbi

However, in the multisample-vcfs folder, the callset-filtered.vcf and graph-filtered.vcf files appear to be problematic. They contain only header lines (lines starting with #) and no variant entries.

After checking the log file generated with:

nohup snakemake -j 50 --use-conda > snakemake.log 2>&1 &

and comparing it with the Snakefile, I found that there were no error messages in the log. However, it seems that the sort_bed and callable_regions steps were not executed.

snakemake.log

I would like to understand what might be causing this issue and how I can resolve it. Additionally, my config.json file is as follows:

{
        "reference": 
                {
                        "filename": "/home/huangyh/pangenome/1_cactus_construction/CN.fasta"
                },
        "assemblies":
                {
                        "sample1" : ["/home/huangyh/pangenome/1_cactus_construction/NA1.fasta", "/home/huangyh/pangenome/1_cactus_construction/NA2.fasta"]
                },
        "trios":
                {

                },
        "scripts": "scripts",
        "outdir" : "/home/huangyh/pangenome/4_pangenie/"
}

Thank you for your time and assistance!

Best regards,
Yihang Huang

[eblerjana]

Hi Yihang,

how were the assemblies generated that you are using as input? Are these haplotype-resolved assemblies, that is, do each of them contain only contigs related to one haplotype? Are there any variants present in the all-haplotypes.vcf and all-haplotypes-callable.vcf?

Best,
Jana

[HuangYihang1222]

Hi Jana,

My assemblies are derived from published genomes (haplotypes) downloaded from other papers, as well as my own assembled genomes (reference). For the input, I removed small contigs and kept only the chromosomes (18 in total).

Variants are present in both the all-haplotypes.vcf and all-haplotypes-callable.vcf, as shown in the figure.

all-haplotypes.vcf

all-haplotypes-callable.vcf

Best,
Yihang

[eblerjana]

Ok I guess the problem is that the current version of the pipeline assumes a human reference genome with chromosomes named as "chr1", "chr2" ... or "1", "2", ... So right now, every entry with a different chromosome name is removed. I'll see if I can fix this and make it more general

[HuangYihang1222]

I have resolved my issue by reviewing and modifying the Snakefile script, but I am not sure whether my logic is correct. My concerns mainly fall into two areas:

1. The issue I mentioned earlier.
2. After generating graph-filtered.vcf, I noticed that it was flagged as containing overlapping variants.

For the first issue, while carefully inspecting the intermediate files, I found that callset.vcf was correctly generated. However, in this section of the Snakefile:

# remove all variants where there is a mendelian conflict in at least one of the trios
# if no trios are given in config, the vcf does not change.
rule check_mendelian_consistency:

It seems that the intended logic described in the comment was not followed. If no trios.ped file was provided, instead of keeping callset.vcf unchanged as callset-filtered.vcf and proceeding to generate graph-filtered.vcf, the script instead produced an empty callset-filtered.vcf file.

To address this, I modified the rule as follows:

rule check_mendelian_consistency:
    input:
        vcf="{outdir}multisample-vcfs/callset.vcf",
        ped="{outdir}trios.ped",
        samples="{outdir}trio-samples.txt"
    output:
        vcf="{outdir}multisample-vcfs/callset-filtered.vcf",
        tsv="{outdir}multisample-vcfs/mendelian-consistency.tsv"
    log:
        "{outdir}multisample-vcfs/mendelian-consistency.log"
    conda:
        "env/calling.yml"
    shell:
        """
        if [ ! -s {input.ped} ]; then
            echo "Warning: trios.ped is empty. Skipping Mendelian consistency check." > {log}
            cp {input.vcf} {output.vcf}
            touch {output.tsv}
        else
            python3 {scripts}/mendelian-consistency.py filter -vcf {input.vcf} -samples {input.samples} -ped {input.ped} -o {output.tsv} 2> {log} 1> {output.vcf}
        fi
        """

After this modification, the pipeline ran successfully. I am not entirely sure if this approach is correct, but I believe it aligns with the intended logic described in the comments.

For the second issue, I followed the recommendations in "What should I do if my input VCF contains overlapping variants?" and used pipelines/run-from-callset to generate a VCF suitable for PanGenie input. However, I encountered a small issue in the Snakefile:

rule pangenie_genotyping:
    input:
        index='{outdir}/pangenie/indexing/',
        reads = lambda wildcards: config['reads'][wildcards.sample]
    output:
        '{outdir}/pangenie/{sample}_graph_genotyping.vcf'
    threads:
        24
    resources:
        mem_total_mb=200000,
        runtime_hrs=4,
        runtime_min=59
    log:
        '{outdir}/pangenie/pangenie-{sample}.log'
    benchmark:
        '{outdir}/benchmarks/pangenie-{sample}.txt'
    params:
        prefix = "{outdir}/pangenie/{sample}_graph"
    shell:
        "{pangenie_genotype} -i {input.reads} -f {input.index}/index -o {params.prefix} -j {threads} -t {threads} &> {log}"

Originally, {input.index}/index was missing a /, which caused subsequent scripts to fail in recognizing the path to index_path_segments.fasta, leading to an error. Adding the missing / resolved the issue, though I am not sure if this was due to my specific hardware or software environment.

Finally, I compared the genotyping results obtained from running the pipeline and the results obtained using pangenome.vcf and callset-biallelic.vcf.gz with PanGenie-index and PanGenie(both applying convert-to-biallelic.py in the end). As expected, I found that the results were completely consistent.

I would like to know whether my modifications to the script are correct. I also hope that my feedback can help improve the software further.

[eblerjana]

I pushed a fixed version of the pipeline. As I explained in my previous comment, the problem was that the mendelian consistency script removed all entries not located on chromosomes named as "chr1", "chr2" etc. independent of the ped file being empty or not. This is fixed now, as well as the issue with overlapping variants.

Note that now the config file has one additional field to specify the names of the sex chromosomes (as it no longer makes the assumption that they are called chrX/chrY/X/Y). If left empty, all chromosomes are treated as diploid.

Please let me know if the new version solves the problems

[HuangYihang1222]

Thank you for the quick fix on this pipeline! It works well overall, but I encountered a few minor issues:

1. When creating the environment, the channel settings cause the environment creation to fail. I had to manually add the channel in the .yml file.
2. Minimap2 in the environment requires a manual upgrade to function properly. You can do this with conda install -c bioconda minimap2=2.26.
3. I encountered an issue with bcftools:

   bcftools --version  
   bcftools: error while loading shared libraries: libgsl.so.25: cannot open shared object file: No such file or directory  
   

   This can be fixed by forcibly installing gsl=2.5, which provides libgsl.so.25.

Other than these issues, everything else seems to be working fine. Thanks again!