<commit_before>#!/usr/bin/perl

## <NAME>, 2017. MPI-CBG and MPI-PKS.
## The CESAR2.0 annotation engine

use strict;
use warnings;
use Scalar::Util::Numeric qw(isint);
use Getopt::Long qw(:config no_ignore_case no_auto_abbrev);
# this perl module must be located in the same dir as this script
use FindBin qw($Bin);
use lib "$Bin/";
use geneAnnotationFunctions;

# input options/parameters:
my $verbose = my $list = 0;
my $maxMem = 30;
GetOptions ("v|verbose" => \$verbose, "maxMemory=i" => \$maxMem, "list" =>\$list);


my $usage = "wrongs # args\n $0 - The core annotation engine that runs CESAR2.0 on a given genePred file and outputs 
coordinates of intact exons in a bed format.
Mandatory positional arguments:
1) input [transript identifier| a file containing transcripts where every line is one transcript]
2) the indexed alignment (.bb) file
3) the genepred file which contains information about all transcripts (output of formatGenePred.pl)
4) reference species --> for example: hg38
5) speciesList --> a comma separated string of species whose exons are being annotated. For example mm10,bosTau7 for mouse and cow
6) output directory --> a file with the same name as the query species is generated in this directory and coordinates of all intact exons are appended to this file
7) the twoBit directory that contains the genomes of the reference and all query species
8) the path to the CESAR binary. This path is expected to contain the extra/ subdirectory that contains donor and acceptor profiles and the substitution matrix
Options:
  -v|verbose
  -maxMemory int    allow CESAR 2.0 to allocate at most that much memory in Gb. Default is 30 (Gb).
  -list [when the input is a file containing one transcript per line] 
";

die $usage if (scalar(@ARGV) < 8);

my $input       = $ARGV[0];
my $mafIndex    = $ARGV[1];
my $genePred    = $ARGV[2];
my $reference   = $ARGV[3];
my $speciesList = $ARGV[4];
my $outDir      = $ARGV[5];
my $twoBitDir   = $ARGV[6];
my $cesarDir    = $ARGV[7];

## Check arguments
die "The mafIndex file '$mafIndex' does not exist\n" if (! -e $mafIndex);
die "The genePred file '$genePred' does not exist\n" if (! -e $genePred);
die "The genes list file '$input' does not exist, you specified that the input is a list\n" if (! -e $input && $list != 0);

my $clade = "human";

my @speciesArray = split(",",$speciesList);
my %speciesOfInterest = map{$_ => 1}@speciesArray;

## Step0: check if two bit files are present for all species of interest
foreach my $species(@speciesArray) {
	my $twoBitFile = "$twoBitDir/$species/$species.2bit";
	die "The two bit file ('$twoBitFile') for species '$species' is not present. Aborting\n" if (! -e $twoBitFile);
}
my $twoBitFileRef = "$twoBitDir/$reference/$reference.2bit";

my %genesOfInterest = ();
if ($list) {
	open(FI,$input) || die "Error opening genesList file '$input'\n";
	while (my $gene = <FI>) {
		chomp $gene;
		$genesOfInterest{$gene} = "T";
	}
	close FI;
} else {
	$genesOfInterest{$input} = "T";
}

## Step1: Get the information about the input transcripts/genes from the genePred file
my %geneInfoHash = ();
open(FI,$genePred) || die "Error opening genePred file\n";
while (my $line = <FI>) {
	$line =~s/\s+$//;
	my ($gene,$chr,$strand,$exons) = (split /\t/,$line)[0,1,2,7];
	if (exists $genesOfInterest{$gene}) {
		$geneInfoHash{$gene} = "$chr#$strand#$exons";
	}
}
close FI;

## check if there is nothing in the geneInfoHash
my @genesInt = keys(%geneInfoHash);
if (scalar(@genesInt == 0)) {
	print "############## ERROR ##############\n";
	print "No information found for the input gene '$input' from the genePred file '$genePred'\n";
	print "Could it be that the input is a file and you forgot to specify the '-list' parameter\nCheck and re-run\n";
	print "###################################\n";
	exit;
}

