rsem v1.1.12, add a script to extract transcript-to-gene-map info from Trinity output

[rsem.git] / extract-transcript-to-gene-map-from-trinity

diff --git a/extract-transcript-to-gene-map-from-trinity b/extract-transcript-to-gene-map-from-trinity
new file mode 100755 (executable)
index 0000000..e04e993
--- /dev/null
+++ b/extract-transcript-to-gene-map-from-trinity
@@ -0,0 +1,34 @@
+#!/usr/bin/perl
+
+use strict;
+
+if (scalar(@ARGV) != 2) {
+    print "Usage: extract-transcript-to-gene-map-from-trinity trinity_fasta_file map_file\n";
+    exit(-1);
+}
+
+open(INPUT, $ARGV[0]);
+open(OUTPUT, ">$ARGV[1]");
+
+my ($tag, $line);
+$tag = <INPUT>; chomp($tag);
+while (substr($tag, 0, 1) eq ">") {
+    $tag = substr($tag, 1);
+    my $cnt = 0;
+    while (($line = <INPUT>) && substr($line, 0, 1) ne ">") {
+       $cnt++;
+    }
+    if ($cnt == 0) { print "Warning: Fasta entry $tag has an empty sequence, it is omitted.\n"; }
+    else {
+       my ($tid, $gid) = split(/ /, $tag);
+       my ($comp, $c, @tmp) = split(/_/, $gid);
+       $gid = join("_", $comp, $c);
+       print OUTPUT "$gid\t$tid\n";
+    }
+    $tag = $line; chomp($tag);
+} 
+
+close(INPUT);
+close(OUTPUT);
+
+