rsem-prepare-reference

   1 #!/usr/bin/perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use File::Basename;
   6 use strict;
   7
   8 my $status;
   9
  10 my $gtfF = "";
  11 my $mappingF = "";
  12 my $polyAChoice = 0; # 0, default, pad polyA tails for all isoforms; 1, --no-polyA; 2, --no-polyA-subset
  13 my $no_polyA = 0; # for option --no-polyA
  14 my $subsetFile = "";
  15 my $polyALen = 125;
  16 my $bowtie_path = "";
  17 my $no_bowtie = 0;
  18 my $no_ntog = 0;
  19 my $quiet = 0;
  20 my $help = 0;
  21
  22 GetOptions("gtf=s" => \$gtfF,
  23            "transcript-to-gene-map=s" => \$mappingF,
  24            "no-polyA" => \$no_polyA,
  25            "no-polyA-subset=s" => \$subsetFile,
  26            "polyA-length=i" => \$polyALen,
  27            "bowtie-path=s" => \$bowtie_path,
  28            "no-bowtie" => \$no_bowtie,
  29            "no-ntog" => \$no_ntog,
  30            "q|quiet" => \$quiet,
  31            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  32
  33 pod2usage(-verbose => 2) if ($help == 1);
  34 pod2usage(-msg => "Set --no-polyA & --no-polyA-subset at the same time!", -exitval => 2, -verbose => 2) if ($no_polyA == 1 && $subsetFile ne '');
  35 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 2);
  36 pod2usage(-msg => "If bowtie is used, --no-ntog cannot be set!", -exitval => 2, -verbose => 2) if (!$no_bowtie && $no_ntog);
  37
  38 if ($no_bowtie && $bowtie_path ne "") { print "Warning: If bowtie is not used, no need to set --bowtie-path option!\n"; }
  39
  40 my $type;
  41
  42 if ($gtfF ne "") { $type = 0; }
  43 else { $type = 1; }
  44
  45 my @list = split(/,/, $ARGV[0]);
  46 my $size = scalar(@list);
  47
  48 if ($size == 1 && (-d $list[0])) {
  49     my $dir = $list[0];
  50     @list = (<$dir/*.fa>, <$dir/*.fasta>);
  51     $size = scalar(@list);
  52 }
  53
  54 if ($no_polyA) { $polyAChoice = 1 }
  55 elsif ($subsetFile ne "") { $polyAChoice = 2; }
  56
  57 if ($bowtie_path ne "") { $bowtie_path .= "/"; }
  58
  59 my ($fn, $dir, $suf) = fileparse($0);
  60 my $command = "";
  61
  62 if ($type == 0) {
  63     $"=" ";
  64     $command = $dir."rsem-extract-reference-transcripts $ARGV[1] $quiet $gtfF";
  65     if ($mappingF ne "") { $command .= " 1 $mappingF"; }
  66     else { $command .= " 0"; }
  67     $command .= " @list";
  68     print "$command\n";
  69     $status = system($command);
  70     if ($status != 0) {
  71         print "rsem-extract-reference-transcripts failed! Please check if you provide correct parameters/options for the pipeline!\n";
  72         exit(-1);
  73     }
  74     print "\n";
  75 }
  76 else {
  77     $"=" ";
  78     $command = $dir."rsem-synthesis-reference-transcripts $ARGV[1] $quiet";
  79     if ($mappingF ne "") { $command .= " 1 $mappingF"; }
  80     else { $command .= " 0"; }
  81     $command .= " @list";
  82     print "$command\n";
  83     $status = system($command);
  84     if ($status != 0) {
  85         print "rsem-synthesis-reference-transcripts failed! Please check if you provide correct parameters/options for the pipeline!\n";
  86         exit(-1);
  87     }
  88     print "\n";
  89 }
  90
  91 $command = $dir."rsem-preref $ARGV[1].transcripts.fa $polyAChoice $ARGV[1] -l $polyALen";
  92 if ($polyAChoice == 2) { $command .= " -f $subsetFile"; }
  93 if ($no_ntog) { $command .= " --no-ntog"; }
  94 if ($quiet) { $command .= " -q"; }
  95
  96 print "$command\n";
  97 $status = system($command);
  98 if ($status != 0) {
  99     print "rsem-preref failed! Please check if you provide correct parameters/options for the pipeline!\n";
 100     exit(-1);
 101 }
 102 print "\n";
 103
 104 if (!$no_bowtie) {
 105     $command = $bowtie_path."bowtie-build -f";
 106     if ($quiet) { $command .= " -q"; }
 107     $command .= " $ARGV[1].idx.fa $ARGV[1]";
 108
 109     print "$command\n";
 110     $status = system($command);
 111     if ($status != 0) {
 112         print "bowtie-build failed! Please check if you provide correct parameters/options for the pipeline!\n";
 113         exit(-1);
 114     }
 115     print "\n";
 116 }
 117
 118 __END__
 119
 120 =head1 NAME
 121
 122 rsem-prepare-reference
 123
 124 =head1 SYNOPSIS
 125
 126 =over
 127
 128  rsem-prepare-reference [options] reference_fasta_file(s) reference_name
 129
 130 =back
 131
 132 =head1 ARGUMENTS
 133
 134 =over
 135
 136 =item B<reference_fasta_file(s)>
 137
 138 Either a comma-separated list of FASTA formatted files OR a directory name. If a directory name is specified, RSEM will read all files with suffix ".fa" or ".fasta" in this directory.  The files should contain either the sequences of transcripts or an entire genome, depending on whether the --gtf option is used.
 139
 140 =item B<reference name>
 141
 142 The name of the reference used. RSEM will generate several reference-related files that are prefixed by this name. This name can contain path information (e.g. /ref/mm9).
 143
 144 =back
 145
 146 =head1 OPTIONS
 147
 148 =over
 149
 150 =item B<--gtf> <file>
 151
 152 If this option is on, RSEM assumes that 'reference_fasta_file(s)' contains the sequence of a genome, and will extract transcript reference sequences using the gene annotations specified in <file>, which should be in GTF format.
 153
 154 If this option is off, RSEM will assume 'reference_fasta_file(s)' contains the reference transcripts. In this case, RSEM assumes that name of each sequence in the FASTA files is its transcript_id.
 155
 156 (Default: off)
 157
 158 =item B<--transcript-to-gene-map> <file>
 159
 160 Use information from <file> to map from transcript (isoform) ids to gene ids.
 161 Each line of <file> should be of the form:
 162
 163 gene_id transcript_id
 164
 165 with the two fields separated by a tab character.
 166
 167 If you are using a GTF file for the "UCSC Genes" gene set from the UCSC Genome Browser, then the "knownIsoforms.txt" file (obtained from the "Downloads" section of the UCSC Genome Browser site) is of this format.
 168
 169 If this option is off, then the mapping of isoforms to genes depends on whether the --gtf option is specified.  If --gtf is specified, then RSEM uses the "gene_id" and "transcript_id" attributes in the GTF file.  Otherwise, RSEM assumes that each sequence in the reference sequence files is a separate gene.
 170
 171 (Default: off)
 172
 173 =item B<--no-polyA>
 174
 175 Do not add poly(A) tails to the end of reference isoforms. (Default: add poly(A) tails to all transcripts)
 176
 177 =item B<--no-polyA-subset> <file>
 178
 179 Add poly(A) tails to all transcripts except those listed in <file>. <file> is a file containing a list of transcript_ids. (Default: add poly(A) tails to all transcripts. This option cannot be used with '--no-polyA'.)
 180
 181 =item B<--polyA-length> <int>
 182
 183 The length of the poly(A) tails to be added. (Default: 125)
 184
 185 =item B<--bowtie-path> <path>
 186
 187 The path to the bowtie executables. (Default: the path to bowtie executables is assumed to be in the user's PATH environment variable)
 188
 189 =item B<--no-bowtie>
 190
 191 Do not build Bowtie indices.  Specify this option if you wish to use an alternative aligner for mapping reads to transcripts.  You should align against the sequences generated in the output file 'reference_name.idx.fa'.  (Default: off)
 192
 193 =item B<--no-ntog>
 194
 195 Disable the conversion of 'N' characters to 'G' characters in the reference sequences.  This conversion is normally desired because it allows some aligners (e.g., Bowtie) to align against all positions in the reference.
 196 (Default: off)
 197
 198 =item B<-q/--quiet>
 199
 200 Suppress the output of logging information. (Default: off)
 201
 202 =item B<-h/--help>
 203
 204 Show help information.
 205
 206 =back
 207
 208 =head1 DESCRIPTION
 209
 210 This program extracts/preprocesses the reference sequences and builds Bowtie indices using default parameters.  If an alternative aligner is to be used, then the '--no-bowtie' option should be specified to disable building Bowtie indices. This program is used in conjunction with the 'rsem-calculate-expression' program.
 211
 212 =head1 OUTPUT
 213
 214 This program will generate 'reference_name.grp', 'reference_name.ti', 'reference_name.transcripts.fa', 'reference_name.seq', 'reference_name.chrlist' (if '--gtf' is on), 'reference_name.idx.fa', and corresponding Bowtie index files (unless '--no-bowtie' is specified).
 215
 216 'reference_name.grp', 'reference_name.ti', 'reference_name.seq', 'reference_name.idx.fa', and 'reference_name.chrlist' are used by RSEM internally.
 217
 218 B<'reference_name.transcripts.fa'> contains the extracted reference transcripts in FASTA format. Poly(A) tails are not added.
 219
 220 =head1 EXAMPLES
 221
 222 1) Suppose we have mouse RNA-Seq data and want to use the UCSC mm9 version of the mouse genome. We have downloaded the UCSC Genes transcript annotations in GTF format (as mm9.gtf) using the Table Browser and the knownIsoforms.txt file for mm9 from the UCSC Downloads. We also have all chromosome files for mm9 in the directory '/data/mm9'.  We want to put the generated reference files under '/ref' with name 'mm9'. We'll add poly(A) tails with length 125. Please note that GTF files generated from UCSC's Table Browser do not contain isoform-gene relationship information.  For the UCSC Genes annotation, this information can be obtained from the knownIsoforms.txt file.  Suppose we want to build Bowtie indices and Bowtie executables are found in '/sw/bowtie'.
 223
 224 There are two ways to write the command:
 225
 226  rsem-prepare-reference --gtf mm9.gtf \
 227                         --transcript-to-gene-map knownIsoforms.txt \
 228                         --bowtie-path /sw/bowtie \
 229                         /data/mm9/chr1.fa,/data/mm9/chr2.fa,...,/data/mm9/chrM.fa \
 230                         /ref/mm9
 231
 232 OR
 233
 234  rsem-prepare-reference --gtf mm9.gtf \
 235                         --transcript-to-gene-map knownIsoforms.txt \
 236                         --bowtie-path /sw/bowtie \
 237                         /data/mm9 \
 238                         /ref/mm9
 239
 240 2) Suppose we only have transcripts from EST tags in 'mm9.fasta'. In addition, we also have isoform-gene information in 'mapping.txt'. We do not want to add any poly(A) tails. The reference_name will be set to 'mm9'. In addition, we do not want to build Bowtie indices, and will use an alternative aligner to align reads against the 'mm9.idx.fa' output file:
 241
 242  rsem-prepare-reference --transcript-to-gene-map mapping.txt \
 243                         --no-polyA \
 244                         --no-bowtie \
 245                         mm9.fasta \
 246                         mm9
 247
 248 =cut