#!/usr/bin/perl -w
use strict;

#Indentify solo LTRs in RM out files
#Shujun Ou (oushujun@msu.edu)
#12-14-2020

my $usage="perl this_script.pl -i RepeatMasker.out -info lib.LTR.info > solo_list";
my $SWscore=300; #minumum Smith-Waterman alignment score
my $minlen=80; #bp, shortest solo LTR length allowed
my $mincov=0.8; #minimum coverage of the library LTR sequence
my $maxcov = 1.2; #maximum coverage
my $mindist=300; #bp, shortest distance between candidate LTR to adjacent LTR internal seq annotation allowed. (not consider as solo if LTR-INT is too close)
my $RMout = '';
my $info = ''; #an extended bed file containing coordinates of all LTR regions in the library used to generate $RMout. Generated by find_LTR.pl

my $i=0;
foreach (@ARGV){
	$SWscore = $ARGV[$i+1] if /^-minsw$/i;
	$minlen = $ARGV[$i+1] if /^-minlen$/i;
	$mincov = $ARGV[$i+1] if /^-mincov$/i;
	$maxcov = $ARGV[$i+1] if /^-maxcov$/i;
	$mindist = $ARGV[$i+1] if /^-mindist$/i;
	$RMout = $ARGV[$i+1] if /^-i$/i;
	$info = $ARGV[$i+1] if /^-info$/i;
	$i++;
	}

open RM, "<$RMout" or die $usage;
open Info, "<$info" or die $usage;

my %info;
while (<Info>){
	my ($id, $srt, $end, $len) = (split);
	$id =~ s/#.*//;
	push @{$info{$id}}, [$srt, $end, $len];
	}
close Info;

my @arr;
while (<RM>){
	chomp;
	s/^\s+//g;
	s/[\(\)]+//g;
	next unless /^[0-9]+/;
	my ($score, $start, $end, $dir, $loc, $supfam, $libs, $libe, $libr)=(split)[0,5,6,8,9,10,11,12,13];
	next if $supfam !~ /LTR|RLX|RLC|RLG/i;
	next unless exists $info{$loc};
	next if $score<=$SWscore or $end-$start<=$minlen; #require alignment score>300 or alignment length>100bp
	($libs, $libr) = ($libr, $libs) if $dir =~ "C";
	# determine if this is an LTR region
	my $is_ltr = 0;
#	print "$_\n";
	$is_ltr = find_ltr ($info{$loc}, $libs, $libe);
	push @arr, [$_, $is_ltr];
	if (@arr>5){
		my ($info2, $is_ltr2) = ($arr[2]->[0], $arr[2]->[1]);
		my ($div2, $chr2, $s2, $e2, $dir2, $loc2, $fam2, $libs2, $libe2, $libr2)=(split /\s+/, $info2)[1,4,5,6,8,9,10,11,12,13];
		my $loc2_len=@{$info{$loc2}}[0]->[2];
		my $cov = abs($e2-$s2)/$loc2_len;

		my $keep=1;
		$keep=0 if $is_ltr2 == 0; #skip this line if this is not an LTR region
		$keep=0 if $cov > $maxcov; #skip this line if the lib sequence span too much on the genome
		foreach ($arr[0], $arr[1], $arr[3], $arr[4]){
			my ($info0, $is_ltr0) = ($_->[0], $_->[1]);
			my ($div0, $chr0, $s0, $e0, $dir0, $loc0, $fam0, $libs0, $libe0, $libr0)=(split /\s+/, $info0)[1,4,5,6,8,9,10,11,12,13];
			$keep=0 if ($chr0 eq $chr2) and ($dir0 eq $dir2) and ($loc0 eq $loc2) and (abs($div2-$div0)<4); #require difference of diversity larger than 4% compared to other LTRs with the same ID matched to be a solo-LTR
			$keep=0 if ($chr0 eq $chr2) and ($is_ltr0 == 0) and (abs($s2-$e0)<$mindist) and (abs($s0-$e2)<$mindist); #if there is int exists in this working window, require it has at least 300bp distance from a solo-LTR
			}
		print "$chr2\t$s2\t$e2\t$chr2:$s2..$e2\t$loc2\t$cov\n" if $keep==1;  #print the actual locus of the solo LTR
		shift @arr;
		}
	}
close RM;


# determine if this is an LTR region
sub find_ltr {
	my ($info, $libs, $libe) = ($_[0], $_[1], $_[2]);
	my @info = @{$info};
	($libs, $libe) = ($libe, $libs) if $libe < $libs;
	my $match_len = $libe - $libs + 1;
	my $is_ltr = 0;
	foreach my $ltr (@info){
		my ($libstr, $libend, $liblen) = ($ltr->[0], $ltr->[1], $ltr->[2]);
		($libstr, $libend) = ($libend, $libstr) if $libstr > $libend;
		my $ltrlen = $libend - $libstr + 1;
		if ($libe <= $libend and $libs >= $libstr){
			$is_ltr = 1 if $match_len/$ltrlen >= $mincov;
			} elsif ($libs < $libend and $libe > $libend){
			$is_ltr = 1 if ($libend - $libs + 1)/$ltrlen >= $mincov and $match_len/$ltrlen <= $maxcov;
			} elsif ($libs < $libstr and $libe > $libstr){
			$is_ltr = 1 if ($libe - $libstr + 1)/$ltrlen >= $mincov and $match_len/$ltrlen <= $maxcov;
			}
#print "$is_ltr, $libstr, $libend, $ltrlen, $liblen, $libs, $libe\n";
		}
	return $is_ltr;
	}
