rsem-plot-transcript-wiggles

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use File::Basename;
   6 use strict;
   7
   8 my $gene_list = 0; # default is 0, means input is a transcript list; 1 means input is a gene list
   9 my $show_unique = 0; # 0, default value, means do not show unique transcript wiggles; 1 means show unique transcript wiggles
  10 my $help = 0;
  11
  12 GetOptions("gene-list" => \$gene_list,
  13            "show-unique" => \$show_unique,
  14            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  15
  16 pod2usage(-verbose => 2) if ($help == 1);
  17 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 3);
  18
  19 my ($fn, $dir, $suf) = fileparse($0);
  20 my $command = "";
  21
  22 unless (-e "$ARGV[0].transcript.sorted.bam") {
  23     $command = $dir."sam/samtools sort $ARGV[0].transcript.bam $ARGV[0].transcript.sorted";
  24     &runCommand($command);
  25 }
  26 unless (-e "$ARGV[0].transcript.readdepth") {
  27     $command = $dir."rsem-bam2readdepth $ARGV[0].transcript.sorted.bam $ARGV[0].transcript.readdepth";
  28     &runCommand($command);
  29 }
  30
  31 if ($show_unique) {
  32     unless (-e "$ARGV[0].uniq.transcript.bam") {
  33         $command = $dir."rsem-get-unique $ARGV[0].transcript.bam $ARGV[0].uniq.transcript.bam";
  34         &runCommand($command);
  35     }
  36     unless (-e "$ARGV[0].uniq.transcript.sorted.bam") {
  37         $command = $dir."sam/samtools sort $ARGV[0].uniq.transcript.bam $ARGV[0].uniq.transcript.sorted";
  38         &runCommand($command);
  39     }
  40     unless (-e "$ARGV[0].uniq.transcript.readdepth") {
  41         $command = $dir."rsem-bam2readdepth $ARGV[0].uniq.transcript.sorted.bam $ARGV[0].uniq.transcript.readdepth";
  42         &runCommand($command);
  43     }
  44 }
  45
  46 $command = $dir."rsem-gen-transcript-plots $ARGV[0] $ARGV[1] $gene_list $show_unique $ARGV[2]";
  47 &runCommand($command);
  48
  49 # command, {err_msg}
  50 sub runCommand {
  51     print $_[0]."\n";
  52     my $status = system($_[0]);
  53     if ($status != 0) {
  54         my $errmsg;
  55         if (scalar(@_) > 1) { $errmsg = $_[1]; }
  56         else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; }
  57         print $errmsg."\n";
  58         exit(-1);
  59     }
  60     print "\n";
  61 }
  62
  63 __END__
  64
  65 =head1 NAME
  66
  67 rsem-plot-transcript-wiggles
  68
  69 =head1 SYNOPSIS
  70
  71 =over
  72
  73  rsem-plot-transcript-wiggles [options] sample_name input_list output_plot_file
  74
  75 =back
  76
  77 =head1 ARGUMENTS
  78
  79 =over
  80
  81 =item B<sample_name>
  82
  83 The name of the sample analyzed.
  84
  85 =item B<input_list>
  86
  87 A list of transcript ids or gene ids. But it cannot be a mixture of transcript & gene ids. Each id occupies one line without extra spaces.
  88
  89 =item B<output_plot_file>
  90
  91 The file name of the pdf file which contains all plots.
  92
  93 =back
  94
  95 =head1 OPTIONS
  96
  97 =over
  98
  99 =item B<--gene-list>
 100
 101 The input-list is a list of gene ids. (Default: off)
 102
 103 =item B<--show-unique>
 104
 105 Show the wiggle plots as stacked bar plots. See description section for details. (Default: off)
 106
 107 =item B<-h/--help>
 108
 109 Show help information.
 110
 111 =back
 112
 113 =head1 DESCRIPTION
 114
 115 This program generates transcript wiggle plots and outputs them in a pdf file. This program can accept either a list of transcript ids or gene ids (if transcript to gene mapping information is provided) and has two modes of showing wiggle plots. If '--show-unique' is not specified, the wiggle plot for each transcript is a histogram where each position has the expected read depth at this position as its height. If '--show-unique' is specified, for each transcript a stacked bar plot is generated. For each position, the read depth of unique reads, which have only one alignment, is showed in black. The read depth of multi-reads, which align to more than one places, is showed in red on top of the read depth of unique reads.This program will use some files RSEM generated previouslly. So please do not delete/move any file 'rsem-calculate-expression' generated.
 116
 117 =head1 OUTPUT
 118
 119 =over
 120
 121 =item B<output_plot_file>
 122
 123 This is a pdf file containing all plots generated. If a list of transcript ids is provided, each page display at most 6 plots in 3 rows and 2 columns. If gene ids are provided, each page display a gene. The gene's id is showed at the top and all its transcripts' wiggle plots are showed in this page. The arrangment of plots is determined automatically. For each transcript wiggle plot, the transcript id is displayed as title. x-axis is position in the transcript and y-axis is read depth.
 124
 125 =item B<sample_name.transcript.sorted.bam and sample_name.transcript.readdepth>
 126
 127 If these files do not exist, 'rsem-plot-transcript-wiggles' will automatically generate them.
 128
 129 =item B<sample_name.uniq.transcript.bam, sample_name.uniq.transcript.sorted.bam and sample_name.uniq.transcript.readdepth>
 130
 131 If '--show-unique' option is specified and these files do not exist, 'rsem-plot-transcript-wiggles' will automatically generate them.
 132
 133 =back
 134
 135 =head1 EXAMPLES
 136
 137 Suppose sample_name and output_plot_file are set to 'mmliver_single_quals' and 'output.pdf' respectively. input_list is set to 'transcript_ids.txt' if transcript ids are provided, and is set to 'gene_ids.txt' if gene ids are provided.
 138
 139 1) Transcript ids are provided and we just want normal wiggle plots:
 140
 141  rsem-plot-transcript-wiggles mmliver_single_quals transcript_ids.txt output.pdf
 142
 143 2) Gene ids are provided and we want to show stacked bar plots:
 144
 145  rsem-plot-transcript-wiggles --gene-list --show-unique mmliver_single_quals gene_ids.txt output.pdf
 146
 147 =cut