New upstream version 0+git20180126.5c6fed5+dfsg

# VarMatch

robust matching of small variant datasets using flexible scoring schemes

# [New] VarMatch for multiple datasets

The GIAB branch of VarMatch can compare complex variants from multiple (>=3) datasets.

To download the GIAB branch, use command ```git clone -b giab https://github.com/medvedevgroup/varmatch.git```

# Authors

- Chen Sun (The Pennsylvania State University)

- Paul Medvedev (The Pennsylvania State University)

# Release Date

### TBA

Any questions about VarMatch, please email to chensun at cse dot psu dot edu.

If you identify a bug in VarMatch, please either reported on 'github Issues' of VarMatch, or email directly to chensun at cse dot psu dot edu.

# Prerequisite

- GCC 4.7 or later for c++11 support

- Python 2.7 or later

> *matplotlib is not a prerequisite if either `-f`, `-G` or `-C` parameter is used

# Installation

**Quick Install Instruction:**

### Quick Install Instruction:

You can build VarMatch from source. 

```

git clone https://github.com/medvedevgroup/varmatch.git

make all

```

### Uninstall

`cd` to the directory of VarMatch

```

make clean

```

# Test Data Set

- Links to a test data set (~15M) : [https://github.com/medvedevgroup/varmatch/blob/master/test_data.txt](https://github.com/medvedevgroup/varmatch/blob/master/test_data.txt)

- Links to the data used for bencharmking in the paper: [https://github.com/medvedevgroup/varmatch/blob/master/data.txt](https://github.com/medvedevgroup/varmatch/blob/master/data.txt)

# Usage

### Quick Usage:

*compare two vcf files to match variants*

```

./varmatch -b baseline.vcf -q query.vcf -g ref.fa -o out -f

./varmatch -b baseline.vcf -q query.vcf -g ref.fa -o output -f

```

- `-b` baseline vcf file

- `-q` query vcf file

- `-o` output file prefix, default value is `out`

- `-f` fast mode*, equivalent to use parameters `-u 0 -m 0 -s 0 -C`

>*fast mode is suggested for ordinary analysis

>*fast mode is suggested for ordinary analysis.

>VarMatch accept baseline and query in VCF file format (e.g. xx.vcf), it does not accept gz file (e.g. xx.vcf.gz) in current version.

>see [Results of VarMatch](#results-of-varmatch) section for intepretation of results in output directory. 

### Detail Usage

  `-G`, `--no_graph`

        disable graphic module

  `-C`, `--disable_curves`

        disable Precision-Recall curves, if use -G or --no_graph, then

        automatically disable these curves

  `-f`, `--fast_mode`

        In this mode, automatically disable graphic module and precision-

        recall curves, only performs one matching criterion.

use `-h/--help` for detailed help message.

# Results of VarMatch

### varmatch (.match) file

You can find varmatch (.match) files in VarMatch output directory, filename is in the format of query`x`.`u`\_`m`\_`s`.match

- `x` is the id of queries.

- `u` is the value of parameter `-u`, `--score_unit`

- `m` is the value of parameter `-m`, `--match_mode`

- `s` is the value of parameter `-s`, `--score_scheme`

For instance, if you use one query VCF file and use `-f` parameter, there is query1.0_0_0.match in your output file.

varmatch file contains the information of matched VCF entries from baseline and query VCF file.

Lines in varmatch file started with `#` are comment lines. They contain general information of baseline and query VCF file, and also general information of varmatch file. The first 2 lines of .match file starts with `###`:

- `###VCF1` is the baseline VCF filename

- `###VCF2` is the query VCF filename

varmatch files contains at least 9 fields:

  1. `CHROM` field represents chromosome ID

  2. `POS` field represents genome position on reference genome sequence

  3. `REF` field represents reference allele

  4. `ALT` field represents alternative alleles, multiple alleles are separated by `/`

  5. `VCF1` field represents variants from baseline. If it is a direct match, this column is `.`. If it is not a direct match, this column contains variants separated by `;`. Each variant contains three information: reference genome position, reference allele, alternative alleles.

  6. `VCF2` field represents variants from query. If it is a direct match, this column is `.`.

  7. `PHASE1` field represents phasing information of variants from baseline. If it is a direct match, this column is `.`.

  8. `PHASE2` field represents phasing information of variants from query. If it is a direct match, this column is `.`.

  9. `SCORE` field represents the total score of variants from baseline and query, calculated based on given score unit, and score scheme.

The meaning of each line in varmatch file:

> variants in `VCF1` column is equivalent to variants in `VCF2` column. If applying them separately on `REF` sequence, which is a substring of reference genome sequence starts at position `POS` of chromosome `CHROM`, can get the same donor sequences in `ALT` column.

The phasing information of variants in `VCF1` and `VCF2` are separately in `PHASE1` and `PHASE2`. 

varmatch (.match) file format gives a standard representation of equivalent variants, especially for complex variants.

If you have any suggestions of improving .match file format, please [contact me](#contact).

### .stat file

It contains some statistical information.

# License

VarMatch is now under Apache2 license!

See [license.txt](https://github.com/medvedevgroup/varmatch/blob/master/license.txt)

# Contact

chensun@cse.psu.edu

You also can report bugs or suggest features using issue tracker at GitHub [https://github.com/medvedevgroup/varmatch](https://github.com/medvedevgroup/varmatch)

# Acknowledgements

If using VarMatch, please cite:

Sun, C., & Medvedev, P. (2016). VarMatch: robust matching of small variant datasets using flexible scoring schemes. Bioinformatics, 33(9), 1301-1308.

Corresponding BiBTex:

```

@article{sun2016varmatch,

  title={VarMatch: robust matching of small variant datasets using flexible scoring schemes},

  author={Sun, Chen and Medvedev, Paul},

  journal={Bioinformatics},

  volume={33},

  number={9},

  pages={1301--1308},

  year={2016},

  publisher={Oxford University Press}

}

```

This project has been supported in part by NSF awards DBI-1356529, CCF-1439057, IIS-1453527, and IIS-1421908.

This file contains raw data needed for benchmarking in the paper.

# Reference Genome Sequence GRCh37

> You can download GRCh37 reference genome sequence in either of the following ftp sites:

### ftp://ftp.ensembl.org/pub/release-75/fasta/homo_sapiens/dna/

### ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/BUILD.37.3/

> The reference genome sequence GRCh37 is stored in these ftp sites by chromosomes. If you want to perform whole genome analysis, a single integrated FASTA file is needed.

# CHM1 Data Set:

### Bowtie2-Freebayes

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/CHM1.bt2.fb.vcf.gz

### Bowtie2-HC

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/CHM1.bt2.hc.vcf.gz

# NA12878 Data Set:

### GIAB benchmark (Version 2.18)

ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/analysis/NIST_union_callsets_06172013/NISTIntegratedCalls_14datasets_131103_allcall_UGHapMerge_HetHomVarPASS_VQSRv2.18_all_nouncert_excludesimplerep_excludesegdups_excludedecoy_excludeRepSeqSTRs_noCNVs.vcf.gz

### If previous link is disabled, go to the following directory for GIAB benchmark history version:

ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/analysis/NIST_union_callsets_06172013/

### Bowtie2-Freebayes

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/NA12878.bt2.fb.vcf.gz

### Bowtie2-UG

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/NA12878.bt2.ug.vcf.gz

### Bowtie2-SAMtools

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/NA12878.bt2.st.vcf.gz

### Bowtie2-Platypus

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/NA12878.bt2.pt.vcf.gz

### BWA_MEM-Freebayes

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20140312_chm1_alignment_heng/vcfs/NA12878.mem.fb.vcf.gz

 No newline at end of file

import numpy as np

import matplotlib

matplotlib.use('Agg')

import matplotlib.pyplot as plt

t = np.arange(0., 5., 0.2)

plt.plot(t,t,'r-o')

plt.plot(t,t/2, 'g-^')

plt.savefig('xx.png')

1   /home/varmatch/human/chr1.fa

2   /home/varmatch/human/chr2.fa

17  /home/varmatch/human/backup/chr17.fa

X   /home/varmatch/human/chrxx.fa

Y   /home/anotherpath/human/chrY/human.y.fa

 No newline at end of file

#!/usr/bin/env python

from sys import argv

import argparse

import math

import scipy.stats as stats

citation = 'Please cite our paper'

parser = argparse.ArgumentParser(epilog=citation)

parser.add_argument('--qu', metavar='N', help='quality number(QUAL) threshold >= N (default: N=30)', default=30)

parser.add_argument('--ab', metavar='N', help='allele balance(AB) threshold <= N%% (default: N=20)', default=20)

parser.add_argument('--fs', metavar='N', help='Fisher strand P-vale <= N (default: N=0.001)', default=0.001)

parser.add_argument('--rd', metavar='N', default=65,

                    help="average read depth=N, maximum read depth(MD) threshold >= N+4*sqrt(N) (default: N=65),"

                         " use --rd 0 to disable MD filter")

parser.add_argument('-i', metavar='input.vcf', help='input VCF file')

parser.add_argument('-o', metavar='output.vcf', help='output VCF file name(default: output.vcf)', default='output.vcf')

parser.add_argument('--homo', action='store_true', help='filter out homozygous variants')

parser.add_argument('--nf', action='store_true', help="no filters used in Heng Li review")

parser.add_argument('--snp', action='store_true', help="only want SNPs")

parser.add_argument('--indel', action='store_true', help='only want INDELs')

args = parser.parse_args()

def main():

    if len(argv) < 2:

        parser.print_help()

        exit()

    filter_homo = args.homo

    if not filter_homo:

        print ('Warning: compulsively filter out homozygous variants :)')

        filter_homo = True

    md = 0  # maximum depth filter

    if args.rd != 0:

        md = args.rd + 4 * math.sqrt(args.rd)

    else:

        print ('Warning: maximum depth(MD) filter is disabled because read depth = 0')

    output_file = open(args.o, 'w')

    with open(args.i) as input_file:

        for line in input_file.readlines():

            qu_fail = False

            ab_fail = False

            fs_fail = False

            md_fail = False

            if line.startswith('#'):

                output_file.write(line)

                continue

            columns = line.split('\t')

            if len(columns) < 8:

                print ('Warning: current variant does not contains enough info for filtering')

                continue

            ab_contain = False

            ab_pass = True

            two_alleles = False

            pv = 1.0

            rd = -1

            srf = -1

            srr = -1

            saf_list = []

            sar_list = []

            alt = columns[4]

            if ',' in alt:

                two_alleles = True

            ref = columns[3]

            is_indel = False

            for a in alt.split(','):

                if len(ref) != len(a):

                    is_indel = True

            if args.snp and is_indel:

                continue

            if args.indel and not is_indel:

                continue

            # Filter out homozygous

            if filter_homo:

                if len(columns) < 10:

                    print('Warning: variant does not contain enough info to filter homozygous variants')

                format_col = columns[8].split(':')

                gt_index = -1

                for i in range(len(format_col)):

                    if format_col[i] == 'GT':

                        gt_index = i

                if gt_index == -1:

                    print ('Warning: variant does not contain genotype info')

                    continue

                val_col = columns[9].split(':')

                gt_val = val_col[gt_index]

                gt_col = []

                if '/' in gt_val:

                    gt_col = gt_val.split('/')

                elif '|' in gt_val:

                    gt_col = gt_val.split('|')

                else:

                    print ('Warning: unrecognized genotype info')

                    continue

                if gt_col[0] == gt_col[1]:

                    continue

            if args.nf:

                output_file.write(line)

                continue

            quality_num = float(columns[5])

            # quality filter(QU)

            if quality_num < args.qu:

                qu_fail = True

            if not qu_fail:

                output_file.write(line)

                continue

            info_col = columns[7].split(';')

            for info in info_col:

                val_col = info.split('=')

                info_name = val_col[0]

                info_val = val_col[1]

                if info_name == 'AB':

                    ab_contain = True

                    if two_alleles:

                        ab_col = info_val.split(',')

                        for ab in ab_col:

                            if float(ab) > args.ab * 0.01:

                                ab_pass = False

                    else:

                        if float(info_val) > args.ab * 0.01:

                            ab_pass = False

                elif info_name == 'DP':

                    rd = int(info_val)

                elif info_name == 'SRF':

                    srf = int(info_val)

                elif info_name == 'SRR':

                    srr = int(info_val)

                elif info_name == 'SAF':

                    if two_alleles:

                        temp_list = info_val.split(',')

                        saf_list = [int(temp_list[0]), int(temp_list[1])]

                    else:

                        saf_list = [int(info_val)]

                elif info_name == 'SAR':

                    if two_alleles:

                        temp_list = info_val.split(',')

                        sar_list = [int(temp_list[0]), int(temp_list[1])]

                    else:

                        sar_list = [int(info_val)]

            # AB filter

            if not ab_contain or not ab_pass:

                ab_fail = True

            if not ab_fail:

                output_file.write(line)

                continue

            # Maximum depth(MD) filter

            if rd == -1:

                print ('Warning: current variant does not contain read depth info')

                continue

            elif rd < md:

                md_fail = True

            if not md_fail:

                output_file.write(line)

                continue

            # Fisher strand filter(FS)

            oddsratio, pv = stats.fisher_exact([[srf, srr], [saf_list[0], sar_list[0]]])

            if pv > args.fs:

                fs_fail = True

            if two_alleles:

                oddsratio, pv = stats.fisher_exact([[srf, srr], [saf_list[1], sar_list[1]]])

                if pv > args.fs:

                    fs_fail = True

            if not fs_fail:

                output_file.write(line)

                continue

    output_file.close()

if __name__ == "__main__":

    main()

 No newline at end of file