src/main/java/org/ncgr/gwas/SNPReadCounter.java
public static void main(String[] args) throws Exception {

        Options options = new Options();
        
        Option geneFileOption = new Option("f", "geneFile", true, "Gene list (e.g. generated by SNPComparer)");
        geneFileOption.setRequired(true);
        options.addOption(geneFileOption);

        Option gffFileOption = new Option("g", "gffFile", true, "GFF file");
        gffFileOption.setRequired(true);
        options.addOption(gffFileOption);

        Option vcfFileOption = new Option("v", "vcfFile", true, "VCF file");
        vcfFileOption.setRequired(true);
        options.addOption(vcfFileOption);

        Option altMinOption = new Option("a", "altMin", true, "minimum ALT depth for SNP to be counted");
        altMinOption.setRequired(true);
        options.addOption(altMinOption);

        Option altFracMinOption = new Option("afm", "altFracMin", true, "minimum ALT fraction for a gene to be output");
        altFracMinOption.setRequired(true);
        options.addOption(altFracMinOption);

        CommandLineParser parser = new DefaultParser();
        HelpFormatter formatter = new HelpFormatter();
        CommandLine cmd;
        try {
            cmd = parser.parse(options, args);
        } catch (ParseException e) {
            System.out.println(e.getMessage());
            formatter.printHelp("SNPReadCounter", options);
            System.exit(1);
            return;
        }

        // filenames
        String geneFilename = cmd.getOptionValue("geneFile");
        String gffFilename = cmd.getOptionValue("gffFile");
        String vcfFilename = cmd.getOptionValue("vcfFile");

        // parameters
        int altMin = Integer.parseInt(cmd.getOptionValue("altMin"));
        double altFracMin = Double.parseDouble(cmd.getOptionValue("altFracMin"));

        // load the GFF file
        GFFLoader gffLoader = new GFFLoader(gffFilename);

        // load the VCF file
        VCFLoader vcfLoader = new VCFLoader(vcfFilename);
        vcfLoader.load();

        // output parameters
        System.out.println("org.ncgr.gwas.SNPReadCounter");
        System.out.println("geneFilename:"+geneFilename);
        System.out.println("vcfFilename:"+vcfFilename);
        System.out.println();

        // output header
        System.out.println("Gene\tSeq\tStart\tEnd\tStrand\tCount\tRF\tRR\tAF\tAR\tAltFrac");

        // run through the genes file, which contains the gene ID in the first column
        String line;
        BufferedReader reader = new BufferedReader(new FileReader(geneFilename));
        // skip the header
        reader.readLine();
        while ((line=reader.readLine())!=null) {
            String[] parts = line.split("\t");
            String gene = parts[0];
            // search for the gene, could be mapped to more than one place
            FeatureList genes = gffLoader.searchID(gene);
            for (FeatureI feature : genes) {
                String seq = feature.seqname();
                Location loc = feature.location();
                int start = loc.start();
                int end = loc.end();
                char strand = '+';
                // if minus strand, indicate with "-" but make start<end
                if (start<0) {
                    strand = '-';
                    int minusStart = -start;
                    int minusEnd = -end;
                    start = minusEnd;
                    end = minusStart;
                }
                int refForward = 0;
                int refReverse = 0;
                int altForward = 0;
                int altReverse = 0;
                int count = 0;
                // run through the variant calls within this gene
                List<VariantContext> vcList = vcfLoader.query(seq, start, end).toList();
                for (VariantContext vc : vcList) {
                    // read depths
                    List<Integer> dp4List = vc.getAttributeAsIntList("DP4", 0);
                    int rF = dp4List.get(0);
                    int rR = dp4List.get(1);
                    int aF = dp4List.get(2);
                    int aR = dp4List.get(3);
                    // include if enough ALT counts
                    if (aF+aR>=altMin) {
                        // increment counts
                        refForward += rF;
                        refReverse += rR;
                        altForward += aF;
                        altReverse += aR;
                        count++;
                    }
                }
                double altFrac = (double)(altForward+altReverse)/(double)(refForward+refReverse+altForward+altReverse);
                if (count>0 && altFrac>=altFracMin) {
                    System.out.println(gene+"\t"+seq+"\t"+start+"\t"+end+"\t"+strand+"\t"+count+"\t"+refForward+"\t"+refReverse+"\t"+altForward+"\t"+altReverse+"\t"+altFrac);
                }
            }
        }
        
    }
