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

use strict;

=head1 NAME

cdna2genome

=head1 SYNOPSIS

    cdna2genome fasta-file-genome list-of-accession-numbers

    cdna2genome -timdb fasta-file-genome-sequence (*)

    -output [gff/embl]
    -nosim4 - do not show sim4 exons/exonsets
    -nocds  - do not show cds 
    -timdb  - switch to turn on ensembl pipeline processing

=head1 DESCRIPTION

This module compares a genomic sequence as fasta file
to a series of cDNAs designated by EMBL accession numbers.
These EMBL accession numbers are extracted and then
aligned using Sim4. Then the EMBL based CDS line is used
to provide coding sequence prediction.

The results are then merged to remove completely identical 
predictions and given out in GFF format

(*) full file specification of a sequence file containing multiple
contigs for an ensembl pipeline analysis directory.  Each contig is
processed in turn and the accessions of vertrna hits to these clones
are pulled out by scanning the hits to each predicted peptide for that
contig.  Output written to file in analysis directory
contig.sim4_vert.gff

=cut

use Bio::Tools::Sim4::Results;
use Bio::AnnSeq;
use Bio::SeqIO;
use Bio::AnnSeqIO;

use Bio::EnsEMBL::Analysis::Log;

use Getopt::Long;


my $log = new Bio::EnsEMBL::Analysis::Log;

my $output = 'gff';
my $nosim4 = 0;
my $nocds  = 0;
my $nologerr = 0;
my $timdb;

&GetOptions( 
	     'output:s' => \$output,
	     'nosim4'   => \$nosim4, 
	     'nocds'    => \$nocds,
             'nolog'    => \$nologerr,
	     'timdb'    => \$timdb,
	     );


my $format = 'Fasta';
my $file  = shift;

my %top_sf;

$file || do { system("perldoc $0"); exit(0) };

# get directory and clonename if timdb
my($dir,$clone,@files);
if($timdb){
    if($file=~/(.*)\/([^\/]+)\.seq/){
	$dir=$1;
	$clone=$2;
	local *TMPDIR;
	opendir(TMPDIR,$dir);
	@files=readdir(TMPDIR);
	closedir(TMPDIR);
    }else{
	die "cannot parse $file as timdb seq file";
    }
}

my $seqin = Bio::SeqIO->new( -format => $format , -file => $file );

# file could have more than sequences, so process them in turn
my $seq;
my $nseq;
while($seq = $seqin->next_seq()){
    $nseq++;

    my @embl_acc;
    my $OUTFILE;
    if($timdb){
	# get list of targets from an unfinished_ana directory if specified

	# get contig name and output file
	my $contig=$seq->id;
	#print "$contig\n";

	# loop over peptide hit files for this contig
	my %acc;
	foreach my $file (@files){
	    if($file=~/$contig\.(\d+)\.tblastn_vert\.msptmp$/){
		#print "  $file\n";
		my $pepid=$1;
		local *TMPFILE;
		open(TMPFILE,"$dir/$file") || die "cannot open $file";
		while(<TMPFILE>){
		    s/^\s+//;
		    my @fields=split(/\s+/,$_);
		    $acc{$fields[8]}.=$pepid.",";
		}
		close(TMPFILE);
	    }
	}
	@embl_acc=(keys %acc);
	
	# next unless msptmp hits
	next unless scalar(@embl_acc);

	my $outfile="$dir/$contig.sim4_vert.$output";
	open(SIM4OUT,">$outfile") || die "cannot open $outfile";
	$OUTFILE=\*SIM4OUT;

	# debug
	#foreach my $acc (@embl_acc){
	#    print $OUTFILE "$contig $acc $acc{$acc}\n";
	#}
 
    }else{
	@embl_acc=@ARGV;
	$OUTFILE=\*STDOUT;
    }

    &_do_cdna2genome($OUTFILE,\@embl_acc);

    if($timdb){
	close($OUTFILE);
    }

}

# if never found a sequence, report
unless($nseq){
    die "no sequence read from $file";
}

# this is the end of the script now
exit 0;

