<commit_before>package org.broadinstitute.sting.playground.gatk.walkers;
//import org.broadinstitute.sting.gatk.iterators.LocusIterator;
import org.broadinstitute.sting.gatk.LocusContext;
import org.broadinstitute.sting.gatk.refdata.ReferenceOrderedDatum;
import org.broadinstitute.sting.gatk.refdata.rodDbSNP;
import org.broadinstitute.sting.gatk.walkers.BasicLociWalker;
import org.broadinstitute.sting.playground.utils.AlleleFrequencyEstimate;
import net.sf.samtools.SAMRecord;

import java.util.List;
import java.util.Arrays;

public class AlleleFrequencyWalker extends BasicLociWalker<AlleleFrequencyEstimate, Integer> {

    int N=2;

    public AlleleFrequencyEstimate map(List<ReferenceOrderedDatum> rodData, char ref, LocusContext context) 
    {
        // Convert context data into bases and 4-base quals
        // Convert context data into bases and 4-base quals
        String bases = getBases(context);
        double quals[][] = getOneBaseQuals(context);

        // Count bases
        int[] base_counts = new int[4];
        for (byte b : bases.getBytes())
            base_counts[nuc2num[b]]++;

        // Find alternate allele - 2nd most frequent non-ref allele
        // (maybe we should check for ties and eval both or check most common including quality scores)
        int altnum = -1; // start with first base numerical identity (0,1,2,3)
        double altcount = -1; // start with first base count
        int refnum = nuc2num[ref];

        for (int b=0; b<4; b++) {
            if ((b != refnum) && (base_counts[b] > altcount)) {
                altnum = b; altcount = base_counts[b];
            }
        }
        assert(altnum != -1);

        AlleleFrequencyEstimate alleleFreq = AlleleFrequencyEstimator(context.getLocation().toString(), N, bases.getBytes(), quals, refnum, altnum, bases.length());

        alleleFreq.notes = String.format("A:%d C:%d G:%d T:%d", 
                                            base_counts[nuc2num['A']],
                                            base_counts[nuc2num['C']],
                                            base_counts[nuc2num['G']],
                                            base_counts[nuc2num['T']]);

        // Print dbSNP data if its there
        if (true) {
            for ( ReferenceOrderedDatum datum : rodData ) {
                if ( datum != null && datum instanceof rodDbSNP) {
                    rodDbSNP dbsnp = (rodDbSNP)datum;
                    //System.out.printf("  DBSNP %s on %s => %s%n", dbsnp.toSimpleString(), dbsnp.strand, Utils.join("/", dbsnp.getAllelesFWD()));
                    alleleFreq.notes += String.format(" ROD: %s ",dbsnp.toMediumString());
                }
            }
        }
        return alleleFreq;
    }

    static public String getBases (LocusContext context) {
        // Convert bases to CharArray
        int numReads = context.getReads().size(); //numReads();
        //byte[] bases = new byte[numReads];
        String base_string = "";
        //int refnum = nuc2num[ref];

        List<SAMRecord> reads = context.getReads();
        List<Integer> offsets = context.getOffsets();
        for (int i =0; i < numReads; i++ ) {
            SAMRecord read = reads.get(i);
            int offset = offsets.get(i);

            //bases[i] = read.getReadBases()[offset];
            base_string += read.getReadString().charAt(offset);
        }
        return base_string;
    }

    static public double[][] getOneBaseQuals (LocusContext context) {
        int numReads = context.getReads().size(); //numReads();
        double[][] quals = new double[numReads][4];

        List<SAMRecord> reads = context.getReads();
        List<Integer> offsets = context.getOffsets();
        for (int i =0; i < numReads; i++ ) {
            SAMRecord read = reads.get(i);
            int offset = offsets.get(i);
            int callednum = nuc2num[read.getReadBases()[offset]];
            quals[i][callednum] = 1 - Math.pow(10.0, (double)read.getBaseQualities()[offset] / -10.0);

            // Set all nonref qual scores to their share of the remaining probality not "used" by the reference base's qual
            double nonref_quals = (1.0 - quals[i][callednum]) / 3;
            for (int b=0; b<4; b++)
                if (b != callednum)
                    quals[i][b] = nonref_quals;
        }

        return quals;
    }

