make sure candidates have sufficient flanking sequence to extended (5…

…0bp)

LTR_retriever

@@ -714,7 +714,7 @@ chomp ($date=`date`);
 print "$date\tModules 2-5: Start to analyze the structure of candidates...
 \t\t\t\tThe terminal motif, TSD, boundary, orientation, age, and superfamily will be identified in this step.\n\n";
 my $index=$_[0];
-`perl $script_path/bin/get_range.pl $index.retriever.scn $index.ltrTE.stg1 -x`;
+`perl $script_path/bin/get_range.pl $index.retriever.scn $index.ltrTE.stg1 -x -strict 1`;
 `cat $index.retriever.scn.extend |sort -fu > $index.retriever.scn.extend.unq; mv $index.retriever.scn.extend.unq $index.retriever.scn.extend`;
 `perl $script_path/bin/call_seq_by_list.pl $index.retriever.scn.extend -C $genome > $index.retriever.scn.extend.fa`; #full TE sequence with 50bp-extended on each side
 

bin/call_seq_by_list.pl

@@ -7,6 +7,7 @@
 eg2. perl call_seq_by_LOC.pl array_list -C your_database up_2000 >result	##call sequence of upper 2000 bp region in the list, from the provided database
 
 Update history:
+	v2.6	add -strict parameter to control inclusiveness of input range (2024/01/07)
 	v2.5	do not output sequence with all Ns (2023/01/30)
 	v2.4	output sequences without headers
 	v2.3	output a list of entirely excluded sequences
@@ -34,6 +35,7 @@
 my $exclude=0; #0 for output sequence specified by list (default); 1 for exclude sequence specified by list
 my $coverage=1; #work with $exclude, if the excluded portion is too long (default 1, [0-1]), discard the entire sequence
 my $purge=0; #work with $exclude, switch on=1/off=0(default) to clean up aligned region and joint unaligned sequences
+my $strict=0; #0=flexible, will adjust coordinate start and stop if they extend outside of the sequence; 1 will skip the entry if exceed.
 my $genome;
 
 my $i=0;
@@ -56,6 +58,7 @@
 	$header=$ARGV[$i+1] if $para=~/^-header$/i;
 	$coverage=$ARGV[$i+1] if $para=~/^-cov$/i;
 	$purge=$ARGV[$i+1] if $para=~/^-purge$/i;
+	$strict=$ARGV[$i+1] if $para=~/^-strict$/i;
 	$exclude=1 if $para=~/^-ex$/i;
 	$i++;
 	}
@@ -123,12 +126,17 @@
 		$stop=$start+$length-1;
 		}
 
-	$start=1 if $start<=0;
 
 	next unless exists $genome{$chr};
-#	next if $genome{$chr}=~/^\s+$/;
 	next if length $genome{$chr} < 10;
-	$stop=length $genome{$chr} if $stop>=length $genome{$chr};
+
+	# adjust boundary if they exceed the sequence range
+	if ($strict == 0){
+		$start=1 if $start<=0;
+		$stop=length $genome{$chr} if $stop>=length $genome{$chr};
+		} else {
+		next if $start<=0 or $stop>=length $genome{$chr};
+		}
 
 	if ($exclude==0){
 		$seq=substr ($genome{$chr}, $start-1, $stop-$start+1) if exists $genome{$chr};

bin/get_range.pl

@@ -160,7 +160,7 @@
 		$lLTR_end+=50;
 		$rLTR_start-=50;
 		$rLTR_end+=50;
-		print Extend "$chr:$element_start..$element_end\t$chr:$lLTR_start..$rLTR_end\n" if defined $chr;
+		print Extend "$chr:$element_start..$element_end\t$chr:$lLTR_start..$rLTR_end\n" if defined $chr and $lLTR_start > 0;
 		}
 
 	if ($full==1){