/**
     * Main class does all the work.
     */
Main class does all the work.
Main class does all the work .
public static void main(String[] args) throws Exception {
        Options options = new Options();
        Option geneFileOption = new Option("f", "geneFile", true, "Gene list (e.g. generated by SNPComparer)");
        geneFileOption.setRequired(true);
        options.addOption(geneFileOption);
        Option gffFileOption = new Option("g", "gffFile", true, "GFF file");
        gffFileOption.setRequired(true);
        options.addOption(gffFileOption);
        Option vcfFileOption = new Option("v", "vcfFile", true, "VCF file");
        vcfFileOption.setRequired(true);
        options.addOption(vcfFileOption);
        Option altMinOption = new Option("a", "altMin", true, "minimum ALT depth for SNP to be counted");
        altMinOption.setRequired(true);
        options.addOption(altMinOption);
        Option altFracMinOption = new Option("afm", "altFracMin", true, "minimum ALT fraction for a gene to be output");
        altFracMinOption.setRequired(true);
        options.addOption(altFracMinOption);
        CommandLineParser parser = new DefaultParser();
        HelpFormatter formatter = new HelpFormatter();
        CommandLine cmd;
        try {
            cmd = parser.parse(options, args);
        } catch (ParseException e) {
            System.out.println(e.getMessage());
            formatter.printHelp("SNPReadCounter", options);
            System.exit(1);
            return;
        }
        String geneFilename = cmd.getOptionValue("geneFile");
        String gffFilename = cmd.getOptionValue("gffFile");
        String vcfFilename = cmd.getOptionValue("vcfFile");
        int altMin = Integer.parseInt(cmd.getOptionValue("altMin"));
        double altFracMin = Double.parseDouble(cmd.getOptionValue("altFracMin"));
        GFFLoader gffLoader = new GFFLoader(gffFilename);
        VCFLoader vcfLoader = new VCFLoader(vcfFilename);
        vcfLoader.load();
        System.out.println("org.ncgr.gwas.SNPReadCounter");
        System.out.println("geneFilename:"+geneFilename);
        System.out.println("vcfFilename:"+vcfFilename);
        System.out.println();
        System.out.println("Gene\tSeq\tStart\tEnd\tStrand\tCount\tRF\tRR\tAF\tAR\tAltFrac");
        String line;
        BufferedReader reader = new BufferedReader(new FileReader(geneFilename));
        reader.readLine();
        while ((line=reader.readLine())!=null) {
            String[] parts = line.split("\t");
            String gene = parts[0];
            FeatureList genes = gffLoader.searchID(gene);
            for (FeatureI feature : genes) {
                String seq = feature.seqname();
                Location loc = feature.location();
                int start = loc.start();
                int end = loc.end();
                char strand = '+';
                if (start<0) {
                    strand = '-';
                    int minusStart = -start;
                    int minusEnd = -end;
                    start = minusEnd;
                    end = minusStart;
                }
                int refForward = 0;
                int refReverse = 0;
                int altForward = 0;
                int altReverse = 0;
                int count = 0;
                List<VariantContext> vcList = vcfLoader.query(seq, start, end).toList();
                for (VariantContext vc : vcList) {
                    List<Integer> dp4List = vc.getAttributeAsIntList("DP4", 0);
                    int rF = dp4List.get(0);
                    int rR = dp4List.get(1);
                    int aF = dp4List.get(2);
                    int aR = dp4List.get(3);
                    if (aF+aR>=altMin) {
                        refForward += rF;
                        refReverse += rR;
                        altForward += aF;
                        altReverse += aR;
                        count++;
                    }
                }
                double altFrac = (double)(altForward+altReverse)/(double)(refForward+refReverse+altForward+altReverse);
                if (count>0 && altFrac>=altFracMin) {
                    System.out.println(gene+"\t"+seq+"\t"+start+"\t"+end+"\t"+strand+"\t"+count+"\t"+refForward+"\t"+refReverse+"\t"+altForward+"\t"+altReverse+"\t"+altFrac);
                }
            }
        }
    }
