bioseq: add --trim-ends

Author: wqiu

bin/bioseq

@@ -27,6 +27,10 @@ GetOptions(
     "count-codons|C",
     "delete|d=s",
     "feat2fas|F",
+    "fetch|f=s", # Retrieve sequence by accession number; broken by NCBI protocol change
+ #   "longest-orf|C",
+ #   "extract|e",
+ #   "dotplot|D=s",
     "hydroB|H",
     "iep",
     "input|i=s",
@@ -36,6 +40,7 @@ GetOptions(
     "linearize|L",
     "mol-wt",
     "no-gaps|g",
+    "no-revcom|Z",
     "num-gaps-dna",
     "num-gaps-aa",
     "num-seq|n",
@@ -45,21 +50,17 @@ GetOptions(
     "remove-stop|X", # for PAML/codeml
 #    "restrict|x=s",
 #    "restrict-coord=s",
+    "rename|N=s",
     "revcom|r",
     "sort=s",
     "split-cdhit=s",
     "subseq|s=s",
     "syn-code",
     "translate|t=i",
-    "fetch|f=s", # Retrieve sequence by accession number; broken by NCBI protocol change
- #   "longest-orf|C",
- #   "extract|e",
- #   "dotplot|D=s",
-    "rename|N=s",
-    "no-revcom|Z",
  #   "slidingwindow|S=i",
  #   "prefix=s",
- #   "split|S=i",
+    #   "split|S=i",
+    "trim-ends",
 ) or pod2usage(2);
 
 use Pod::Usage;

lib/Bio/BPWrapper/SeqManipulations.pm

@@ -51,7 +51,7 @@ pick_by_id del_by_id find_by_re
 pick_by_re del_by_re
 pick_by_file del_by_file
 find_by_ambig pick_by_ambig del_by_ambig find_by_length
-del_by_length codon_sim codon_info);
+del_by_length codon_sim codon_info trim_ends);
 
 # Package global variables
 my ($in, $out, $seq, %opts, $filename, $in_format, $out_format, $guesser);
@@ -99,7 +99,7 @@ my %opt_dispatch = (
 #	'prefix' => \&anonymize,
 	'rename' => \&rename_id,
 #	'slidingwindow' => \&sliding_window,
-#	'split' => \&split_seqs,
+	'trim-ends' => \&trim_ends,
   );
 
 my %filter_dispatch = (
@@ -193,6 +193,28 @@ sub write_out {
 
 ################### subroutines ########################
 
+# from vapid trim python code
+sub trim_ends {
+    while ($seq = $in->next_seq()) {
+	my $xip = 0;
+	my $len = $seq->length();
+	my @mono = split //, $seq->seq();
+	for(my $i = 0; $i < $len; $i++) {
+	    $xip ++ if $mono[$i] =~ /[NX\-\?]/i;
+	    last unless $mono[$i] =~ /[NX\-\?]/i;
+	}
+
+	my $y = $len;
+	for(my $i = $len - 1; $i >=0; $i--) {
+	    $y-- if $mono[$i] =~ /[NX\-\?]/i;
+	    last unless $mono[$i] =~ /[NX\-\?]/i;
+	}
+	$seq->seq($seq->subseq($xip+1, $y));
+	$out->write_seq($seq);
+    }
+}
+
+
 sub codon_sim {
 #    my $seq = $in->next_seq(); # only the first sequence used
 #    if (&_internal_stop_or_x($seq->translate()->seq())) {

test-files/test-bioseq-trim-ends.fas

@@ -0,0 +1,7 @@
+>N40:1:101:1716:ORFK00021
+-----?----------NN---nn-----atgtataaacaacaatattttatttctggc
+aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
+aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
+aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
+aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
+tattat-----?----------NN---NN