my @tempFiles = ();  ## create all temp files upfront and store them in @tempFiles array
my $mafFile      = `mktemp /dev/shm/exon.maf.XXXXXXXXX`; chomp $mafFile;
my $speciesBed   = `mktemp /dev/shm/species.bed.XXXXXX`; chomp $speciesBed;
my $speciesFasta = `mktemp /dev/shm/species.fa.XXXXXX`; chomp $speciesFasta;
my $cesarInput   = `mktemp /dev/shm/cesarIn.fa.XXXXXX`; chomp $cesarInput;
my $cesarOutput  = `mktemp /dev/shm/cesarOut.aln.XXXXXX`; chomp $cesarOutput;
push(@tempFiles,$mafFile,$speciesBed,$speciesFasta,$cesarInput,$cesarOutput);
	
foreach my $gene(keys(%geneInfoHash)) {
	my($chrRef,$strandRef,$allExons) = (split /#/,$geneInfoHash{$gene});

	print "The gene '$gene' comes on chromosome '$chrRef', strand '$strandRef'. Exons are '$allExons'\n" if ($verbose);
	my @exonsList = split(",",$allExons);
	@exonsList = reverse(@exonsList) if ($strandRef eq "-");

	## Step2: Get exon phases: i.e the split codon bases for each exons on both the 5' and the 3' end
	my $noe = scalar(@exonsList);
	my $ref2PhaseList = getExonPhaseList(\@exonsList);
	my @exonPhasesList = @$ref2PhaseList;

	## Step3: Run mafExtract on every exon and store the cdsStart, cdsStop, the chr and the strand for every query species present in the maf
	## At this stage also identify species-exon pairs which should  be excluded because:
	## 1) exon comes on more than one scaffold|strand, 
	## 2) nothing  but e-Lines 
	## 3) large insertions in the exon sequence of the query species

	my %ignoreExonsHash = ();
	my %cdsStartHash = my %cdsStopHash = my %strandHash = my %chrHash = ();
	my %exonsLengthHash = my %exonNumberHash = ();
	my %refSeqHash = my %referenceSpliceSites = ();

	my %hashInfo = ();
	## Initialize every entry in the hashInfo to NA at the start
	foreach my $species(@speciesArray) {
		my $k = 0;
		while ($k < $noe) {
			$hashInfo{$species}{$k}     = "NA#NA#NA#0";
			$k++;
		}
	}
	print "Processing gene '$gene'\n";

	my %exonHasStop = ();
	print "Fetching the alignment coordinates for each species and for each exon from '$gene' from the indexed maf\n###\n" if ($verbose); 
	my $exonCt = 1;
	foreach my $e(@exonsList) {
		print "Processing exon '$e'\n" if ($verbose);
		$exonNumberHash{$e} = $exonCt;
		my($refStart,$refStop) = (split /-/,$e)[0,1];
		my $lengthExon = abs($refStart - $refStop);
		$exonsLengthHash{$exonCt} = $lengthExon;
		
		## process reference sequence
		my ($phase5Prime,$phase3Prime) = (split /-/,$exonPhasesList[$exonCt-1])[0,1];
		my $refStartSS = $refStart - 2;  # the coordinates of the splice site
		my $refStopSS  = $refStop + 2;
		
		my $seqRefPrint = getSequenceFrom2Bit($twoBitFileRef,$chrRef,$refStartSS,$refStopSS);
		$seqRefPrint = revComp($seqRefPrint) if ($strandRef eq "-");
		
		## Get splice site nucleotides and find out if the acceptor is AC or the donor is AT
		print "The reference sequence is '$seqRefPrint'\n" if ($verbose);
		
		my $accSite = lc(substr($seqRefPrint,0,2));
		my $donorSite = lc(substr($seqRefPrint,length($seqRefPrint)-2,2));
		$referenceSpliceSites{$exonCt}{"acc"} = $accSite;
		$referenceSpliceSites{$exonCt}{"do"}  = $donorSite;
		
		## Get reference sequence that will be used as input to CESAR
		$seqRefPrint = substr($seqRefPrint,2,length($seqRefPrint)-4);  ## Get rid of the splice sites from the refSequence
		$seqRefPrint = fixPhase($seqRefPrint,$phase5Prime,$phase3Prime);
		
		## check if exon has stop codon:
		my $exonSeqInRef = $seqRefPrint;
		$exonSeqInRef = substr($exonSeqInRef,$phase5Prime,length($exonSeqInRef)-$phase5Prime-$phase3Prime);
		$exonSeqInRef =~s/\|//g;
		
		my $s= 0;
		while ($s < length($exonSeqInRef)) {
			my $codon = substr($exonSeqInRef,$s,1).substr($exonSeqInRef,$s+1,1).substr($exonSeqInRef,$s+2,1);

			if ($codon eq "TAA" || $codon eq "TAG" || $codon eq "TGA") {	 
				if ( ! ($s == length($exonSeqInRef) - 3) ) { ## This should not be the last codon at the end of the sequence though:
					$exonHasStop{$exonCt} = "T";
					print "'$exonSeqInRef' has stop codon at position '$s'\n" if ($verbose);
				}
			}
			$s = $s+3;
		}
		
		$refSeqHash{$exonCt} = $seqRefPrint;
		
		my $mafExtractCall = "set -o pipefail; mafExtract -region=$chrRef:$refStart-$refStop $mafIndex stdout|mafSpeciesSubset stdin speciesList=$speciesList,$reference species.lst=NULL $mafFile";
		system($mafExtractCall) == 0 || die "Error running '$mafExtractCall'\n";
		my($ref2cdsStartQuery,$ref2cdsStopQuery,$ref2strandQuery,$ref2chrQuery,$ref2moreThanOneStrand,$ref2moreThanOneChr,$ref2SpeciesList,$ref2cdsPrint,$ref2seqLength) = getCoordinatesFromMaf($mafFile,$reference);
		
		my @speciesPresentInMaf = @$ref2SpeciesList;
		foreach my $species(@speciesPresentInMaf) {
			
			my $cdsStart = $ref2cdsStartQuery->{$species};
			my $cdsStop  = $ref2cdsStopQuery->{$species};
			
			if ($ref2moreThanOneStrand->{$species} eq "T" || $ref2moreThanOneChr->{$species} eq "T" || ($cdsStart eq "NA" && $cdsStop eq "NA") || ($cdsStart > $cdsStop) || (abs($cdsStart-$cdsStop) - $lengthExon > 10000) ) {
				$ignoreExonsHash{$exonCt}{$species} = "T";<commit_msg>fix related to a bash specific feature in annotateGenesViaCESAR.pl<commit_after>#!/usr/bin/perl

## <NAME>, 2017. MPI-CBG and MPI-PKS.
## The CESAR2.0 annotation engine

use strict;
use warnings;
use Scalar::Util::Numeric qw(isint);
use Getopt::Long qw(:config no_ignore_case no_auto_abbrev);
# this perl module must be located in the same dir as this script
use FindBin qw($Bin);
use lib "$Bin/";
use geneAnnotationFunctions;

# input options/parameters:
my $verbose = my $list = 0;
my $maxMem = 30;
GetOptions ("v|verbose" => \$verbose, "maxMemory=i" => \$maxMem, "list" =>\$list);


my $usage = "wrongs # args\n $0 - The core annotation engine that runs CESAR2.0 on a given genePred file and outputs 
coordinates of intact exons in a bed format.
Mandatory positional arguments:
1) input [transript identifier| a file containing transcripts where every line is one transcript]
2) the indexed alignment (.bb) file
3) the genepred file which contains information about all transcripts (output of formatGenePred.pl)
4) reference species --> for example: hg38
5) speciesList --> a comma separated string of species whose exons are being annotated. For example mm10,bosTau7 for mouse and cow
6) output directory --> a file with the same name as the query species is generated in this directory and coordinates of all intact exons are appended to this file
7) the twoBit directory that contains the genomes of the reference and all query species
8) the path to the CESAR binary. This path is expected to contain the extra/ subdirectory that contains donor and acceptor profiles and the substitution matrix
Options:
  -v|verbose
  -maxMemory int    allow CESAR 2.0 to allocate at most that much memory in Gb. Default is 30 (Gb).
  -list [when the input is a file containing one transcript per line] 
