<commit_before>#!/usr/bin/env perl
#-------------------------------------------------------------------------------
# Purpose:	Calculate pi, theta, and Tajima's D and print to STDOUT
# Author:	<NAME>
# Created:	9/21/15
#-------------------------------------------------------------------------------

use warnings;
use strict;
use Getopt::Long;
use Data::Dumper;
use File::Path qw(make_path);
use File::Spec::Functions qw(catdir catfile);
use File::Touch;
use CoGe::Algos::PopGen::Diversity;
use CoGe::Algos::PopGen::FileFormats;
use Tabix;

my ($VCF_FILE, $GFF_FILE, $FASTA_FILE, $CHR, $GENE_NAME, $FEAT_TYPE, 
    $OUTPUT_PATH, $DEBUG, $DEBUGFH);
GetOptions(
    "vcf=s"    => \$VCF_FILE,     # required VCF input file
    "gff=s"    => \$GFF_FILE,     # required GFF input file
    "fasta=s"  => \$FASTA_FILE,   # required reference FASTA input file
    "output=s" => \$OUTPUT_PATH,  # optional output path for results files
    "gene=s"   => \$GENE_NAME,    # optional gene name for testing
    "type=s"   => \$FEAT_TYPE,    # optional GFF type for testing
    "chr=s"    => \$CHR,          # optional chromosome name for testing
    "debug"    => \$DEBUG         # optional debug output
);

my $NUM_INDIVIDUALS = 12 * 2;

# Create output path if specified
if (defined $OUTPUT_PATH) {
    make_path($OUTPUT_PATH) unless (-r $OUTPUT_PATH);
    unless (-r $OUTPUT_PATH) {
        print STDERR "Error: couldn't create output path '$OUTPUT_PATH'\n";
        exit(-1);
    }
}
else {
    $OUTPUT_PATH = './';
}

# Create debug file
if ($DEBUG) {
    my $logfile = catfile($OUTPUT_PATH, 'sumstats.log');
    open($DEBUGFH, ">$logfile") or die "Cannot open log output file '$logfile'\n";
}

# Load GFF
print STDERR "Loading GFF\n";
my $pAnnot = loadGFF(file => $GFF_FILE);

# Load FASTA
print STDERR "Loading FASTA\n";
my $pSeq = loadFASTA($FASTA_FILE);

# Open Tabix indexed VCF
print STDERR "Loading VCF\n";
my $tabix = Tabix->new(-data => $VCF_FILE);

# Create results file
my $resultsfile = catfile($OUTPUT_PATH, 'sumstats.tsv');
open(my $fh, ">$resultsfile") or die "Cannot open results output file '$resultsfile'\n";

