search

kevlar-dev / kevlar

Reference-free variant discovery in large eukaryotic genomes
https://kevlar.readthedocs.io
MIT License
40 stars 9 forks source link

> 30,000 variants reported from trio-analysis with FPR 0.001? #386

Open moldach opened 3 years ago

moldach commented 3 years ago

I'm seeing a very high number of reported variants from kevlar. Even after restricting the max_fpr values to 0.001 I'm seeing > 30,000 variants?

config.json

{
    "ksize": 31,
    "recountmem": "950G",
    "numsplit": 16,
    "samples": {
        "casemin": 6,
        "ctrlmax": 1,
        "case": {
            "fastx": [
                "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018P_R1.fastq.gz",
                "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018P_R2.fastq.gz"
            ],
            "memory": "950G",
            "label": "Proband",
            "max_fpr": 0.001
        },
        "controls": [
            {
                "fastx": [
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018M_R1.fastq.gz",
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018M_R2.fastq.gz"
                ],
                "memory": "950G",
                "label": "Mother",
                "max_fpr": 0.001
            },
            {
                "fastx": [
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018F_R1.fastq.gz",
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018F_R2.fastq.gz"
                ],
                "memory": "950G",
                "label": "Father",
                "max_fpr": 0.001
            }
        ],
        "coverage": {
            "mean": 30.0,
            "stdev": 10.0
        }
    },
    "mask": {
        "fastx": [
            "/tiered/mtgraovac/indexes/GRCh37/Homo_sapiens.GRCh37.dna.toplevel.fa"
        ],
        "memory": "100G",
        "max_fpr": 0.005
    },
    "reference": {
        "fasta": "/tiered/mtgraovac/indexes/GRCh37/Homo_sapiens.GRCh37.dna.toplevel.fa",
        "memory": "100G",
        "max_fpr": 0.025
    },
    "localize": {
        "seedsize": 51,
        "delta": 50,
        "seqpattern": ".",
        "maxdiff": 10000
    },
    "varfilter": null
}

Counting lines in output

## snakemake submission
$ snakemake --snakefile Snakefile --configfile config.json --cores 32 --directory ./ -p calls

## counting lines in the vcf
$ bgzip -d calls.scored.sorted.vcf.gz
$ wc -l calls.scored.sorted.vcf 
31561 calls.scored.sorted.vcf
standage commented 3 years ago

Hi @moldach. How many of those calls have a PASS value in the FILTER column? If that number is still high, my first guess is that many of these calls are going to be associated with common variants or repetitive regions such as segmental duplications. I'd suggest providing a BED file with the locations of common variants and segdups to the varfilter config value so that associated variant calls can be marked appropriately. See https://kevlar.readthedocs.io/en/latest/tutorial.html for more info.

Given the amount of memory you're using for each sample, I'm wondering what your coverage is. I tested mostly trios with 30x-50x average coverage. You're using much more memory than I ever did, and if that is due to high coverage you might need to fiddle with the case/control abundance thresholds a bit as well.

I hope this helps!

moldach commented 3 years ago

Didn't think of that! I think two of the GATK callers I tried before reported everything in the VCFand you had to filter for PASS but not with the tools I'm using most often so I forgot about checking that.

I tried using bcftools to get that information from kevlar output but I get an error:

$ bcftools view -i 'ID=PASS' calls.scored.sorted.vcf > output.vcf

[E::bcf_hdr_parse_line] Could not parse the header line: "##INFO=<ID=IKMERS,Number=1,Type=Integer,Description="number of "interesting" (novel) k-mers spanning the variant alternate allele">"
[W::bcf_hdr_parse] Could not parse header line: ##INFO=<ID=IKMERS,Number=1,Type=Integer,Description="number of "interesting" (novel) k-mers spanning the variant alternate allele">
No such INFO field: PASS
standage commented 3 years ago

I can't say that I'm familiar with bcftools, but the FILTER column is the last of the seven "core" columns in VCF, not one of the subsequent variable number of INFO columns described in the headers. For a less-sophisticated approach, I'd suggest the following shell one-liners.