";

die $usage if (scalar(@ARGV) < 8);

my $input       = $ARGV[0];
my $mafIndex    = $ARGV[1];
my $genePred    = $ARGV[2];
my $reference   = $ARGV[3];
my $speciesList = $ARGV[4];
my $outDir      = $ARGV[5];
my $twoBitDir   = $ARGV[6];
my $cesarDir    = $ARGV[7];

## Check arguments
die "The mafIndex file '$mafIndex' does not exist\n" if (! -e $mafIndex);
die "The genePred file '$genePred' does not exist\n" if (! -e $genePred);
die "The genes list file '$input' does not exist, you specified that the input is a list\n" if (! -e $input && $list != 0);

my $clade = "human";

my @speciesArray = split(",",$speciesList);
my %speciesOfInterest = map{$_ => 1}@speciesArray;

## Step0: check if two bit files are present for all species of interest
foreach my $species(@speciesArray) {
	my $twoBitFile = "$twoBitDir/$species/$species.2bit";
	die "The two bit file ('$twoBitFile') for species '$species' is not present. Aborting\n" if (! -e $twoBitFile);
}
my $twoBitFileRef = "$twoBitDir/$reference/$reference.2bit";

my %genesOfInterest = ();
if ($list) {
	open(FI,$input) || die "Error opening genesList file '$input'\n";
	while (my $gene = <FI>) {
		chomp $gene;
		$genesOfInterest{$gene} = "T";
	}
	close FI;
} else {
	$genesOfInterest{$input} = "T";
}

