rsem-generate-ngvector

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use File::Basename;
   6 use strict;
   7
   8 my $k = 25;
   9 my $help = 0;
  10
  11 GetOptions("k=i" => \$k,
  12            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  13
  14 pod2usage(-verbose => 2) if ($help == 1);
  15 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 2);
  16
  17 my ($fn, $dir, $suf) = fileparse($0);
  18 my $command = "";
  19
  20 $command = $dir."rsem-for-ebseq-calculate-clustering-info $k $ARGV[0] $ARGV[1].ump";
  21 &runCommand($command);
  22
  23 $command = $dir."rsem-for-ebseq-generate-ngvector-from-clustering-info $ARGV[1].ump $ARGV[1].ngvec";
  24 &runCommand($command);
  25
  26 # command, {err_msg}
  27 sub runCommand {
  28     print $_[0]."\n";
  29     my $status = system($_[0]);
  30     if ($status != 0) {
  31         my $errmsg = "";
  32         if (scalar(@_) > 1) { $errmsg .= $_[1]."\n"; }
  33         $errmsg .= "\"$_[0]\" failed! Plase check if you provide correct parameters/options for the pipeline!\n";
  34         print $errmsg;
  35         exit(-1);
  36     }
  37     print "\n";
  38 }
  39
  40 __END__
  41
  42 =head1 NAME
  43
  44 rsem-generate-ngvector
  45
  46 =head1 SYNOPSIS
  47
  48 =over
  49
  50  rsem-generate-ngvector [options] input_fasta_file output_name
  51
  52 =back
  53
  54 =head1 ARGUMENTS
  55
  56 =over
  57
  58 =item B<input_fasta_file>
  59
  60 The fasta file containing all reference transcripts. The transcripts must be in the same order as those in expression value files. Thus, 'reference_name.transcripts.fa' generated by 'rsem-prepare-reference' should be used.
  61
  62 =item B<output_name>
  63
  64 The name of all output files. The Ng vector will be stored as 'output_name.ngvec'.
  65
  66 =back
  67
  68 =head1 OPTIONS
  69
  70 =over
  71
  72 =item B<-k> <int>
  73
  74 k mer length. See description section. (Default: 25)
  75
  76 =item B<-h/--help>
  77
  78 Show help information.
  79
  80 =back
  81
  82 =head1 DESCRIPTION
  83
  84 This program generates the Ng vector required by EBSeq for isoform level differential expression analysis based on reference sequences only. EBSeq can take variance due to read mapping ambiguity into consideration by grouping isoforms with parent gene's number of isoforms. However, for de novo assembled transcriptome, it is hard to obtain an accurate gene-isoform relationship. Instead, this program groups isoforms by using measures on read mappaing ambiguity directly. First, it calcualtes the 'unmappability' of each transcript. The 'unmappability' of a transcript is the ratio between the number of k mers with at least one perfect match to other transcripts and the total number of k mers of this transcript, where k is a parameter. Then, Ng vector is generated by applying Kmeans algorithm to the 'unmappability' values with number of clusters set as 3. 'rsem-generate-ngvector' will make sure the mean 'unmappability' scores for clusters are in ascending order. All transcripts whose lengths are less than k are assigned to cluster 3.
  85
  86 If your reference is a de novo assembled transcript set, you should run 'rsem-generate-ngvector' first. Then load the resulting 'output_name.ngvec' into R. For example, you can use
  87
  88  NgVec <- scan(file="output_name.ngvec", what=0, sep="\n")
  89
  90 . After that, replace 'IsoNgTrun' with 'NgVec' in the second line of section 3.2.5 (Page 10) of EBSeq's vignette:
  91
  92  IsoEBres=EBTest(Data=IsoMat, NgVector=NgVec, ...)
  93
  94 This program only needs to run once per RSEM reference.
  95
  96 =head1 OUTPUT
  97
  98 =over
  99
 100 =item B<output_name.ump>
 101
 102 'unmappability' scores for each transcript. This file contains two columns. The first column is transcript name and the second column is 'unmappability' score.
 103
 104 =item B<output_name.ngvec>
 105
 106 Ng vector generated by this program.
 107
 108 =back
 109
 110 =head1 EXAMPLES
 111
 112 Suppose the reference sequences file is '/ref/mouse_125/mouse_125.transcripts.fa' and we set the output_name as 'mouse_125':
 113
 114  rsem-generate-ngvector /ref/mouse_125/mouse_125.transcripts.fa mouse_125
 115
 116 =cut