    public AlleleFrequencyEstimate AlleleFrequencyEstimator(String location, int N, byte[] bases, double[][] quals, int refnum, int altnum, int depth)
    {

        // q = hypothetical %nonref
        // qstar = true %nonref
        // G = N, qstar, alt
        // N = number of chromosomes
        // alt = alternate hypothesis base (most frequent nonref base)
        // ref = reference base

        // b = number of bases at locus

        double epsilon = 0; //  1e-2;
        double qstar;
        int qstar_N;

        double qstep = 0.001;
        double qend = 1.0 + qstep / 10; // makes sure we get to 1.0 even with rounding error of qsteps accumulating

        double best_qstar     = Math.log10(0);
        double best_qhat      = Math.log10(0);<commit_msg>AlleleFrequencyWalker now can parse 4-base probs<commit_after>package org.broadinstitute.sting.playground.gatk.walkers;
//import org.broadinstitute.sting.gatk.iterators.LocusIterator;
import org.broadinstitute.sting.gatk.LocusContext;
import org.broadinstitute.sting.gatk.refdata.ReferenceOrderedDatum;
import org.broadinstitute.sting.gatk.refdata.rodDbSNP;
import org.broadinstitute.sting.gatk.walkers.BasicLociWalker;
import org.broadinstitute.sting.playground.utils.AlleleFrequencyEstimate;
import net.sf.samtools.SAMRecord;

import java.util.List;

public class AlleleFrequencyWalker extends BasicLociWalker<AlleleFrequencyEstimate, Integer> {

    int N=2;

    public AlleleFrequencyEstimate map(List<ReferenceOrderedDatum> rodData, char ref, LocusContext context) 
    {
        // Convert context data into bases and 4-base quals
        // Convert context data into bases and 4-base quals
        String bases = getBases(context);
        double quals[][] = getOneBaseQuals(context);

        // Count bases
        int[] base_counts = new int[4];
        for (byte b : bases.getBytes())
            base_counts[nuc2num[b]]++;

        // Find alternate allele - 2nd most frequent non-ref allele
        // (maybe we should check for ties and eval both or check most common including quality scores)
        int altnum = -1; // start with first base numerical identity (0,1,2,3)
        double altcount = -1; // start with first base count
        int refnum = nuc2num[ref];

        for (int b=0; b<4; b++) {
            if ((b != refnum) && (base_counts[b] > altcount)) {
                altnum = b; altcount = base_counts[b];
            }
        }
        assert(altnum != -1);

        AlleleFrequencyEstimate alleleFreq = AlleleFrequencyEstimator(context.getLocation().toString(), N, bases.getBytes(), quals, refnum, altnum, bases.length());

        alleleFreq.notes = String.format("A:%d C:%d G:%d T:%d", 
                                            base_counts[nuc2num['A']],
                                            base_counts[nuc2num['C']],
                                            base_counts[nuc2num['G']],
                                            base_counts[nuc2num['T']]);

        // Print dbSNP data if its there
        if (true) {
            for ( ReferenceOrderedDatum datum : rodData ) {
                if ( datum != null && datum instanceof rodDbSNP) {
                    rodDbSNP dbsnp = (rodDbSNP)datum;
                    //System.out.printf("  DBSNP %s on %s => %s%n", dbsnp.toSimpleString(), dbsnp.strand, Utils.join("/", dbsnp.getAllelesFWD()));
                    alleleFreq.notes += String.format(" ROD: %s ",dbsnp.toMediumString());
                }
            }
        }
        return alleleFreq;
    }

    static public String getBases (LocusContext context) {
        // Convert bases to CharArray
        int numReads = context.getReads().size(); //numReads();
        //byte[] bases = new byte[numReads];
        String base_string = "";
        //int refnum = nuc2num[ref];

        List<SAMRecord> reads = context.getReads();
        List<Integer> offsets = context.getOffsets();
        for (int i =0; i < numReads; i++ ) {
            SAMRecord read = reads.get(i);
            int offset = offsets.get(i);

            //bases[i] = read.getReadBases()[offset];
            base_string += read.getReadString().charAt(offset);
        }
        return base_string;
    }

    static public double[][] getOneBaseQuals (LocusContext context) {
        int numReads = context.getReads().size(); //numReads();
        double[][] quals = new double[numReads][4];

        List<SAMRecord> reads = context.getReads();
        List<Integer> offsets = context.getOffsets();
        for (int i =0; i < numReads; i++ ) {
            SAMRecord read = reads.get(i);
            int offset = offsets.get(i);

            // First, set the called base to the correct quality
            int callednum = nuc2num[read.getReadBases()[offset]]; // number of the called base (0,1,2,3)
            quals[i][callednum] = 1 - Math.pow(10.0, (double)read.getBaseQualities()[offset] / -10.0);

            // For other 3 bases, check if 2-base probs exist
            assert (reads.size() > 0);
            Object SQ_field = reads.get(0).getAttribute("SQ");
            if (SQ_field == null) {
                // Set all nonref qual scores to their share of the remaining probality not "used" by the reference base's qual
                double nonref_quals = (1.0 - quals[i][callednum]) / 3;
                for (int b=0; b<4; b++)
                    if (b != callednum)
                        quals[i][b] = nonref_quals;
            }else{
                assert (SQ_field instanceof byte[]);
                byte[] hex_quals = (byte[]) SQ_field;
                System.out.printf("SQ field (hex): %s\n", bytesToHexString(hex_quals));
                System.out.printf("SAM record: %s\n", read.format());

                int hex_qual = hex_quals[offset];
                int called2num = hex_qual & 0x3;
                double qual2 = (double)(hex_qual >> 2) / 100.0;
                System.out.printf("2ND %x %d %f\n", hex_qual, called2num, qual2);
                quals[i][called2num] = qual2;

                // 
                double nonref_quals = (1.0 - quals[i][callednum] - quals[i][called2num]) / 3;
                for (int b=0; b<4; b++)
                    if (b != callednum && b != called2num)
                        quals[i][b] = nonref_quals;
            }
        }

        //print_base_qual_matrix(quals);
        return quals;
    }

    // Code pulled from SAMUtils for debugging 
    static String bytesToHexString(final byte[] data) {
        final char[] chars = new char[2 * data.length];
        for (int i = 0; i < data.length; i++) {
            final byte b = data[i];
            chars[2*i] = toHexDigit((b >> 4) & 0xF);
            chars[2*i+1] = toHexDigit(b & 0xF);
        }
        return new String(chars);
    }

    private static char toHexDigit(final int value) {
        return (char) ((value < 10) ? ('0' + value) : ('A' + value - 10));
    }

    public AlleleFrequencyEstimate AlleleFrequencyEstimator(String location, int N, byte[] bases, double[][] quals, int refnum, int altnum, int depth)
    {

        // q = hypothetical %nonref
        // qstar = true %nonref
        // G = N, qstar, alt
        // N = number of chromosomes
        // alt = alternate hypothesis base (most frequent nonref base)
        // ref = reference base

        // b = number of bases at locus

        double epsilon = 0; //  1e-2;
        double qstar;
        int qstar_N;

        double qstep = 0.001;
        double qend = 1.0 + qstep / 10; // makes sure we get to 1.0 even with rounding error of qsteps accumulating

        double best_qstar     = Math.log10(0);
        double best_qhat      = Math.log10(0);
{'fraction_non_alphanumeric': 0.12042885638297872, 'fraction_numerical': 0.012882313829787235, 'mean_word_length': 2.0135236664162286, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4827621', 'n_tokens_mistral': 3809, 'n_tokens_neox': 3579, 'n_words': 932}