## Step1: Get the information about the input transcripts/genes from the genePred file
my %geneInfoHash = ();
open(FI,$genePred) || die "Error opening genePred file\n";
while (my $line = <FI>) {
	$line =~s/\s+$//;
	my ($gene,$chr,$strand,$exons) = (split /\t/,$line)[0,1,2,7];
	if (exists $genesOfInterest{$gene}) {
		$geneInfoHash{$gene} = "$chr#$strand#$exons";
	}
}
close FI;

## check if there is nothing in the geneInfoHash
my @genesInt = keys(%geneInfoHash);
if (scalar(@genesInt == 0)) {
	print "############## ERROR ##############\n";
	print "No information found for the input gene '$input' from the genePred file '$genePred'\n";
	print "Could it be that the input is a file and you forgot to specify the '-list' parameter\nCheck and re-run\n";
	print "###################################\n";
	exit;
}

my @tempFiles = ();  ## create all temp files upfront and store them in @tempFiles array
my $mafFile      = `mktemp /dev/shm/exon.maf.XXXXXXXXX`; chomp $mafFile;
my $speciesBed   = `mktemp /dev/shm/species.bed.XXXXXX`; chomp $speciesBed;
my $speciesFasta = `mktemp /dev/shm/species.fa.XXXXXX`; chomp $speciesFasta;
my $cesarInput   = `mktemp /dev/shm/cesarIn.fa.XXXXXX`; chomp $cesarInput;
my $cesarOutput  = `mktemp /dev/shm/cesarOut.aln.XXXXXX`; chomp $cesarOutput;
push(@tempFiles,$mafFile,$speciesBed,$speciesFasta,$cesarInput,$cesarOutput);
	