sub _do_cdna2genome{
    my($OUTFILE,$rembl_acc)=@_;

    my $aseq = Bio::AnnSeq->new();
    $aseq->seq($seq);

    my $OUTFILE;

  ACC: 
    foreach my $embl_acc ( @$rembl_acc ) {

	# get out the cdna.
	my $aseqio = Bio::AnnSeqIO->new( -format => 'EMBL', 
					 -file => "getz -e '[embl-acc:$embl_acc]' |");    
	my $cdna;
	eval {
	    $cdna = $aseqio->next_annseq();
	};

	# skip if we can't load it or no cdna.
	if( $@ || !$cdna) {
	    $log->store("Unable to read $embl_acc. Skipping. (Apologies)\n$@\n");
	    next ACC;
	}

	# get out CDS coordinates.
	
	my $cds_start;
	my $cds_end;
	
	foreach my $sf ( $cdna->top_SeqFeatures ) {
	    if( $sf->primary_tag() eq 'CDS' ) {
		$cds_start = $sf->start;
		$cds_end   = $sf->end;
	    }
	}
	
	if( !defined $cds_start ) {
	    print STDERR "EMBL $embl_acc has not single CDS line. Skipping\n";
	}<commit_msg>cdna edits<commit_after>#!/usr/local/bin/perl

use strict;

#  perl ../scripts/cdna2genome.pl -timdb ~/th/unfinished_ana/sanger/data/dJ718J7/dJ718J7.seq

=head1 NAME

cdna2genome

=head1 SYNOPSIS

    cdna2genome fasta-file-genome list-of-accession-numbers

    cdna2genome -timdb fasta-file-genome-sequence (*)

    -output [gff/embl]
    -nosim4 - do not show sim4 exons/exonsets
    -nocds  - do not show cds 
    -timdb  - switch to turn on ensembl pipeline processing

=head1 DESCRIPTION

This module compares a genomic sequence as fasta file
to a series of cDNAs designated by EMBL accession numbers.
These EMBL accession numbers are extracted and then
aligned using Sim4. Then the EMBL based CDS line is used
to provide coding sequence prediction.

The results are then merged to remove completely identical 
predictions and given out in GFF format

(*) full file specification of a sequence file containing multiple
contigs for an ensembl pipeline analysis directory.  Each contig is
processed in turn and the accessions of vertrna hits to these clones
are pulled out by scanning the hits to each predicted peptide for that
contig.  Output written to file in analysis directory
contig.sim4_vert.gff

=cut

use Bio::Tools::Sim4::Results;
use Bio::AnnSeq;
use Bio::SeqIO;
use Bio::AnnSeqIO;

use Bio::EnsEMBL::Analysis::Log;

use Getopt::Long;


my $log = new Bio::EnsEMBL::Analysis::Log;

my $output = 'gff';
my $nosim4 = 0;
my $nocds  = 0;
my $nologerr = 0;
my $timdb;

&GetOptions( 
	     'output:s' => \$output,
	     'nosim4'   => \$nosim4, 
	     'nocds'    => \$nocds,
             'nolog'    => \$nologerr,
	     'timdb'    => \$timdb,
	     );


my $format = 'Fasta';
my $file  = shift;

my %top_sf;

$file || do { system("perldoc $0"); exit(0) };

# get directory and clonename if timdb
my($dir,$clone,@files);
if($timdb){
    if($file=~/(.*)\/([^\/]+)\.seq/){
	$dir=$1;
	$clone=$2;
	local *TMPDIR;
	opendir(TMPDIR,$dir);
	@files=readdir(TMPDIR);
	closedir(TMPDIR);
    }else{
	die "cannot parse $file as timdb seq file";
    }
}

my $seqin = Bio::SeqIO->new( -format => $format , -file => $file );

# file could have more than sequences, so process them in turn
my $seq;
my $nseq;
while($seq = $seqin->next_seq()){
    $nseq++;

    my @embl_acc;
    my $OUTFILE;
    if($timdb){
	# get list of targets from an unfinished_ana directory if specified

	# get contig name and output file
	my $contig=$seq->id;
	#print "$contig\n";

	# loop over peptide hit files for this contig
	my %acc;
	foreach my $file (@files){
	    if($file=~/$contig\.(\d+)\.tblastn_vert\.msptmp$/){
		#print "  $file\n";
		my $pepid=$1;
		local *TMPFILE;
		open(TMPFILE,"$dir/$file") || die "cannot open $file";
		while(<TMPFILE>){
		    s/^\s+//;
		    my @fields=split(/\s+/,$_);
		    $acc{$fields[8]}.=$pepid.",";
		}
		close(TMPFILE);
	    }
	}
	@embl_acc=(keys %acc);
	
	# next unless msptmp hits
	next unless scalar(@embl_acc);

	my $outfile="$dir/$contig.sim4_vert.$output";
	open(SIM4OUT,">$outfile") || die "cannot open $outfile";
	$OUTFILE=\*SIM4OUT;

	# debug
	#foreach my $acc (@embl_acc){
	#    print $OUTFILE "$contig $acc $acc{$acc}\n";
	#}
 
    }else{
	@embl_acc=@ARGV;
	$OUTFILE=\*STDOUT;
    }

    &_do_cdna2genome($OUTFILE,\@embl_acc);

    if($timdb){
	close($OUTFILE);
    }

}

# if never found a sequence, report
unless($nseq){
    die "no sequence read from $file";
}

# this is the end of the script now
exit 0;

sub _do_cdna2genome{
    my($OUTFILE,$rembl_acc)=@_;

    my $aseq = Bio::AnnSeq->new();
    $aseq->seq($seq);

    #my $OUTFILE;

  ACC: 
    foreach my $embl_acc ( @$rembl_acc ) {

	# get out the cdna.
	my $aseqio = Bio::AnnSeqIO->new( -format => 'EMBL', 
					 -file => "efetch -a em:$embl_acc |");    
	my $cdna;
	eval {
	    $cdna = $aseqio->next_annseq();
	};

	# skip if we can't load it or no cdna.
	if( $@ || !$cdna) {
	    $log->store("Unable to read $embl_acc. Skipping. (Apologies)\n$@\n");
	    next ACC;
	}

	# get out CDS coordinates.
	
	my $cds_start;
	my $cds_end;
	
	foreach my $sf ( $cdna->top_SeqFeatures ) {
	    if( $sf->primary_tag() eq 'CDS' ) {
		$cds_start = $sf->start;
		$cds_end   = $sf->end;
	    }
	}
	
	if( !defined $cds_start ) {
	    print STDERR "EMBL $embl_acc has not single CDS line. Skipping\n";
	}
{'fraction_non_alphanumeric': 0.1549112426035503, 'fraction_numerical': 0.007928994082840236, 'mean_word_length': 3.266027258960121, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 49, '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': '4001162', 'n_tokens_mistral': 3453, 'n_tokens_neox': 3208, 'n_words': 913}