$ # Just count the variants passing filters
$ grep -c $'\tPASS\t' calls.scored.sorted.vcf
$
$ # Actually retrieve the passing variants
$ grep $'\tPASS\t' calls.scored.sorted.vcf > output.vcf
moldach commented 3 years ago

Okay thank you.

But, riddle me this, kevlar run with max_fpr of 0.05 and max_fpr of 0.001 are returning exactly the same number of putative variants when filtered for PASS. Why would this be?

max_fpr=0.05

{
    "ksize": 31,
    "recountmem": "250G",
    "numsplit": 16,
    "samples": {
        "casemin": 6,
        "ctrlmax": 1,
        "case": {
            "fastx": [
                "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018P_R1.fastq.gz",
                "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018P_R2.fastq.gz"
            ],
            "memory": "250G",
            "label": "Proband",
            "max_fpr": 0.05
        },
        "controls": [
            {
                "fastx": [
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018M_R1.fastq.gz",
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018M_R2.fastq.gz"
                ],
                "memory": "250G",
                "label": "Mother",
                "max_fpr": 0.05
            },
            {
                "fastx": [
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018F_R1.fastq.gz",
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018F_R2.fastq.gz"
                ],
                "memory": "250G",
                "label": "Father",
                "max_fpr": 0.05
            }
        ],
        "coverage": {
            "mean": 30.0,
            "stdev": 10.0
        }
    },
    "mask": {
        "fastx": [
            "/tiered/mtgraovac/indexes/GRCh37/Homo_sapiens.GRCh37.dna.toplevel.fa"
        ],
        "memory": "40G",
        "max_fpr": 0.005
    },
    "reference": {
        "fasta": "/tiered/mtgraovac/indexes/GRCh37/Homo_sapiens.GRCh37.dna.toplevel.fa",
        "memory": "80G",
        "max_fpr": 0.025
    },
    "localize": {
        "seedsize": 51,
        "delta": 50,
        "seqpattern": ".",
        "maxdiff": 10000
    },
    "varfilter": null
}
$ grep -c $'\tPASS\t' calls.scored.sorted.vcf
3233

max_fpr=0.001

{
    "ksize": 31,
    "recountmem": "950G",
    "numsplit": 16,
    "samples": {
        "casemin": 6,
        "ctrlmax": 1,
        "case": {
            "fastx": [
                "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018P_R1.fastq.gz",
                "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018P_R2.fastq.gz"
            ],
            "memory": "950G",
            "label": "Proband",
            "max_fpr": 0.001
        },
        "controls": [
            {
                "fastx": [
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018M_R1.fastq.gz",
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018M_R2.fastq.gz"
                ],
                "memory": "950G",
                "label": "Mother",
                "max_fpr": 0.001
            },
            {
                "fastx": [
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018F_R1.fastq.gz",
                    "/tiered/mtgraovac/CombGene/CG00018/FASTQ/CG00018F_R2.fastq.gz"
                ],
                "memory": "950G",
                "label": "Father",
                "max_fpr": 0.001
            }
        ],
        "coverage": {
            "mean": 30.0,
            "stdev": 10.0
        }
    },
    "mask": {
        "fastx": [
            "/tiered/mtgraovac/indexes/GRCh37/Homo_sapiens.GRCh37.dna.toplevel.fa"
        ],
        "memory": "100G",
        "max_fpr": 0.005
    },
    "reference": {
        "fasta": "/tiered/mtgraovac/indexes/GRCh37/Homo_sapiens.GRCh37.dna.toplevel.fa",
        "memory": "100G",
        "max_fpr": 0.025
    },
    "localize": {
        "seedsize": 51,
        "delta": 50,
        "seqpattern": ".",
        "maxdiff": 10000
    },
    "varfilter": null
}
$ grep -c $'\tPASS\t' calls.scored.sorted.vcf
3233
standage commented 3 years ago

Most likely explanation: while the lower FPR offers better theoretical accuracy, both FPRs offer identical effective accuracy. A more thorough investigation of the 3k PASSing variants may reveal another more interesting story though.

standage commented 3 years ago

To be clear, I meant another additional story, not another alternative story.