<commit_before>#!/usr/local/bin/perl

##############################################################################
# Used to obtain full data from logistic regression model
# loops through reference genome and outputs 1 line per base, as long as
# valid covariate data exists
##############################################################################

use strict;
use warnings;
use POSIX;
use Getopt::Long;
use Pod::Usage;
use File::Basename;
use File::Path qw(make_path);
use List::Util qw(first max maxstr min minstr reduce shuffle sum);
use List::MoreUtils 'pairwise';
use Cwd;
use Benchmark;
use Tie::File;

# Set options and inputs
my $wdir=getcwd;
my $parentdir="/net/bipolar/jedidiah/mutation";

my $baseopt;
my $chr;
my $categ;
my $bw = 100;
my $adj=1;

GetOptions ('b=s'=> \$baseopt,
'chr=s'=> \$chr,
'categ=s' => \$categ,
'bw=i' => \$bw,
'adj=i' => \$adj) or pod2usage(1);

my $f_covs = "$parentdir/output/logmod_data/${bw}kb_mut_cov2.txt";

make_path("$parentdir/output/logmod_data/chr${chr}/");

my $b1;
my $b2;
if($baseopt eq "AT"){
	$b1="A";
	$b2="T";
} else {
	$b1="C";
	$b2="G";
}

my $mask_flag=0;
my $binwidth=$bw*1000;

my $nextchr;
if ($chr<22) {
	$nextchr=$chr+1;
} elsif ($chr==22) {
	$nextchr="X";
}

my $subseq=2;
if ($adj!=0) {
	$subseq = $adj*2+1;
}

# initialize output file
my $outfile = "$parentdir/output/logmod_data/chr${chr}_${categ}_sites.txt";
my $OUT;
open($OUT, '>', $outfile) or die "can't write to $outfile: $!\n";

# initialize covariate data
open my $covs, '<', $f_covs or die "can't open $f_covs: $!";

# initialize singleton file
my $f_positions = "$parentdir/output/logmod_data/chr${chr}_${categ}_pos_examples.txt"; #main line for full processing
open my $positions, '<', $f_positions or die "can't open $f_positions: $!";

# Index motif file names
my $f_mlist = "$parentdir/output/7bp_1000k_rates.txt";
open my $mlist, '<', $f_mlist or die "can't open $f_mlist: $!";

our %fhash=();
my @fn;
while(<$mlist>){
  chomp;
  my @line=split(/\t/, $_);
  my $seq=$line[1];
  my $cat=$line[2];

  if($cat eq $categ){
    #my $key=join("\t", @line[0 .. 1]);
    #my $pcs=join("\t", @line[2 .. $#line]);

    my $filename="$parentdir/output/logmod_data/chr${chr}/chr${chr}_${categ}_$seq.txt";
    push(@fn, $filename);
    $fhash{$seq}=$filename;
    # print "$hash{$_}\n";
  }
}

my %handles = get_write_handles(@fn);

# initialize phastCons data
# my $f_cons = "$parentdir/reference_data/chr$chr.phastCons46way.primates.wigFix";
# open my $cons, '<', $f_cons or die "can't open $f_cons: $!";

# Get reference sequence
my $seq=&getRef();
my $altseq=$seq;
$altseq =~ tr/ACGT/TGCA/;

my $seqlength=length($seq);
# print "seqlength of chr$chr: $max\n"; #<-used to validate that getRef() returns correct seq length

my $printheader=0;
if($printheader==1){
	print $OUT "CHR \t POS \t BIN \t Sequence \t mut \n"; #<-add header to output, if needed
}