public static void main ( String [ ] args ) throws Exception { Options options = new Options ( ) ; Option geneFileOption = new Option ( "f" , "geneFile" , true , "Gene list (e.g. generated by SNPComparer)" ) ; geneFileOption . setRequired ( true ) ; options . addOption ( geneFileOption ) ; Option gffFileOption = new Option ( "g" , "gffFile" , true , "GFF file" ) ; gffFileOption . setRequired ( true ) ; options . addOption ( gffFileOption ) ; Option vcfFileOption = new Option ( "v" , "vcfFile" , true , "VCF file" ) ; vcfFileOption . setRequired ( true ) ; options . addOption ( vcfFileOption ) ; Option altMinOption = new Option ( "a" , "altMin" , true , "minimum ALT depth for SNP to be counted" ) ; altMinOption . setRequired ( true ) ; options . addOption ( altMinOption ) ; Option altFracMinOption = new Option ( "afm" , "altFracMin" , true , "minimum ALT fraction for a gene to be output" ) ; altFracMinOption . setRequired ( true ) ; options . addOption ( altFracMinOption ) ; CommandLineParser parser = new DefaultParser ( ) ; HelpFormatter formatter = new HelpFormatter ( ) ; CommandLine cmd ; try { cmd = parser . parse ( options , args ) ; } catch ( ParseException e ) { System . out . println ( e . getMessage ( ) ) ; formatter . printHelp ( "SNPReadCounter" , options ) ; System . exit ( 1 ) ; return ; } String geneFilename = cmd . getOptionValue ( "geneFile" ) ; String gffFilename = cmd . getOptionValue ( "gffFile" ) ; String vcfFilename = cmd . getOptionValue ( "vcfFile" ) ; int altMin = Integer . parseInt ( cmd . getOptionValue ( "altMin" ) ) ; double altFracMin = Double . parseDouble ( cmd . getOptionValue ( "altFracMin" ) ) ; GFFLoader gffLoader = new GFFLoader ( gffFilename ) ; VCFLoader vcfLoader = new VCFLoader ( vcfFilename ) ; vcfLoader . load ( ) ; System . out . println ( "org.ncgr.gwas.SNPReadCounter" ) ; System . out . println ( "geneFilename:" + geneFilename ) ; System . out . println ( "vcfFilename:" + vcfFilename ) ; System . out . println ( ) ; System . out . println ( "Gene\tSeq\tStart\tEnd\tStrand\tCount\tRF\tRR\tAF\tAR\tAltFrac" ) ; String line ; BufferedReader reader = new BufferedReader ( new FileReader ( geneFilename ) ) ; reader . readLine ( ) ; while ( ( line = reader . readLine ( ) ) != null ) { String [ ] parts = line . split ( "\t" ) ; String gene = parts [ 0 ] ; FeatureList genes = gffLoader . searchID ( gene ) ; for ( FeatureI feature : genes ) { String seq = feature . seqname ( ) ; Location loc = feature . location ( ) ; int start = loc . start ( ) ; int end = loc . end ( ) ; char strand = '+' ; if ( start < 0 ) { strand = '-' ; int minusStart = - start ; int minusEnd = - end ; start = minusEnd ; end = minusStart ; } int refForward = 0 ; int refReverse = 0 ; int altForward = 0 ; int altReverse = 0 ; int count = 0 ; List < VariantContext > vcList = vcfLoader . query ( seq , start , end ) . toList ( ) ; for ( VariantContext vc : vcList ) { List < Integer > dp4List = vc . getAttributeAsIntList ( "DP4" , 0 ) ; int rF = dp4List . get ( 0 ) ; int rR = dp4List . get ( 1 ) ; int aF = dp4List . get ( 2 ) ; int aR = dp4List . get ( 3 ) ; if ( aF + aR >= altMin ) { refForward += rF ; refReverse += rR ; altForward += aF ; altReverse += aR ; count ++ ; } } double altFrac = ( double ) ( altForward + altReverse ) / ( double ) ( refForward + refReverse + altForward + altReverse ) ; if ( count > 0 && altFrac >= altFracMin ) { System . out . println ( gene + "\t" + seq + "\t" + start + "\t" + end + "\t" + strand + "\t" + count + "\t" + refForward + "\t" + refReverse + "\t" + altForward + "\t" + altReverse + "\t" + altFrac ) ; } } } }
Main class does all the work.
Main class does all the work .
// filenames // parameters // load the GFF file // load the VCF file // output parameters // output header // run through the genes file, which contains the gene ID in the first column // skip the header // search for the gene, could be mapped to more than one place // if minus strand, indicate with "-" but make start<end // run through the variant calls within this gene // read depths // include if enough ALT counts // increment counts