<commit_before>/*
 *  Groovy NGS Utils - Some simple utilites for processing Next Generation Sequencing data.
 *
 *  Copyright (C) 2018 <NAME>, <EMAIL><at><EMAIL>
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
package gngs.tools

import gngs.*
import groovy.transform.CompileStatic
import groovy.transform.Memoized
import groovy.util.logging.Log
import htsjdk.variant.variantcontext.Genotype
import htsjdk.variant.variantcontext.GenotypesContext
import htsjdk.variant.variantcontext.VariantContext

import static gngs.Sex.*

/**
 * A utility class to contain some functionality for computing allele numbers / counts
 * 
 * @author <NAME>
 */
@CompileStatic
class AlleleNumber {
   
    static boolean isX(String chr) {
        chr.endsWith('X')
    }
    
    static boolean isY(String chr) {
        chr.endsWith('Y')
    }
    
    static int getAlleleNumber(gngs.Sex sex, String chr) {
        
        if(sex == MALE) {
            if(isX(chr) || isY(chr))
                return 1
            else
                return 2
        }
        else
        if(sex == FEMALE) {
            if(isY(chr))
                return 0
            else
                return 2
        }
        else
            assert false : "Only MALE and FEMALE sexes are supported"
    }
}

/**
 * A tool to create a "population statistics VCF" that represents the frequency of
 * variants from a given set of samples by entries in the info field.
 * <p>
 * The following population metrics are calculated:
 * <li>AN   -   the total number of possible alleles that could have been observed at the site
 * <li>AC   -   the total number of alleles that were observed at the site
 * <li>GTC  -   the total number of individuals that were sampled to generate AN and AC
 * <p>
 * The current implementation infers the sexes of the samples so they do not need to be specified.
 * However this does rely on having sufficient coverage of the autosomes and X chromosome to determine
 * sex. If this is not true, an error will be reported.
 * <p> 
 * <b>Limitations:</b>
 * The current implementation has a number of limitations:
 * 
 * <li>Only one entry is genererated for each unique site in the input VCFs. Therefore this will
 *      <i>incorrect</i> when there is more than one allele observed at a site (multi-allelic sites).
 * <li>If a site appears twice in one of the source VCFs it may be counted incorrectly (eg: multi-allelic sites that are 
 *     decomposed)
 * <li>Only VCF information is considered, not gVCF, which is to say, an allele is counted as 
 *     unobserved at a site even if there was no coverage and therefore was not ascertained at that site.
 * 
 * These limitations should be corrected to generate truly accurate counts.
 * 
 * @author <NAME>
 */
@Log
class CreatePopulationStatisticsVCF extends ToolBase {
    
    int numSamples
    
    /**
     * The inferred sexes as a list
     */
    private List<gngs.Sex> sexes
    
    /**
     * The inferred sexes of the samples in the VCFs
     */
    Map<String,gngs.Sex> sampleSexes = [:]
    
    PrintWriter out = null

    @Override
    public void run() {
        
        if(opts.arguments().isEmpty()) {
            parser.usage()
            System.exit(0)
        }
        
        if(out != null) {
            // Use it as is
        }
        else
        if(opts.o) {
            out = new PrintWriter(new File(opts.o))
        }
        else {
            out = new PrintWriter(System.out)
        }
        
        printVCFHeader()
        
        inferSexes()
        
        // Since we end up with a sex inferred for every sample ...
        numSamples = sampleSexes.size()
        
        VariantContext lastVariant = null
        ProgressCounter progress = new ProgressCounter(withRate: true, withTime:true, extra: {
            "chr: $lastVariant.contig pos: $lastVariant.start (${lastVariant.genotypes[0].sampleName})"
        })
        VCFSiteWalker locusWalker = new VCFSiteWalker(opts.arguments())
        locusWalker.walk { List<VariantContext> variants ->
            VariantContext v0 = variants[0]
            lastVariant = v0
            int ac = computeAlleleCount(variants)
            int an = computeAlleleNumber(v0.contig)
            int gtc = variants.size()*numSamples // todo: check this is right
            printVCFSite(v0, ac, an, gtc)
            progress.count()
        }
        progress.end()
    }
    