# Create hash keyed by Chr/Bin pairs, with row of PCs as value
print "Indexing chr${chr} covariate data...\n";
our %hash=();
while (<$covs>){
	chomp;
	my @line=split(/\t/, $_);
	my $key=join("\t", @line[0 .. 1]);
	my $pcs=join("\t", @line[2 .. $#line]);

	$hash{$key}=$pcs;
}

# my $key=join("\t", 20, 100);
# print "$hash{$key}\n";

# Create hash keyed by singleton positions, with input line as value
print "Indexing chr${chr}: ${categ} singleton positions...\n";
my @POS;
my %poshash=();
while (<$positions>) {
	chomp;
	my @line=split(/\t/, $_);
	my $key=$line[2];
	push (@POS, $key);

	$poshash{$key}=$_;
}

# Write data files
print "Writing chr${chr}: ${categ} data file...\n";
for my $strpos (0 .. $seqlength){
	my $base = substr($seq, $strpos, 1);
	my $pos = $strpos+1;

	my $bin = ceil($pos/$binwidth);
	my $key2=join("\t", $chr, $bin);

	if(defined $hash{$key2}){
		if(($base =~ /$b1|$b2/) & (!exists $poshash{$pos})){
			# push (@POS, $pos); # add position to exclusion list
			my $localseq = substr($seq, $pos-$adj-1, $subseq);
			my $altlocalseq = reverse substr($altseq, $pos-$adj-1, $subseq);

			# Coerce local sequence info to format used in R
			my $sequence;
			if(substr($localseq,$adj,1) lt substr($altlocalseq,$adj,1)){
				$sequence = $localseq . '(' . $altlocalseq . ')';
			} else {
				$sequence = $altlocalseq . '(' . $localseq . ')';
			}

			# write line if site has non-N context
			if ($sequence !~ /N/) {
				my $covs=&updateCovs($chr, $bin, $pos);

				my $file=$fhash{$sequence};
				# my $mref={$handles{$file}};

				print $OUT "$chr\t$bin\t$pos\t$sequence\t 0 \t$covs\n";
				print {$handles{$file}} "$chr\t$bin\t$pos\t$sequence\t 0 \t$covs\n";
			}
		}elsif(exists $poshash{$pos}){
			my $covs=&updateCovs($chr, $bin, $pos);
			my @line=split(/\t/, $poshash{$pos});
			my $sequence=$line[3];

			my $file=$fhash{$sequence};
			# my $mref={$handles{$file}};

			print $OUT "$poshash{$pos}\t$covs\n";
			print {$handles{$file}} "$poshash{$pos}\t$covs\n";<commit_msg>Fix noncanonical motif error on chr3<commit_after>#!/usr/local/bin/perl

##############################################################################
# Used to obtain full data from logistic regression model
# loops through reference genome and outputs 1 line per base, as long as
# valid covariate data exists
##############################################################################

use strict;
use warnings;
use POSIX;
use Getopt::Long;
use Pod::Usage;
use File::Basename;
use File::Path qw(make_path);
use List::Util qw(first max maxstr min minstr reduce shuffle sum);
use List::MoreUtils 'pairwise';
use Cwd;
use Benchmark;
use Tie::File;

# Set options and inputs
my $wdir=getcwd;
my $parentdir="/net/bipolar/jedidiah/mutation";

my $baseopt;
my $chr;
my $categ;
my $bw = 100;
my $adj=1;

GetOptions ('b=s'=> \$baseopt,
'chr=s'=> \$chr,
'categ=s' => \$categ,
'bw=i' => \$bw,
'adj=i' => \$adj) or pod2usage(1);

my $f_covs = "$parentdir/output/logmod_data/${bw}kb_mut_cov2.txt";

make_path("$parentdir/output/logmod_data/chr${chr}/");

my $b1;
my $b2;
if($baseopt eq "AT"){
	$b1="A";
	$b2="T";
} else {
	$b1="C";
	$b2="G";
}

my $mask_flag=0;
my $binwidth=$bw*1000;

my $nextchr;
if ($chr<22) {
	$nextchr=$chr+1;
} elsif ($chr==22) {
	$nextchr="X";
}

my $subseq=2;
if ($adj!=0) {
	$subseq = $adj*2+1;
}

# initialize output file
my $outfile = "$parentdir/output/logmod_data/chr${chr}_${categ}_sites.txt";
my $OUT;
open($OUT, '>', $outfile) or die "can't write to $outfile: $!\n";

# initialize covariate data
open my $covs, '<', $f_covs or die "can't open $f_covs: $!";

# initialize singleton file
my $f_positions = "$parentdir/output/logmod_data/chr${chr}_${categ}_pos_examples.txt"; #main line for full processing
open my $positions, '<', $f_positions or die "can't open $f_positions: $!";

# Index motif file names
my $f_mlist = "$parentdir/output/7bp_1000k_rates.txt";
open my $mlist, '<', $f_mlist or die "can't open $f_mlist: $!";

our %fhash=();
my @fn;
while(<$mlist>){
  chomp;
  my @line=split(/\t/, $_);
  my $seq=$line[1];
  my $cat=$line[2];

  if($cat eq $categ){
    #my $key=join("\t", @line[0 .. 1]);
    #my $pcs=join("\t", @line[2 .. $#line]);

    my $filename="$parentdir/output/logmod_data/chr${chr}/chr${chr}_${categ}_$seq.txt";
    push(@fn, $filename);
    $fhash{$seq}=$filename;
    # print "$hash{$_}\n";
  }
}

my %handles = get_write_handles(@fn);

# initialize phastCons data
# my $f_cons = "$parentdir/reference_data/chr$chr.phastCons46way.primates.wigFix";
# open my $cons, '<', $f_cons or die "can't open $f_cons: $!";

# Get reference sequence
my $seq=&getRef();
my $altseq=$seq;
$altseq =~ tr/ACGT/TGCA/;

my $seqlength=length($seq);
# print "seqlength of chr$chr: $max\n"; #<-used to validate that getRef() returns correct seq length

my $printheader=0;
if($printheader==1){
	print $OUT "CHR \t POS \t BIN \t Sequence \t mut \n"; #<-add header to output, if needed
}

# Create hash keyed by Chr/Bin pairs, with row of PCs as value
print "Indexing chr${chr} covariate data...\n";
our %hash=();
while (<$covs>){
	chomp;
	my @line=split(/\t/, $_);
	my $key=join("\t", @line[0 .. 1]);
	my $pcs=join("\t", @line[2 .. $#line]);

	$hash{$key}=$pcs;
}

# my $key=join("\t", 20, 100);
# print "$hash{$key}\n";

# Create hash keyed by singleton positions, with input line as value
print "Indexing chr${chr}: ${categ} singleton positions...\n";
my @POS;
my %poshash=();
while (<$positions>) {
	chomp;
	my @line=split(/\t/, $_);
	my $key=$line[2];
	push (@POS, $key);

	$poshash{$key}=$_;
}

# Write data files
print "Writing chr${chr}: ${categ} data file...\n";
for my $strpos (0 .. $seqlength){
	my $base = substr($seq, $strpos, 1);
	my $pos = $strpos+1;

	my $bin = ceil($pos/$binwidth);
	my $key2=join("\t", $chr, $bin);

	if(defined $hash{$key2}){
		if(($base =~ /$b1|$b2/) & (!exists $poshash{$pos})){
			# push (@POS, $pos); # add position to exclusion list
			my $localseq = substr($seq, $pos-$adj-1, $subseq);
			my $altlocalseq = reverse substr($altseq, $pos-$adj-1, $subseq);

			# Coerce local sequence info to format used in R
			my $sequence;
			if(substr($localseq,$adj,1) lt substr($altlocalseq,$adj,1)){
				$sequence = $localseq . '(' . $altlocalseq . ')';
			} else {
				$sequence = $altlocalseq . '(' . $localseq . ')';
			}

			# write line if site has non-N context
			if ($sequence !~ /MNSW/) {
				my $covs=&updateCovs($chr, $bin, $pos);

				my $file=$fhash{$sequence};
				# my $mref={$handles{$file}};

				print $OUT "$chr\t$bin\t$pos\t$sequence\t 0 \t$covs\n";
				print {$handles{$file}} "$chr\t$bin\t$pos\t$sequence\t 0 \t$covs\n";
			}
		}elsif(exists $poshash{$pos}){
			my $covs=&updateCovs($chr, $bin, $pos);
			my @line=split(/\t/, $poshash{$pos});
			my $sequence=$line[3];

			my $file=$fhash{$sequence};
			# my $mref={$handles{$file}};

			print $OUT "$poshash{$pos}\t$covs\n";
			print {$handles{$file}} "$poshash{$pos}\t$covs\n";
{'fraction_non_alphanumeric': 0.24066135946111453, 'fraction_numerical': 0.013370075525617474, 'mean_word_length': 4.505056179775281, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 21, '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': '1941585', 'n_tokens_mistral': 4148, 'n_tokens_neox': 3839, 'n_words': 966}