Metadata-Version: 1.2
Name: mitty
Version: 2.23.1
Summary: Simulator for genomic data
Home-page: https://github.com/sbg/Mitty
Author: Seven Bridges Genomics
Author-email: kaushik.ghose@sbgdinc.com
License: UNKNOWN
Download-URL: https://github.com/sbg/Mitty/archive/development.zip
Description: Mitty is a data simulator meant to help debug aligners and variant
        callers.
        
        It requires Python 3.4 or later. It is released under the
        `Apache <LICENSE.txt>`__ license.
        
        |DOI|
        
        .. figure:: docs/images/genome-simulation.png?raw=true
           :alt: Genome simulation
        
           Genome simulation
        
        .. figure:: docs/images/read-simulation.png?raw=true
           :alt: Read simulation
        
           Read simulation
        
        .. figure:: docs/images/reports.png?raw=true
           :alt: Report generation
        
           Report generation
        
        Features
        ========
        
        -  Generate reads from the whole genome, a single tiny region or from a
           set of regions as desired
        -  Handles X, Y chromosomes and polyploidy IF the VCF GT field is
           properly set
        -  Read qname stores correct POS, CIGAR and the sizes of variants
           covered by the read
        -  Name of sample included in read allowing us to mix FASTQs from
           different simulations/sources
        
           -  Can mix viral contamination into reads
           -  Can do tumor/normal mixes
        
        -  Corruption module adds sequencing errors to reads
        
           -  Read models can be sampled from existing BAM files
        
        -  "God aligner" writes out a BAM with perfect alignments which can be
           used for BAM comparisons
        -  Simple genome simulator to generate VCFs with SNPs and different
           sizes of Insertions and Deletions for aligner/caller testing
        
        It is also informative to browse the `release
        notes <release_notes.txt>`__ and poke around the detailed documentation
        under the `/docs <docs/>`__ folder
        
        Quickstart
        ==========
        
        Install
        -------
        
        ::
        
            conda create -n mymitty python=3.5
        
        -  Using a virtual env is recommended, but it's a personal choice
        -  Mitty requires Python3 to run
        
        Install Mitty from the public github repository. **Requires pip 9.0 or
        later**
        
        ::
        
            pip install --upgrade setuptools pip
            pip install git+https://github.com/sbg/Mitty.git
        
        Run tests
        ---------
        
        ::
        
            nosetests mitty.test -v
        
        Program help
        ------------
        
        Help is available from the command line. Simply invoking the base
        program with no arguments will list all the sub programs with some short
        help:
        
        ::
        
            mitty
        
        Detailed help on particular commands is also available:
        
        ::
        
            mitty generate-reads --help
        
        Running commands with the verbose option allows you to tune what
        messages (ranging from Errors to Debug) you get.
        
        ::
        
            mitty -v{1,2,3,4} <command>
        
        Detailed tutorial with commentary
        =================================
        
        Example scripts and data
        ------------------------
        
        If you want to follow along with this tutorial you will find the
        relevant example scripts and test data in the
        https://github.com/kghosesbg/mitty-demo-data project.
        
        *The separate project has been created to avoid making the main source
        tree bulky by filling it with binary files not needed for program
        operation.*
        
        (Some of the examples require other tools such as ``samtools``, ``bwa``
        and ``vcftools``)
        
        Generating reads
        ----------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/generating-reads/gen-reads.sh>`__)
        
        Aliases used
        ~~~~~~~~~~~~
        
        ::
        
            FASTA=../data/human_g1k_v37.fa.gz
            SAMPLEVCF=../data/1kg.20.22.vcf.gz
            SAMPLENAME=HG00119
            REGION_BED=region.bed
            FILTVCF=HG00119-filt.vcf.gz
            READMODEL=1kg-pcr-free.pkl
            COVERAGE=30
            READ_GEN_SEED=7
            FASTQ_PREFIX=HG00119-reads
            READ_CORRUPT_SEED=8
        
        Prepare VCF file for taking reads from
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        The read simulator can not properly generate ground truth reads from
        overlapping variants (e.g. a deletion that overlaps a SNP on the same
        chromosome copy), complex variants (variant calls where the REF and ALT
        are both larger than 1 bp) or any entry that does not precisely lay out
        the ALT sequence in the ALT column, such as angle bracket ID references
        and variants breakend notation. Such calls must be filtered from a VCF
        file before it can be used as a sample to generate reads from.
        
        ::
        
            mitty -v4 filter-variants \
              ${SAMPLEVCF} \
              ${SAMPLENAME} \
              ${REGION_BED} \
              - \
              2> vcf-filter.log | bgzip -c > ${FILTVCF}
        
            tabix -p vcf ${FILTVCF}
        
        The ``region.bed`` file is used to select out just the parts of the VCF
        we would like to generate reads from. The input VCF can be a population
        VCF (say from the 1000G project). The ``${SAMPLENAME}`` indicates which
        sample to extract from the VCF.
        
        ``vcf-filter.log`` looks like::
        
        ::
        
            $ cat vcf-filter.log 
            DEBUG:mitty.lib.vcfio:Starting filtering ...
            DEBUG:mitty.lib.vcfio:Filtering ('1', 1000000, 2000000)
            DEBUG:mitty.lib.vcfio:Complex variant 1:1827835 CTT -> ('CT', 'C')
            DEBUG:mitty.lib.vcfio:Filtering ('3', 60830384, 61830384)
            DEBUG:mitty.lib.vcfio:Complex variant 3:60835995 TTCTCTCTCTCTCTCTCTCTCTC -> ('TTCTCTCTCTCTCTCTCTCTCTCTCTCTC', 'T')
            DEBUG:mitty.lib.vcfio:Complex variant 3:60897726 TAC -> ('TACAC', 'T')
            DEBUG:mitty.lib.vcfio:Complex variant 3:60970457 GGTGTGT -> ('GGT', 'G')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61205001 ATG -> ('ATGTG', 'A')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61309628 GGTGTGTGTGTGTGTGT -> ('GGTGTGTGTGTGTGTGTGTGT', 'G')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61360782 AGTGTGTGTGT -> ('AGTGTGTGTGTGT', 'A')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61469726 CAAAAAAAAAAAAAAA -> ('CA', 'C')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61488707 TTGTG -> ('TTG', 'T')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61509647 TCACACACACA -> ('TCACACACACACACA', 'T')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61522251 TACACAC -> ('TACAC', 'T')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61541525 CA -> ('CAAA', 'C')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61633465 CAAAAAAA -> ('CAAAAAAAAAAAAA', 'C')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61718383 AAAATAAAT -> ('AAAAT', 'A')
            DEBUG:mitty.lib.vcfio:Complex variant 3:61731724 CTTT -> ('CTT', 'C')
            DEBUG:mitty.lib.vcfio:Processed 2807 variants
            DEBUG:mitty.lib.vcfio:Sample had 2807 variants
            DEBUG:mitty.lib.vcfio:Discarded 15 variants
            DEBUG:mitty.lib.vcfio:Took 0.28023505210876465 s
        
        **NOTE: If the BED file is not sorted, the output file needs to be
        sorted again.**
        
        ``${FILTVCF}`` can now be used to generate reads and will serve as a
        truth VCF for VCF benchmarking.
        
        Listing and inspecting read models
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        ``mitty list-read-models`` will list built in read models
        
        ``mitty list-read-models -d ./mydir`` will list additional custom models
        stored in ``mydir``
        
        ``mitty describe-read-model 1kg-pcr-free.pkl model.png`` prints a panel
        of plots showing read characteristics
        
        See later for a list of read models supplied with Mitty and their
        characteristics
        
        Creating custom read models
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        **This is only needed if none of the existing read models match your
        requirements**
        
        Prepare a Illumina type read model from a BAM
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/read-models/sampled-model.sh>`__)
        
        ::
        
            BAM=../data/sample.bam
            MODELNAME=sampled-model.pkl
        
            mitty -v4 create-read-model bam2illumina \
              --every 10 \
              --min-mq 20 \
              -t 2 \
              --max-bp 300 \
              --max-tlen 1000 \
              ${BAM} \
              ${MODELNAME} \
              'Sampled model created for the demo'
        
            mitty describe-read-model ${MODELNAME} ${MODELNAME}.png
        
        Create arbitrary Illumina type read models
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        
        The read model file is just a Python pickle file of a dictionary
        carrying specifications for the model. You can create arbitrary models
        by specifying your own parameters. Please see the `read model
        documentation <docs/readmodelformat.md>`__ for a description of all the
        parameters.
        
        Easily prepare completely synthetic Illumina type read model
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/read-models/synthetic-model.sh>`__)
        
        Mitty also supplies a model generator (``synth-illumina``) to generate
        custom Illumina like read models with template sizes and base quality
        patterns following simple mathematical distributions. This model
        generator allows us to quickly create reads models with a wide variety
        of independently variable parameters.
        
        ::
        
            MODELNAME=synthetic-model.pkl
        
            mitty create-read-model synth-illumina \
              ${MODELNAME} \
              --read-length 121 \
              --mean-template-length 400 \
              --std-template-length 20 \
              --bq0 30 \
              --k 200 \
              --sigma 5 \
              --comment 'Model created for the demo' \
        
        The generated read model looks like:
        
        .. figure:: docs/images/synthetic-model.pkl.png?raw=true
           :alt: 
        
        Generating perfect reads
        ~~~~~~~~~~~~~~~~~~~~~~~~
        
        ::
        
            mitty -v4 generate-reads \
              ${FASTA} \
              ${FILTVCF} \
              ${SAMPLENAME} \
              ${REGION_BED} \
              ${READMODEL} \
              ${COVERAGE} \
              ${READ_GEN_SEED} \
              >(gzip > ${FASTQ_PREFIX}1.fq.gz) \
               ${FASTQ_PREFIX}-lq.txt \
               --fastq2 >(gzip > ${FASTQ_PREFIX}2.fq.gz) \
               --threads 2
        
        *(When you supply a model file name to the read generator it will first
        look among the builtin read models to see if the file name is a match
        (typically these are in the ``mitty/data/readmodels`` folder). It will
        then treat the model file name as a path and try and load that from your
        file system - which is the case in this particular example.)*
        
        The produced FASTQs have qnames encoding the correct read alignments for
        each read in the template. qnames may exceed the SAM specification limit
        (nominally 254 characters, but there
        `are <https://github.com/samtools/htslib/issues/520>`__
        `caveats <https://github.com/pysam-developers/pysam/issues/447>`__). In
        such cases the qname in the FASTQ is truncated and the complete qname is
        printed in the side-car file ``${FASTQ_PREFIX}-lq.txt``.
        
        The qname format can be obtained by executing ``mitty qname``
        
        Reference reads
        ^^^^^^^^^^^^^^^
        
        As you might expect, by passing a VCF with no variants we can generate
        reads with no variants, representing the reference genome. One Mitty
        feature to be aware of is ploidy inference from the VCF: if there are no
        variants in a contig, Mitty assumes that contig is diploid, otherwise
        Mitty uses the GT (genotype) tag to infer ploidy.
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/reference-reads/ref-reads-m.sh>`__)
        
        Hence for a human male we set up the VCF (``human-m-ref.vcf``) as:
        
        ::
        
            ##fileformat=VCFv4.2
            ##FILTER=<ID=PASS,Description="All filters passed">
            ##fileDate=20160824
            ##contig=<ID=1,length=249250621>
            ##contig=<ID=2,length=243199373>
            ...
            ##contig=<ID=X,length=155270560>
            ##contig=<ID=Y,length=59373566>
            ##contig=<ID=MT,length=16569>
            ##FORMAT=<ID=GT,Number=1,Type=String,Description="Consensus Genotype across all datasets with called genotype">
            #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ref
            X       1       .       A       G       50      PASS    .       GT      0
            Y       1       .       A       G       50      PASS    .       GT      0
        
        When we generate reads from this VCF say from chromosome 1, X and Y we
        see the following trace:
        
        ::
        
            DEBUG:mitty.lib.vcfio:Contig: 1, ploidy: 2 (Assumed. Contig was empty)
            DEBUG:mitty.lib.vcfio:Contig: X, ploidy: 1 
            DEBUG:mitty.lib.vcfio:Contig: Y, ploidy: 1 
        
        which tells us that contig 1 has been assumed to be diploid, whereas X
        and Y are inferred to be haploid because of the dummy entries we set.
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/reference-reads/ref-reads-f.sh>`__)
        
        Technically the human female VCF could be set up as a VCF with just the
        header (Mitty would infer all contigs to be diploid), but it turns out
        that some tools (including pysam) operate incorrectly when the VCF is
        completely empty, so we supply one dummy line for the VCF
        (``human-f-ref.vcf``) as:
        
        ::
        
            ##fileformat=VCFv4.2
            ##FILTER=<ID=PASS,Description="All filters passed">
            ##fileDate=20160824
            ##contig=<ID=1,length=249250621>
            ##contig=<ID=2,length=243199373>
            ...
            ##contig=<ID=X,length=155270560>
            ##contig=<ID=Y,length=59373566>
            ##contig=<ID=MT,length=16569>
            ##FORMAT=<ID=GT,Number=1,Type=String,Description="Consensus Genotype across all datasets with called genotype">
            #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ref
            X       1       .       A       G       50      PASS    .       GT      0/0
        
        Truncating reads
        ^^^^^^^^^^^^^^^^
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/generating-reads/truncated-reads.sh>`__)
        
        For some experiments you might want to generate custom sized reads.
        ``generate-reads`` allows you to do this with the ``--truncate-to``
        argument
        
        ::
        
            FASTQ_PREFIX=HG00119-truncated-reads
            mitty -v4 generate-reads \
              ${FASTA} \
              ${FILTVCF} \
              ${SAMPLENAME} \
              ${REGION_BED} \
              ${READMODEL} \
              ${COVERAGE} \
              ${READ_GEN_SEED} \
              >(gzip > ${FASTQ_PREFIX}1.fq.gz) \
              ${FASTQ_PREFIX}-lq.txt \
              --fastq2 >(gzip > ${FASTQ_PREFIX}2.fq.gz) \
              --truncate-to 60 \
              --threads 2
        
        This generates the same kind of reads as before, but all the reads are
        60bp long, instead of their usual length. You can not make reads longer
        than what the model originally specifies (This to ensure that the read
        corruption code will work seamlessly with such truncated reads.)
        
        Un-pairing reads
        ^^^^^^^^^^^^^^^^
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/generating-reads/unpaired-reads.sh>`__)
        
        For some experiments you might want to use the existing Illumina or
        other model that normally produces paired-end reads to generate
        single-end reads instead. The ``--unpair`` argument allows you to do
        this. Note that in this case you should not pass in a second output
        FASTQ file.
        
        ::
        
            FASTQ_PREFIX=HG00119-unpaired-reads
            mitty -v4 generate-reads \
              ${FASTA} \
              ${FILTVCF} \
              ${SAMPLENAME} \
              ${REGION_BED} \
              ${READMODEL} \
              ${COVERAGE} \
              ${READ_GEN_SEED} \
              >(gzip > ${FASTQ_PREFIX}.fq.gz) \
               ${FASTQ_PREFIX}-lq.txt \
               --unpair \
               --threads 2
        
        Corrupting reads
        ~~~~~~~~~~~~~~~~
        
        The reads generated using the previous command have no base call errors.
        Base call errors can be introduced into the reads using the following
        command.
        
        ::
        
            mitty -v4 corrupt-reads \
              ${READMODEL} \
              ${FASTQ_PREFIX}1.fq.gz >(gzip > ${FASTQ_PREFIX}-corrupt1.fq.gz) \
              ${FASTQ_PREFIX}-lq.txt \
              ${FASTQ_PREFIX}-corrupt-lq.txt \
              ${READ_CORRUPT_SEED} \
              --fastq2-in ${FASTQ_PREFIX}2.fq.gz \
              --fastq2-out >(gzip > ${FASTQ_PREFIX}-corrupt2.fq.gz) \
              --threads 2
        
        *Using a read model different to that used to generate the reads
        originally can lead to undefined behavior, including crashes.*
        
        As mentioned, the side-car file ``${FASTQ_PREFIX}-lq.txt`` carries
        qnames longer than 254 characters. The output side-car file
        ``${FASTQ_PREFIX}-corrupt-lq.txt`` similarly carries longer qnames from
        the corrupted reads. The qname for each corrupted template is identical
        to the original, uncorrupted template, except for the addition of an
        MD-like tag that allows recovery of the original bases before sequencing
        errors were introduced.
        
        The BQ profile of the sample FASTQ generated by this command (generated
        by FASTQC) looks like
        
        Mate 1: |FASTQC screenshot showing BQ distribution|
        
        Mate 2: |FASTQC screenshot showing BQ distribution|
        
        can be compared with the empirical model profile shown in the Appendix
        for the ``1kg-pcr-free.pkl`` model.
        
        After passing these reads through an aligner and viewing them on a
        genome browser, such as IGV one can make a quick inspection of the
        alignments.
        
        .. figure:: docs/images/igv-alignment-qname.png?raw=true
           :alt: IGV screenshot showing read qname and one het variant
        
           IGV screenshot showing read qname and one het variant
        
        Since the qname carries the correct alignment and CIGAR string you can
        match that against the actual alignment and CIGAR string for spot
        checks.
        
        Perfect BAM (God aligner)
        ~~~~~~~~~~~~~~~~~~~~~~~~~
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/god-aligner/god-aligner.sh>`__)
        
        Passing the simulated FASTQ through the god aligner produces a "perfect
        BAM" which can be used as a truth BAM for comparing alignments from
        different aligners. This truth BAM can also be used to test variant
        callers by removing one moving part (the aligner) from the analysis
        chain.
        
        ::
        
            FASTA=../data/human_g1k_v37.fa.gz
            FASTQ_PREFIX=../generating-reads/HG00119-reads
            GODBAM=HG00119-god.bam
        
            mitty -v4 god-aligner \
              ${FASTA} \
              ${FASTQ_PREFIX}-corrupt1.fq.gz \
              ${FASTQ_PREFIX}-corrupt-lq.txt \
              ${GODBAM} \
              --fastq2 ${FASTQ_PREFIX}-corrupt2.fq.gz \
              --threads 2
        
        Analysis
        --------
        
        Mitty supplies some tools to help with benchmarking and debugging of
        aligner/caller pipelines.
        
        Alignment accuracy
        ------------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/alignment-accuracy/alignment-accuracy.sh>`__)
        
        *(Assumes bwa and samtools are installed)*
        
        ::
        
            BAM=hg001-bwa.bam
        
            bwa mem \
              ${FASTA} \
              ${FASTQ_PREFIX}-corrupt1.fq.gz \
              ${FASTQ_PREFIX}-corrupt2.fq.gz | samtools view -bSho temp.bam
            samtools sort temp.bam > ${BAM}
            samtools index ${BAM}
        
            mitty -v4 debug alignment-analysis process\
              ${BAM} \
              ${FASTQ_PREFIX}-corrupt-lq.txt \
              ${BAM}.alignment.npy \
              --fig-prefix ${BAM}.alignment \
              --max-d 200 \
              --max-size 50 \
              --plot-bin-size 10
        
        This invocation will process ``${BAM}`` and summarize the alignment
        performance in a numpy data file (``${BAM}.alignment.npy``). It will
        plot them in a set of figures named with the prefix
        ``${BAM}.alignment``. The alignment error will be assessed upto a
        maximum of 200bp. The program will check variants from 50bp deletions to
        50bp insertions, putting them into 10bp size bins. SNPs are always
        counted and placed in their own spearate bin.
        
        .. figure:: docs/images/aligner-report-example-1.png?raw=true
           :alt: MQ plots
        
           MQ plots
        
        .. figure:: docs/images/aligner-report-example-2.png?raw=true
           :alt: Alignment accuracy plots
        
           Alignment accuracy plots
        
        Subset a BAM for detailed analysis
        ----------------------------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/subset-bam/subset-bam.sh>`__)
        
        The ``subset-bam`` debug subtool allows us to select out reads from a
        BAM based on whether they contain variants and whether they fall within
        a certain d\_err range
        
        Ex 1: Extract only reads from SNPs
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        ::
        
            BAMOUT=HG00119-bwa-snps.bam
            mitty -v4 debug subset-bam \
              ${BAMIN} \
              ${FASTQ_PREFIX}-corrupt-lq.txt \
              ${BAMOUT} \
              --v-range 0 0 \
              --reject-reference-reads \
              --processes 2
        
        .. figure:: docs/images/igv-subsetbam-snp-example.png?raw=true
           :alt: Reads under SNPs
        
           Reads with SNPs only
        
        Ex 2: Extract only correctly aligned reads from deletions
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        ::
        
            BAMOUT=HG00119-bwa-good-del.bam
            mitty -v4 debug subset-bam \
              ${BAMIN} \
              ${FASTQ_PREFIX}-corrupt-lq.txt \
              ${BAMOUT} \
              --v-range -10000 -1 \
              --d-range -5 5 \
              --reject-reference-reads \
              --processes 2
        
        Ex 3: Extract only reads from deletions mis-aligned by more than 5 bp
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        ::
        
            BAMOUT=HG00119-bwa-poor-del.bam
            mitty -v4 debug subset-bam \
              ${BAMIN} \
              ${FASTQ_PREFIX}-corrupt-lq.txt \
              ${BAMOUT} \
              --v-range -10000 -1 \
              --d-range -5 5 \
              --reject-d-range \
              --reject-reference-reads \
              --processes 2
        
        .. figure:: docs/images/igv-subsetbam-del-example.png?raw=true
           :alt: Reads undel DELs
        
           Reads undel DELs
        
        Variant calling accuracy, parametrized by variant size
        ------------------------------------------------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/variant-call-analysis/call-analysis.sh>`__)
        
        We can use a set of tools developed by the GA4GH consortium to compare a
        VCF produced by a pipeline with a truth VCF. One of the outputs of the
        comparator tools is a VCF (called an evaluation VCF) where each call is
        annotated with by whether it is a TP, FN, FP or GT error.
        
        ``variant-call-analysis`` is a program that summarizes the data in such
        an evaluation VCF in terms of variant size.
        
        ::
        
            EVCF=../data/0.9.29.eval.vcf.gz
            CSV=0.9.29.eval.data.csv
            FIG=caller-report-example.png
        
            mitty -v4 debug variant-call-analysis process \
             ${EVCF} \
             ${CSV} \
             --max-size 75 \
             --fig-file ${FIG} \
             --plot-bin-size 5 \
             --title 'Example call analysis plot'
        
        This invocation will process ``0.9.29.eval.vcf.gz`` produced by vcfeval,
        write the results as a comma separated file (``0.9.29.eval.data.csv``)
        and then plot them in ``caller-report-example.png``.
        
        .. figure:: docs/images/caller-report-example.png?raw=true
           :alt: P/R plots
        
           P/R plots
        
        To replot already processed data use the ``plot`` subcommand instead of
        the ``process`` subcommand
        
        ::
        
            mitty -v4 debug variant-call-analysis plot \
              ${CSV} \
              caller-report-example2.png \
              --plot-bin-size 10 \
              --plot-range 50 \
              --title 'Example call analysis plot'
        
        Improvements, regressions in variant calling
        --------------------------------------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/call-fate/call-fate.sh>`__)
        
        When comparing different versions of a pipeline, or two different
        pipelines the Precision and Recall curves and summary tables give some
        information about the improvements and regressions introduced, but only
        at a very coarse level. The ``call-fate`` tool compares two evaluation
        VCFs and tracks the transitions of variant calls between different call
        categories (TP, FN, GT, FP) which allow us to see in greater detail the
        improvements and regressions going from one pipeline to the other.
        
        In the examples below we are comparing two evaluation VCF files
        ``{0.9.29, 0.9.32}.eval.vcf.gz``
        
        ::
        
            EVCF1=../data/0.9.29.eval.vcf.gz
            EVCF2=../data/0.9.32.eval.vcf.gz
            OUTPREFIX=fate-29-32
        
            mitty -v4 debug call-fate \
             ${EVCF1} \
             ${EVCF2} \
             - \
             ${OUTPREFIX}-summary.txt | vcf-sort | bgzip -c > ${OUTPREFIX}.vcf.gz
        
        (This assumes we have vcf-tools available so we can sort the VCF)
        
        The program produces a summary table output:
        
        ::
        
            Improved    SNP INDEL
            --------------------
            FN->TP:     391 587
            FN->GT:     35  134
            GT->TP:     342 429
            FP->N:      23710   8680
        
            Unchanged   SNP INDEL
            --------------------
            TP->TP:     3504751 680531
            FN->FN:     4571    12259
            GT->GT:     1851    12639
            FP->FP:     379274  211967
        
            Regressed   SNP INDEL
            --------------------
            TP->FN:     351 488
            TP->GT:     134 305
            GT->FN:     27  129
            N->FP:      13674   7180
        
        And a VCF file with 12 samples, corresponding to the 12 categories
        above. For each variant the GT field is 0/0 for all samples except the
        one corresponding to the transition category it belongs to. This allows
        us to easily visualize the fate of individual variants using, for
        example, IGV.
        
        .. figure:: docs/images/call-fate-igv.png?raw=true
           :alt: call-fate VCF overlay on IGV
        
           call-fate VCF overlay on IGV
        
        Set differences of two or more BAM files derived from the same FASTQ(s)
        -----------------------------------------------------------------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/partition-bams/partition-bams.sh>`__)
        
        One way of making a detailed examination of the effects of changes to
        alignment algorithms is to track how read alignments from the same FASTQ
        change. The ``partition-bams`` subtool allows us to take 2 or more BAMs
        and apply a membership criterion to each read and classify the reads
        according to how they fared in each of the BAMs.
        
        In the `associated
        example <https://github.com/kghosesbg/mitty-demo-data/blob/master/partition-bams/partition-bams.sh>`__
        ``bwa mem`` is run with three different values of the ``-r`` parameter
        on a small FASTQ. We then apply the membership criterion \|d\_err\| < 10
        and analyze the three BAMs.
        
        ::
        
            FASTQ_PREFIX=../generating-reads/HG00119-reads
        
            mitty -v4 debug partition-bams \
              myderr \
              d_err \
              --threshold 10 \
              --sidecar_in ${FASTQ_PREFIX}-corrupt-lq.txt \
              --bam HG00119-bwa1.bam \
              --bam HG00119-bwa2.bam \
              --bam HG00119-bwa3.bam
        
        This tool produces a summary file ``myderr_summary.txt`` that looks
        like:
        
        ::
        
            (A)(B)(C)   576
            (A)(B)C 0
            (A)B(C) 0
            (A)BC   57
            A(B)(C) 84
            A(B)C   0
            AB(C)   0
            ABC 359793
        
        In this nomenclature A is the set and (A) is the complement of this set.
        The set labels A, B, C ... (upto a maximum of 10) refer to the BAM files
        in sequence they were supplied.
        
        Thus, ABC means all the reads which have a \|d\_err\| < 10 in all the
        three files. AB(C) means all the reads which have a \|d\_err\| < 10 in A
        and B but not C, and so on. A reader familiar with Venn diagrams is
        referred to the chart below for a translation of the three dimensional
        case to a three way Venn diagram. Higher dimensions are harder to
        visualize as Venn diagrams.
        
        .. figure:: docs/images/sets.png?raw=true
           :alt: Sets to Venn diagram
        
           Sets to Venn diagram
        
        The tool also produces a set of files following the naming convention:
        
        ::
        
            myderr_(A)(B)(C)_A.bam
            myderr_(A)(B)(C)_B.bam
            myderr_(A)(B)(C)_C.bam
            myderr_(A)(B)C_A.bam
            myderr_(A)(B)C_B.bam
            myderr_(A)(B)C_C.bam
            ...
        
        The first part of the name follows the convention outlined above. The
        trailing A, B, C refer to the orginal source BAM of the reads. So
        ``myderr_(A)(B)(C)_B.bam`` carries reads from bam B that have \|d\_err\|
        >= 10 in all the three BAMs.
        
        An example of throwing these files up on a genome browser and inspecting
        them is given below
        
        .. figure:: docs/images/igv-sets.png?raw=true
           :alt: BAM partitions on IGV
        
           IGV Bam Partitions
        
        The criteria the ``partition-bam`` tool can be run on can be obtained by
        passing it the ``--criteria`` option.
        
        Generating samples (genomes)
        ----------------------------
        
        Mitty also has features to generate simulated genomes in the form of VCF
        files.
        
        Simulated variants
        ~~~~~~~~~~~~~~~~~~
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/simulating-variants/simulate-variants.sh>`__)
        
        The ``simulate-variants`` command generates a VCF with simulated
        variants. The program carries three basic models for variant simulation
        - SNPs, insertions and deletions and is invoked as follows:
        
        ::
        
            FASTA=../data/human_g1k_v37.fa.gz
            SAMPLENAME=S0
            BED=region.bed
            VCF=sim.vcf.gz
        
            mitty -v4 simulate-variants \
              - \
              ${FASTA} \
              ${SAMPLENAME} \
              ${BED} \
              7 \
              --p-het 0.6 \
              --model SNP 0.001 1 1 \
              --model INS 0.0001 10 100 \
              --model DEL 0.0001 10 100 | bgzip -c > ${VCF}
        
            tabix ${VCF}
        
        The model parameters are given by
        
        .. raw:: html
        
           <P>
        
         refers to the variant model to use
        
        .. raw:: html
        
           <P>
        
        is the probability of a variant being placed on any given base indicate
        the size ranges of the variants produced. These are ignored for SNP
        
        This VCF should be run through the ``filter-variants`` program as usual
        before taking reads. This is especially important because the simulation
        can produce illegaly overlapping variants which will be taken out by
        this step.
        
        Invoking ``mitty simulate-variants --list-models`` will list available
        models
        
        Miscellaneous utilities
        -----------------------
        
        Bam to truth
        ------------
        
        (`Example
        script <https://github.com/kghosesbg/mitty-demo-data/blob/master/bam-to-truth/bam_to_truth.sh>`__)
        
        Sometimes we want to treat the alignment from one aligner (e.g. BWA) as
        the truth and then check how other aligners do relative to that. An
        ideal tool would do a read by read comparison, and we have some other
        tools that do this, however such comparisons, because they need to
        matchup read qnames, can become expensive. This is a compromise method.
        
        ``bam-to-truth`` creates FASTQ file(s) from a BAM file, changing the
        qname to encode the alignment of the read. The FASTQ files can then be
        used like any other simulated FASTQ, to analyze alignment performance
        for other aligners relative to the original aligner. The code only
        writes out reads for which both mates are mapped and for which both
        mates have MQ greater than the supplied threshold.
        
        Variant size distribution
        -------------------------
        
        Plot variant size distribution in VCF file:
        
        ::
        
            mitty -v4 debug variant-by-size \
              hg001.vcf.gz hg001.variant.size.csv \
              --plot-bin-size 5 \
              --max-size 100 \
              --title "HG001" \
              --fig-file hg001.variant.png 
        
        Appendix
        ========
        
        Qname format
        ------------
        
        Read alignment and simulation metadata are stored in the qname in the
        following format.
        
        ::
        
            $ mitty qname
        
                  @index|sn|chrom:copy|strand:pos:cigar:v1,v2,...:MD|strand:pos:cigar:v1,v2,...:MD*
        
            index:  unique code for each template. If you must know, the index is simply a monotonic counter
                    in base-36 notation
            sn:     sample name. Useful if simulation is an ad-mixture of sample + contaminants or multiple samples
            chrom:  chromosome id of chromosome the read is taken from
            copy:   copy of chromosome the read is taken from (0, 1, 2, ...) depends on ploidy
            strand: 0: forward strand, 1: reverse strand
            pos:    Position of first (left-most) reference matching base in read (Same definition as for POS in BAM)
                    One based
                    For reads coming from completely inside an insertion this is, however, the POS for the insertion
            cigar:  CIGAR string for correct alignment.
                    For reads coming from completely inside an insertion, however, the CIGAR string is:
                    '>p+nI' where:
                       '>' is the unique key that indicates a read inside a long insertion
                       'p' is how many bases into the insertion branch the read starts
                       'n' is simply the length of the read
            v1,v2,..: Comma separated list of variant sizes covered by this read
                        0 = SNP
                        + = INS
                        - = DEL
            MD:     Read corruption MD tag. This is empty for perfect reads, but is filled with an MD formatted
                    string for corrupted reads. The MD string is referenced to the original, perfect read, not
                    a reference sequence.
        
            Notes:
        
            - The alignment information is repeated for every read in the template.
              The example shows what a PE template would look like.
            - The pos value is are one based to make comparing qname info in genome browser easier
            - qnames longer than N characters are stored in a side-car file alongside the simulated FASTQs.
              The qname in the FASTQ file itself is truncated to N characters.
        
              Nominally, N=254 according to the SAM spec, but due to bugs in some versions of htslib
              this has been set shorter to 240.
              A truncated qname is detected when the last character is not *
        
        Example from a perfect read:
        
        ``@8|INTEGRATION|1:1|1:1930067:27=1X63=1I158=:0,1:|0:1929910:184=1X63=1I1=:0,1:*``
        
        Corresponding read after passing through read corruption model:
        
        ``@8|INTEGRATION|1:1|1:1930067:27=1X63=1I158=:0,1:126T29T16G15G28C0T2T8C9T0C6T|0:1929910:184=1X63=1I1=:0,1:9G25G3T1C36T87G10G0C0A4A1G4C0A0G0G0A1A0T0C3A1C0A0T0C0C0T0G2T0A0A1A0G1G0T0G0A0A0A0C1T1A0T0C0T0C0T0A1T1A0A1A0T0A0C0A0A*``
        
        ::
        
            @8           - simulated reads serial number
            INTEGRATION  - name of sample (from VCF) reads are generated from
            1:1          - read is from chromosome 1, copy 1 (copy numbers start 0, 1, ...)
            1:1930067:27=1X63=1I158=:0,1:126T29T16G15G28C0T2T8C9T0C6T 
                         - read one metadata with fields as:
        
                1:1930067       - read is from reverse strand (1, as opposed to 0 - see the metdata for the mate) at pos 1930067
                27=1X63=1I158=  - this is the correct CIGAR for this read
                0,1             - this read carries two variants a SNP (0) and an insertion (1)
                126T29T16G15G28C0T2T8C9T0C6T  
                                - This string follows the conventions for an MD tag and indicates the sequencing error
                                  relative to the uncorrupted read
        
            0:1929910:184=1X63=1I1=:0,1:9G25G3T1C36T87G10G0C0A4A1G4C0A0G0G0A1A0T0C3A1C0A0T0C0C0T0G2T0A0A1A0G1G0T0G0A0A0A0C1T1A0T0C0T0C0T0A1T1A0A1A0T0A0C0A0A 
                         - read two metadata in same format
        
        Built-in read models
        --------------------
        
        .. figure:: docs/images/1kg-pcr-free.png?raw=true
           :alt: 
        
        .. figure:: docs/images/hiseq-X-v2.5-Garvan.png?raw=true
           :alt: 
        
        .. figure:: docs/images/old-Garvan.png?raw=true
           :alt: 
        
        .. figure:: docs/images/hiseq-2500-v1-pcr-free.png?raw=true
           :alt: 
        
        .. figure:: docs/images/hiseq-X-v1-HLI.png?raw=true
           :alt: 
        
        .. |DOI| image:: https://zenodo.org/badge/65392390.svg
           :target: https://zenodo.org/badge/latestdoi/65392390
        .. |FASTQC screenshot showing BQ distribution| image:: docs/images/1kg-pcr-free-corrupt-fastqc-r1.png?raw=true
        .. |FASTQC screenshot showing BQ distribution| image:: docs/images/1kg-pcr-free-corrupt-fastqc-r2.png?raw=true
        
Keywords: simulator,genomics,ngs,read mapper,aligner,variant caller
Platform: UNKNOWN
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3.4
Requires-Python: >=3.4