    @CompileStatic
    @Memoized
    int computeAlleleNumber(String contig) {
       (int)sexes.sum { gngs.Sex sex ->  AlleleNumber.getAlleleNumber(sex, contig) } <commit_msg>add flush to fix truncated VCF output<commit_after>/*
 *  Groovy NGS Utils - Some simple utilites for processing Next Generation Sequencing data.
 *
 *  Copyright (C) 2018 <NAME>, ssadedin<at>gmail.com
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
package gngs.tools

import gngs.*
import groovy.transform.CompileStatic
import groovy.transform.Memoized
import groovy.util.logging.Log
import htsjdk.variant.variantcontext.Genotype
import htsjdk.variant.variantcontext.GenotypesContext
import htsjdk.variant.variantcontext.VariantContext

import static gngs.Sex.*

/**
 * A utility class to contain some functionality for computing allele numbers / counts
 * 
 * @author <NAME>
 */
@CompileStatic
class AlleleNumber {
   
    static boolean isX(String chr) {
        chr.endsWith('X')
    }
    
    static boolean isY(String chr) {
        chr.endsWith('Y')
    }
    
    static int getAlleleNumber(gngs.Sex sex, String chr) {
        
        if(sex == MALE) {
            if(isX(chr) || isY(chr))
                return 1
            else
                return 2
        }
        else
        if(sex == FEMALE) {
            if(isY(chr))
                return 0
            else
                return 2
        }
        else
            assert false : "Only MALE and FEMALE sexes are supported"
    }
}

/**
 * A tool to create a "population statistics VCF" that represents the frequency of
 * variants from a given set of samples by entries in the info field.
 * <p>
 * The following population metrics are calculated:
 * <li>AN   -   the total number of possible alleles that could have been observed at the site
 * <li>AC   -   the total number of alleles that were observed at the site
 * <li>GTC  -   the total number of individuals that were sampled to generate AN and AC
 * <p>
 * The current implementation infers the sexes of the samples so they do not need to be specified.
 * However this does rely on having sufficient coverage of the autosomes and X chromosome to determine
 * sex. If this is not true, an error will be reported.
 * <p> 
 * <b>Limitations:</b>
 * The current implementation has a number of limitations:
 * 
 * <li>Only one entry is genererated for each unique site in the input VCFs. Therefore this will
 *      <i>incorrect</i> when there is more than one allele observed at a site (multi-allelic sites).
 * <li>If a site appears twice in one of the source VCFs it may be counted incorrectly (eg: multi-allelic sites that are 
 *     decomposed)
 * <li>Only VCF information is considered, not gVCF, which is to say, an allele is counted as 
 *     unobserved at a site even if there was no coverage and therefore was not ascertained at that site.
 * 
 * These limitations should be corrected to generate truly accurate counts.
 * 
 * @author <NAME>
 */
@Log
class CreatePopulationStatisticsVCF extends ToolBase {
    
    int numSamples
    
    /**
     * The inferred sexes as a list
     */
    private List<gngs.Sex> sexes
    
    /**
     * The inferred sexes of the samples in the VCFs
     */
    Map<String,gngs.Sex> sampleSexes = [:]
    
    PrintWriter out = null

    @Override
    public void run() {
        
        if(opts.arguments().isEmpty()) {
            parser.usage()
            System.exit(0)
        }
        
        if(out != null) {
            // Use it as is
        }
        else
        if(opts.o) {
            out = new PrintWriter(new File(opts.o))
        }
        else {
            out = new PrintWriter(System.out)
        }
        
        printVCFHeader()
        
        out.flush()
        
        inferSexes()
        
        // Since we end up with a sex inferred for every sample ...
        numSamples = sampleSexes.size()
        
        VariantContext lastVariant = null
        ProgressCounter progress = new ProgressCounter(withRate: true, withTime:true, extra: {
            "chr: $lastVariant.contig pos: $lastVariant.start (${lastVariant.genotypes[0].sampleName})"
        })
        VCFSiteWalker locusWalker = new VCFSiteWalker(opts.arguments())
        locusWalker.walk { List<VariantContext> variants ->
            VariantContext v0 = variants[0]
            lastVariant = v0
            int ac = computeAlleleCount(variants)
            int an = computeAlleleNumber(v0.contig)
            int gtc = variants.size()*numSamples // todo: check this is right
            printVCFSite(v0, ac, an, gtc)
            progress.count()
        }
        progress.end()
    }
    
    @CompileStatic
    @Memoized
    int computeAlleleNumber(String contig) {
       (int)sexes.sum { gngs.Sex sex ->  AlleleNumber.getAlleleNumber(sex, contig) } 
{'fraction_non_alphanumeric': 0.07641757416447616, 'fraction_numerical': 0.005257228689447991, 'mean_word_length': 2.3426419830561658, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4554650', 'n_tokens_mistral': 3088, 'n_tokens_neox': 2916, 'n_words': 1139}