foreach my $gene(keys(%geneInfoHash)) {
	my($chrRef,$strandRef,$allExons) = (split /#/,$geneInfoHash{$gene});

	print "The gene '$gene' comes on chromosome '$chrRef', strand '$strandRef'. Exons are '$allExons'\n" if ($verbose);
	my @exonsList = split(",",$allExons);
	@exonsList = reverse(@exonsList) if ($strandRef eq "-");

	## Step2: Get exon phases: i.e the split codon bases for each exons on both the 5' and the 3' end
	my $noe = scalar(@exonsList);
	my $ref2PhaseList = getExonPhaseList(\@exonsList);
	my @exonPhasesList = @$ref2PhaseList;

	## Step3: Run mafExtract on every exon and store the cdsStart, cdsStop, the chr and the strand for every query species present in the maf
	## At this stage also identify species-exon pairs which should  be excluded because:
	## 1) exon comes on more than one scaffold|strand, 
	## 2) nothing  but e-Lines 
	## 3) large insertions in the exon sequence of the query species

	my %ignoreExonsHash = ();
	my %cdsStartHash = my %cdsStopHash = my %strandHash = my %chrHash = ();
	my %exonsLengthHash = my %exonNumberHash = ();
	my %refSeqHash = my %referenceSpliceSites = ();

	my %hashInfo = ();
	## Initialize every entry in the hashInfo to NA at the start
	foreach my $species(@speciesArray) {
		my $k = 0;
		while ($k < $noe) {
			$hashInfo{$species}{$k}     = "NA#NA#NA#0";
			$k++;
		}
	}
	print "Processing gene '$gene'\n";

	my %exonHasStop = ();
	print "Fetching the alignment coordinates for each species and for each exon from '$gene' from the indexed maf\n###\n" if ($verbose); 
	my $exonCt = 1;
	foreach my $e(@exonsList) {
		print "Processing exon '$e'\n" if ($verbose);
		$exonNumberHash{$e} = $exonCt;
		my($refStart,$refStop) = (split /-/,$e)[0,1];
		my $lengthExon = abs($refStart - $refStop);
		$exonsLengthHash{$exonCt} = $lengthExon;
		
		## process reference sequence
		my ($phase5Prime,$phase3Prime) = (split /-/,$exonPhasesList[$exonCt-1])[0,1];
		my $refStartSS = $refStart - 2;  # the coordinates of the splice site
		my $refStopSS  = $refStop + 2;
		
		my $seqRefPrint = getSequenceFrom2Bit($twoBitFileRef,$chrRef,$refStartSS,$refStopSS);
		$seqRefPrint = revComp($seqRefPrint) if ($strandRef eq "-");
		
		## Get splice site nucleotides and find out if the acceptor is AC or the donor is AT
		print "The reference sequence is '$seqRefPrint'\n" if ($verbose);
		
		my $accSite = lc(substr($seqRefPrint,0,2));
		my $donorSite = lc(substr($seqRefPrint,length($seqRefPrint)-2,2));
		$referenceSpliceSites{$exonCt}{"acc"} = $accSite;
		$referenceSpliceSites{$exonCt}{"do"}  = $donorSite;
		
		## Get reference sequence that will be used as input to CESAR
		$seqRefPrint = substr($seqRefPrint,2,length($seqRefPrint)-4);  ## Get rid of the splice sites from the refSequence
		$seqRefPrint = fixPhase($seqRefPrint,$phase5Prime,$phase3Prime);
		
		## check if exon has stop codon:
		my $exonSeqInRef = $seqRefPrint;
		$exonSeqInRef = substr($exonSeqInRef,$phase5Prime,length($exonSeqInRef)-$phase5Prime-$phase3Prime);
		$exonSeqInRef =~s/\|//g;
		
		my $s= 0;
		while ($s < length($exonSeqInRef)) {
			my $codon = substr($exonSeqInRef,$s,1).substr($exonSeqInRef,$s+1,1).substr($exonSeqInRef,$s+2,1);

			if ($codon eq "TAA" || $codon eq "TAG" || $codon eq "TGA") {	 
				if ( ! ($s == length($exonSeqInRef) - 3) ) { ## This should not be the last codon at the end of the sequence though:
					$exonHasStop{$exonCt} = "T";
					print "'$exonSeqInRef' has stop codon at position '$s'\n" if ($verbose);
				}
			}
			$s = $s+3;
		}
		
		$refSeqHash{$exonCt} = $seqRefPrint;
		
		my $mafExtractCall = "/bin/bash -c 'set -o pipefail; mafExtract -region=$chrRef:$refStart-$refStop $mafIndex stdout|mafSpeciesSubset stdin NULL $mafFile -speciesList=$speciesList,$reference'";
		print "Running $mafExtractCall\n" if ($verbose);
		system($mafExtractCall) == 0 || die "Error running '$mafExtractCall'\n";
		my($ref2cdsStartQuery,$ref2cdsStopQuery,$ref2strandQuery,$ref2chrQuery,$ref2moreThanOneStrand,$ref2moreThanOneChr,$ref2SpeciesList,$ref2cdsPrint,$ref2seqLength) = getCoordinatesFromMaf($mafFile,$reference);
		
		my @speciesPresentInMaf = @$ref2SpeciesList;
		foreach my $species(@speciesPresentInMaf) {
			
			my $cdsStart = $ref2cdsStartQuery->{$species};
			my $cdsStop  = $ref2cdsStopQuery->{$species};
			
			if ($ref2moreThanOneStrand->{$species} eq "T" || $ref2moreThanOneChr->{$species} eq "T" || ($cdsStart eq "NA" && $cdsStop eq "NA") || ($cdsStart > $cdsStop) || (abs($cdsStart-$cdsStop) - $lengthExon > 10000) ) {
				$ignoreExonsHash{$exonCt}{$species} = "T";
{'fraction_non_alphanumeric': 0.15093384338433843, 'fraction_numerical': 0.012376237623762377, 'mean_word_length': 4.645284217211813, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, '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': '6414380', 'n_tokens_mistral': 6683, 'n_tokens_neox': 6047, 'n_words': 1912}