rsem-prepare-reference

   1 #!/usr/bin/env perl
   2
   3 use Getopt::Long;
   4 use Pod::Usage;
   5 use FindBin;
   6 use lib $FindBin::RealBin;
   7 use rsem_perl_utils;
   8
   9 use Env qw(@PATH);
  10 push(@PATH, $FindBin::RealBin);
  11
  12 use strict;
  13
  14 my $status;
  15
  16 my $gtfF = "";
  17 my $mappingF = "";
  18 my $polyAChoice = 0; # 0, default, pad polyA tails for all isoforms; 1, --no-polyA; 2, --no-polyA-subset
  19 my $no_polyA = 0; # for option --no-polyA
  20 my $subsetFile = "";
  21 my $polyALen = 125;
  22 my $bowtie_path = "";
  23 my $no_bowtie = 0;
  24 my $no_ntog = 0;
  25 my $bowtie2 = 0;
  26 my $bowtie2_path = "";
  27 my $quiet = 0;
  28 my $help = 0;
  29
  30 my $alleleMappingF = "";
  31
  32 GetOptions("gtf=s" => \$gtfF,
  33            "transcript-to-gene-map=s" => \$mappingF,
  34            "allele-to-gene-map=s" => \$alleleMappingF,
  35            "no-polyA" => \$no_polyA,
  36            "no-polyA-subset=s" => \$subsetFile,
  37            "polyA-length=i" => \$polyALen,
  38            "bowtie-path=s" => \$bowtie_path,
  39            "no-bowtie" => \$no_bowtie,
  40            "no-ntog" => \$no_ntog,
  41            "bowtie2" => \$bowtie2,
  42            "bowtie2-path=s" => \$bowtie2_path,
  43            "q|quiet" => \$quiet,
  44            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  45
  46 pod2usage(-verbose => 2) if ($help == 1);
  47 pod2usage(-msg => "Set --no-polyA & --no-polyA-subset at the same time!", -exitval => 2, -verbose => 2) if ($no_polyA == 1 && $subsetFile ne '');
  48 pod2usage(-msg => "--transcript-to-gene-map and --allele-to-gene-map are mutually exclusive!", -exitval => 2, -verbose => 2) if (($mappingF ne "") && ($alleleMappingF ne ""));
  49 pod2usage(-msg => "--gtf and --allele-to-gene-map are mutually exclusive!", -exitval => 2, -verbose => 2) if (($gtfF ne "") && ($alleleMappingF ne ""));
  50 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 2);
  51
  52 if ($bowtie2) { $no_bowtie = 1; $no_ntog = 1; }
  53
  54 pod2usage(-msg => "If bowtie is used, --no-ntog cannot be set!", -exitval => 2, -verbose => 2) if (!$no_bowtie && $no_ntog);
  55
  56 if ($no_bowtie && ($bowtie_path ne "")) { print "Warning: If bowtie is not used, no need to set --bowtie-path option!\n"; }
  57 if (!$bowtie2 && ($bowtie2_path ne "")) { print "Warning: If bowtie2 is not used, no need to set --bowtie2-path option!\n"; }
  58
  59 my $type;
  60
  61 if ($gtfF ne "") { $type = 0; }
  62 else { $type = 1; }
  63
  64 my @list = split(/,/, $ARGV[0]);
  65 my $size = scalar(@list);
  66
  67 if ($size == 1 && (-d $list[0])) {
  68     my $dir = $list[0];
  69     @list = (<$dir/*.fa>, <$dir/*.fasta>);
  70     $size = scalar(@list);
  71 }
  72
  73 pod2usage(-msg => "reference_fasta_file(s) is empty! Please check if you provide the correct folder name or file suffixes!", -exitval => 2, -verbose => 2) if ($size <= 0);
  74
  75 if ($no_polyA) { $polyAChoice = 1 }
  76 elsif ($subsetFile ne "") { $polyAChoice = 2; }
  77
  78 if ($bowtie_path ne "") { $bowtie_path .= "/"; }
  79 if ($bowtie2_path ne "") { $bowtie2_path .= "/"; }
  80
  81 my $command = "";
  82
  83 if ($type == 0) {
  84     $"=" ";
  85     $command = "rsem-extract-reference-transcripts $ARGV[1] $quiet $gtfF";
  86     if ($mappingF ne "") { $command .= " 1 $mappingF"; }
  87     else { $command .= " 0"; }
  88     $command .= " @list";
  89     &runCommand($command);
  90 }
  91 else {
  92     $"=" ";
  93     $command = "rsem-synthesis-reference-transcripts $ARGV[1] $quiet";
  94     if ($mappingF ne "") { $command .= " 1 $mappingF"; }
  95     elsif ($alleleMappingF ne "") { $command .= " 2 $alleleMappingF"; }
  96     else { $command .= " 0"; }
  97     $command .= " @list";
  98     &runCommand($command);
  99 }
 100
 101 $command = "rsem-preref $ARGV[1].transcripts.fa $polyAChoice $ARGV[1] -l $polyALen";
 102 if ($polyAChoice == 2) { $command .= " -f $subsetFile"; }
 103 if ($no_ntog) { $command .= " --no-ntog"; }
 104 if ($quiet) { $command .= " -q"; }
 105
 106 &runCommand($command);
 107
 108 if (!$no_bowtie) {
 109     $command = $bowtie_path."bowtie-build -f";
 110     if ($quiet) { $command .= " -q"; }
 111     $command .= " $ARGV[1].idx.fa $ARGV[1]";
 112
 113     &runCommand($command);
 114 }
 115
 116 if ($bowtie2) {
 117     $command = $bowtie2_path."bowtie2-build -f";
 118     if ($quiet) { $command .= " -q"; }
 119     $command .= " $ARGV[1].idx.fa $ARGV[1]";
 120
 121     &runCommand($command);
 122 }
 123
 124 __END__
 125
 126 =head1 NAME
 127
 128 rsem-prepare-reference
 129
 130 =head1 SYNOPSIS
 131
 132 rsem-prepare-reference [options] reference_fasta_file(s) reference_name
 133
 134 =head1 ARGUMENTS
 135
 136 =over
 137
 138 =item B<reference_fasta_file(s)>
 139
 140 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.
 141
 142 =item B<reference name>
 143
 144 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).
 145
 146 =back
 147
 148 =head1 OPTIONS
 149
 150 =over
 151
 152 =item B<--gtf> <file>
 153
 154 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.
 155
 156 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.
 157
 158 (Default: off)
 159
 160 =item B<--transcript-to-gene-map> <file>
 161
 162 Use information from <file> to map from transcript (isoform) ids to gene ids.
 163 Each line of <file> should be of the form:
 164
 165 gene_id transcript_id
 166
 167 with the two fields separated by a tab character.
 168
 169 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.
 170
 171 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.
 172
 173 (Default: off)
 174
 175 =item B<--allele-to-gene-map> <file>
 176
 177 Use information from <file> to provide gene_id and transcript_id information for each allele-specific transcript.
 178 Each line of <file> should be of the form:
 179
 180 gene_id transcript_id allele_id
 181
 182 with the fields separated by a tab character.
 183
 184 This option is designed for quantifying allele-specific expression. It is only valid if '--gtf' option is not specified. allele_id should be the sequence names presented in the FASTA-formatted files.
 185
 186 (Default: off)
 187
 188 =item B<--no-polyA>
 189
 190 Do not add poly(A) tails to the end of reference isoforms. (Default: adding poly(A) tails to all transcripts)
 191
 192 =item B<--no-polyA-subset> <file>
 193
 194 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'.)
 195
 196 =item B<--polyA-length> <int>
 197
 198 The length of the poly(A) tails to be added. (Default: 125)
 199
 200 =item B<--bowtie-path> <path>
 201
 202 The path to the Bowtie executables. (Default: the path to Bowtie executables is assumed to be in the user's PATH environment variable)
 203
 204 =item B<--no-bowtie>
 205
 206 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)
 207
 208 =item B<--no-ntog>
 209
 210 Disable the conversion of 'N' characters to 'G' characters in the reference sequences prepared for aligners.  This conversion is in particular desired for Bowtie aligner to align against all positions in the reference. (Default: off)
 211
 212 =item B<--bowtie2>
 213
 214 Build Bowtie 2 indices instead of Bowtie indices. Turn on this option will automatically turn on '--no-bowtie' and '--no-ntog' options. (Default: off)
 215
 216 =item B<--bowtie2-path>
 217
 218 The path to the Bowtie 2 executables. (Default: the path to Bowtie 2 executables is assumed to be in the user's PATH environment variable)
 219
 220 =item B<-q/--quiet>
 221
 222 Suppress the output of logging information. (Default: off)
 223
 224 =item B<-h/--help>
 225
 226 Show help information.
 227
 228 =back
 229
 230 =head1 DESCRIPTION
 231
 232 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.
 233
 234 =head1 OUTPUT
 235
 236 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).
 237
 238 'reference_name.grp', 'reference_name.ti', 'reference_name.seq', and 'reference_name.chrlist' are used by RSEM internally.
 239
 240 B<'reference_name.transcripts.fa'> contains the extracted reference transcripts in FASTA format. Poly(A) tails are added unless '--no-polyA' is set.
 241
 242 B<'reference_name.idx.fa'> is used by aligners to build their own indices. If '--no-ntog' is set, this file should be identical to 'reference_name.transcripts.fa'.
 243
 244 =head1 EXAMPLES
 245
 246 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 'mouse_125'. 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'.
 247
 248 There are two ways to write the command:
 249
 250  rsem-prepare-reference --gtf mm9.gtf \
 251                         --transcript-to-gene-map knownIsoforms.txt \
 252                         --bowtie-path /sw/bowtie \
 253                         /data/mm9/chr1.fa,/data/mm9/chr2.fa,...,/data/mm9/chrM.fa \
 254                         /ref/mouse_125
 255
 256 OR
 257
 258  rsem-prepare-reference --gtf mm9.gtf \
 259                         --transcript-to-gene-map knownIsoforms.txt \
 260                         --bowtie-path /sw/bowtie \
 261                         /data/mm9 \
 262                         /ref/mouse_125
 263
 264 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 'mouse_0'. In addition, we do not want to build Bowtie indices, and will use an alternative aligner to align reads against the 'mouse_0.idx.fa' output file:
 265
 266  rsem-prepare-reference --transcript-to-gene-map mapping.txt \
 267                         --no-polyA \
 268                         --no-bowtie \
 269                         mm9.fasta \
 270                         mouse_0
 271
 272 =cut