rsem-plot-transcript-wiggles

   1 #!/usr/bin/env perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5
   6 use FindBin;
   7 use lib $FindBin::RealBin;
   8 use rsem_perl_utils;
   9
  10 use Env qw(@PATH);
  11 @PATH = ($FindBin::RealBin, "$FindBin::RealBin/sam", @PATH);
  12
  13 use strict;
  14
  15
  16 my $gene_list = 0; # default is 0, means input is not a gene list
  17 my $transcript_list = 0; # default is 0, this option can only be turned on if allele-specific expression is calculated
  18 my $show_unique = 0; # 0, default value, means do not show unique transcript wiggles; 1 means show unique transcript wiggles
  19 my $help = 0;
  20
  21 GetOptions("gene-list" => \$gene_list,
  22            "transcript-list" => \$transcript_list,
  23            "show-unique" => \$show_unique,
  24            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  25
  26 pod2usage(-verbose => 2) if ($help == 1);
  27 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 3);
  28
  29 my $alleleS = 0;
  30 if (-e "$ARGV[0].alleles.results") { $alleleS = 1; }
  31
  32 pod2usage(-msg => "--transcript-list cannot be set if allele-specific reference is not built!", -exitval => 2, -verbose => 2) if (!$alleleS && $transcript_list);
  33 pod2usage(-msg => "--gene-list and --transcript-list cannot be set at the same time!", -exitval => 2, -verbose => 2) if ($gene_list && $transcript_list);
  34
  35 my $command = "";
  36
  37 unless (-e "$ARGV[0].transcript.sorted.bam") {
  38     $command = "samtools sort $ARGV[0].transcript.bam $ARGV[0].transcript.sorted";
  39     &runCommand($command);
  40 }
  41 unless (-e "$ARGV[0].transcript.readdepth") {
  42     $command = "rsem-bam2readdepth $ARGV[0].transcript.sorted.bam $ARGV[0].transcript.readdepth";
  43     &runCommand($command);
  44 }
  45
  46 if ($show_unique) {
  47     unless (-e "$ARGV[0].uniq.transcript.bam") {
  48         $command = "rsem-get-unique $ARGV[0].transcript.bam $ARGV[0].uniq.transcript.bam";
  49         &runCommand($command);
  50     }
  51     unless (-e "$ARGV[0].uniq.transcript.sorted.bam") {
  52         $command = "samtools sort $ARGV[0].uniq.transcript.bam $ARGV[0].uniq.transcript.sorted";
  53         &runCommand($command);
  54     }
  55     unless (-e "$ARGV[0].uniq.transcript.readdepth") {
  56         $command = "rsem-bam2readdepth $ARGV[0].uniq.transcript.sorted.bam $ARGV[0].uniq.transcript.readdepth";
  57         &runCommand($command);
  58     }
  59 }
  60
  61 my $id_type;
  62
  63 if ($alleleS) {
  64     $id_type = 0;
  65     if ($transcript_list) { $id_type = 1; }
  66     if ($gene_list) { $id_type = 2; }
  67 }
  68 else {
  69     $id_type = 1;
  70     if ($gene_list) { $id_type = 2; }
  71 }
  72
  73 $command = "rsem-gen-transcript-plots $ARGV[0] $ARGV[1] $alleleS $id_type $show_unique $ARGV[2]";
  74 &runCommand($command);
  75
  76 __END__
  77
  78 =head1 NAME
  79
  80 rsem-plot-transcript-wiggles
  81
  82 =head1 SYNOPSIS
  83
  84 rsem-plot-transcript-wiggles [options] sample_name input_list output_plot_file
  85
  86 =head1 ARGUMENTS
  87
  88 =over
  89
  90 =item B<sample_name>
  91
  92 The name of the sample analyzed.
  93
  94 =item B<input_list>
  95
  96 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.
  97
  98 =item B<output_plot_file>
  99
 100 The file name of the pdf file which contains all plots.
 101
 102 =back
 103
 104 =head1 OPTIONS
 105
 106 =over
 107
 108 =item B<--gene-list>
 109
 110 The input-list is a list of gene ids. (Default: off)
 111
 112 =item B<--transcript-list>
 113
 114 The input-list is a list of transcript ids. This option can only be turned on if allele-specific expression is calculated. (Default: off)
 115
 116 =item B<--show-unique>
 117
 118 Show the wiggle plots as stacked bar plots. See description section for details. (Default: off)
 119
 120 =item B<-h/--help>
 121
 122 Show help information.
 123
 124 =back
 125
 126 =head1 DESCRIPTION
 127
 128 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. If allele-specific expression is calculated, the basic unit for plotting is an allele-specific transcript and plots can be grouped by either transcript ids (--transcript-list) or gene ids (--gene-list).
 129
 130 =head1 OUTPUT
 131
 132 =over
 133
 134 =item B<output_plot_file>
 135
 136 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. If allele-specific expression is calculated, the basin unit becomes an allele-specific transcript and transcript ids and gene ids can be used to group allele-specific transcripts.
 137
 138 =item B<sample_name.transcript.sorted.bam and sample_name.transcript.readdepth>
 139
 140 If these files do not exist, 'rsem-plot-transcript-wiggles' will automatically generate them.
 141
 142 =item B<sample_name.uniq.transcript.bam, sample_name.uniq.transcript.sorted.bam and sample_name.uniq.transcript.readdepth>
 143
 144 If '--show-unique' option is specified and these files do not exist, 'rsem-plot-transcript-wiggles' will automatically generate them.
 145
 146 =back
 147
 148 =head1 EXAMPLES
 149
 150 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.
 151
 152 1) Transcript ids are provided and we just want normal wiggle plots:
 153
 154  rsem-plot-transcript-wiggles mmliver_single_quals transcript_ids.txt output.pdf
 155
 156 2) Gene ids are provided and we want to show stacked bar plots:
 157
 158  rsem-plot-transcript-wiggles --gene-list --show-unique mmliver_single_quals gene_ids.txt output.pdf
 159
 160 =cut