# Iterate through GFF entities and calculate summary stats
print STDERR "Calculating pi, theta, and Tajima's D\n";
my %allele_cache;
foreach my $type (sort keys %$pAnnot) {
    next if (defined $FEAT_TYPE and $type ne $FEAT_TYPE);
    
    foreach my $chr (sort keys %{$pAnnot->{$type}}) {
        next if (defined $CHR and $chr ne $CHR);
    
        # Print header line
        print $fh join("\t", "#$type $chr", 'GENE NAME', 'START', 'END', 
            'TOTAL SITES', 'TOTAL SEG. SITES', 'TOTAL PI', 'TOTAL THETA', "TOTAL TAJIMA'S D");
        print $fh "\t", join("\t", 
            '0-FOLD SITES', '0-FOLD SEG. SITES', '0-FOLD PI', '0-FOLD THETA', "0-FOLD TAJIMA'S D", 
            '4-FOLD SITES', '4-FOLD SEG. SITES', '4-FOLD PI', '4-FOLD THETA', "4-FOLD TAJIMA'S D")
            if ($type eq 'cds');
        print $fh "\n";
            
        foreach my $id (sort keys %{$pAnnot->{$type}{$chr}}) {
            next if (defined $GENE_NAME and $id ne $GENE_NAME);
            
            my %stats;
            #my ($sites, $segregating, $compared, $mismatch) = (0, 0, 0, 0);
            my ($codingOffset, $codingSeq, $featStart, $featEnd, $featChr);
            $codingOffset = 0;<commit_msg>Move chr into separate header<commit_after>#!/usr/bin/env perl
#-------------------------------------------------------------------------------
# Purpose:	Calculate pi, theta, and Tajima's D and print to STDOUT
# Author:	<NAME>
# Created:	9/21/15
#-------------------------------------------------------------------------------

use warnings;
use strict;
use Getopt::Long;
use Data::Dumper;
use File::Path qw(make_path);
use File::Spec::Functions qw(catdir catfile);
use File::Touch;
use CoGe::Algos::PopGen::Diversity;
use CoGe::Algos::PopGen::FileFormats;
use Tabix;

my ($VCF_FILE, $GFF_FILE, $FASTA_FILE, $CHR, $GENE_NAME, $FEAT_TYPE, 
    $OUTPUT_PATH, $DEBUG, $DEBUGFH);
    
GetOptions(
    "vcf=s"    => \$VCF_FILE,     # required VCF input file
    "gff=s"    => \$GFF_FILE,     # required GFF input file
    "fasta=s"  => \$FASTA_FILE,   # required reference FASTA input file
    "output=s" => \$OUTPUT_PATH,  # optional output path for results files
    "gene=s"   => \$GENE_NAME,    # optional gene name for testing
    "type=s"   => \$FEAT_TYPE,    # optional GFF type for testing
    "chr=s"    => \$CHR,          # optional chromosome name for testing
    "debug"    => \$DEBUG         # optional debug output
);

my $NUM_INDIVIDUALS = 12 * 2;

# Create output path if specified
if (defined $OUTPUT_PATH) {
    make_path($OUTPUT_PATH) unless (-r $OUTPUT_PATH);
    unless (-r $OUTPUT_PATH) {
        print STDERR "Error: couldn't create output path '$OUTPUT_PATH'\n";
        exit(-1);
    }
}
else {
    $OUTPUT_PATH = './';
}

# Create debug file
if ($DEBUG) {
    my $logfile = catfile($OUTPUT_PATH, 'sumstats.log');
    open($DEBUGFH, ">$logfile") or die "Cannot open log output file '$logfile'\n";
}

# Load GFF
print STDERR "Loading GFF\n";
my $pAnnot = loadGFF(file => $GFF_FILE);

# Load FASTA
print STDERR "Loading FASTA\n";
my $pSeq = loadFASTA($FASTA_FILE);

# Open Tabix indexed VCF
print STDERR "Loading VCF\n";
my $tabix = Tabix->new(-data => $VCF_FILE);

# Create results file
my $resultsfile = catfile($OUTPUT_PATH, 'sumstats.tsv');
open(my $fh, ">$resultsfile") or die "Cannot open results output file '$resultsfile'\n";

# Iterate through GFF entities and calculate summary stats
print STDERR "Calculating pi, theta, and Tajima's D\n";
my %allele_cache;
foreach my $type (sort keys %$pAnnot) {
    next if (defined $FEAT_TYPE and $type ne $FEAT_TYPE);
    
    # Print header line
    print $fh join("\t", "#$type", 'GENE NAME', 'START', 'END', 
        'TOTAL SITES', 'TOTAL SEG. SITES', 'TOTAL PI', 'TOTAL THETA', "TOTAL TAJIMA'S D");
    print $fh "\t", join("\t", 
        '0-FOLD SITES', '0-FOLD SEG. SITES', '0-FOLD PI', '0-FOLD THETA', "0-FOLD TAJIMA'S D", 
        '4-FOLD SITES', '4-FOLD SEG. SITES', '4-FOLD PI', '4-FOLD THETA', "4-FOLD TAJIMA'S D")
        if ($type eq 'cds');
    print $fh "\n";
    
    foreach my $chr (sort keys %{$pAnnot->{$type}}) {
        next if (defined $CHR and $chr ne $CHR);
    
        print $fh "#$chr\n";
            
        foreach my $id (sort keys %{$pAnnot->{$type}{$chr}}) {
            next if (defined $GENE_NAME and $id ne $GENE_NAME);
            
            my %stats;
            #my ($sites, $segregating, $compared, $mismatch) = (0, 0, 0, 0);
            my ($codingOffset, $codingSeq, $featStart, $featEnd, $featChr);
            $codingOffset = 0;
{'fraction_non_alphanumeric': 0.19796040791841632, 'fraction_numerical': 0.007198560287942412, 'mean_word_length': 2.5586979722518675, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1522755', 'n_tokens_mistral': 2445, 'n_tokens_neox': 2320, 'n_words': 689}