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 @PATH = ($FindBin::RealBin, @PATH);
  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 = 0;
  23 my $bowtie_path = "";
  24 my $bowtie2 = 0;
  25 my $bowtie2_path = "";
  26 my $quiet = 0;
  27 my $help = 0;
  28
  29 my $alleleMappingF = "";
  30
  31 GetOptions("gtf=s" => \$gtfF,
  32            "transcript-to-gene-map=s" => \$mappingF,
  33            "allele-to-gene-map=s" => \$alleleMappingF,
  34            "no-polyA" => \$no_polyA,
  35            "no-polyA-subset=s" => \$subsetFile,
  36            "polyA-length=i" => \$polyALen,
  37            "bowtie" => \$bowtie,
  38            "bowtie-path=s" => \$bowtie_path,
  39            "bowtie2" => \$bowtie2,
  40            "bowtie2-path=s" => \$bowtie2_path,
  41            "q|quiet" => \$quiet,
  42            "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
  43
  44 pod2usage(-verbose => 2) if ($help == 1);
  45 pod2usage(-msg => "Set --no-polyA & --no-polyA-subset at the same time!", -exitval => 2, -verbose => 2) if ($no_polyA == 1 && $subsetFile ne '');
  46 pod2usage(-msg => "--transcript-to-gene-map and --allele-to-gene-map are mutually exclusive!", -exitval => 2, -verbose => 2) if (($mappingF ne "") && ($alleleMappingF ne ""));
  47 pod2usage(-msg => "--gtf and --allele-to-gene-map are mutually exclusive!", -exitval => 2, -verbose => 2) if (($gtfF ne "") && ($alleleMappingF ne ""));
  48 pod2usage(-msg => "Invalid number of arguments!", -exitval => 2, -verbose => 2) if (scalar(@ARGV) != 2);
  49
  50 if (!$bowtie && ($bowtie_path ne "")) { print "Warning: If Bowtie is not used, no need to set --bowtie-path option!\n"; }
  51 if (!$bowtie2 && ($bowtie2_path ne "")) { print "Warning: If Bowtie 2 is not used, no need to set --bowtie2-path option!\n"; }
  52
  53 my $type;
  54
  55 if ($gtfF ne "") { $type = 0; }
  56 else { $type = 1; }
  57
  58 my @list = split(/,/, $ARGV[0]);
  59 my $size = scalar(@list);
  60
  61 if ($size == 1 && (-d $list[0])) {
  62     my $dir = $list[0];
  63     @list = (<$dir/*.fa>, <$dir/*.fasta>);
  64     $size = scalar(@list);
  65 }
  66
  67 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);
  68
  69 if ($no_polyA) { $polyAChoice = 1 }
  70 elsif ($subsetFile ne "") { $polyAChoice = 2; }
  71
  72 if ($bowtie_path ne "") { $bowtie_path .= "/"; }
  73 if ($bowtie2_path ne "") { $bowtie2_path .= "/"; }
  74
  75 my $command = "";
  76
  77 if ($type == 0) {
  78     $"=" ";
  79     $command = "rsem-extract-reference-transcripts $ARGV[1] $quiet $gtfF";
  80     if ($mappingF ne "") { $command .= " 1 $mappingF"; }
  81     else { $command .= " 0"; }
  82     $command .= " @list";
  83     &runCommand($command);
  84 }
  85 else {
  86     $"=" ";
  87     $command = "rsem-synthesis-reference-transcripts $ARGV[1] $quiet";
  88     if ($mappingF ne "") { $command .= " 1 $mappingF"; }
  89     elsif ($alleleMappingF ne "") { $command .= " 2 $alleleMappingF"; }
  90     else { $command .= " 0"; }
  91     $command .= " @list";
  92     &runCommand($command);
  93 }
  94
  95 $command = "rsem-preref $ARGV[1].transcripts.fa $polyAChoice $ARGV[1] -l $polyALen";
  96 if ($polyAChoice == 2) { $command .= " -f $subsetFile"; }
  97 if ($quiet) { $command .= " -q"; }
  98
  99 &runCommand($command);
 100
 101 if ($bowtie) {
 102     $command = $bowtie_path."bowtie-build -f";
 103     if ($quiet) { $command .= " -q"; }
 104     $command .= " $ARGV[1].n2g.idx.fa $ARGV[1]";
 105
 106     &runCommand($command);
 107 }
 108
 109 if ($bowtie2) {
 110     $command = $bowtie2_path."bowtie2-build -f";
 111     if ($quiet) { $command .= " -q"; }
 112     $command .= " $ARGV[1].idx.fa $ARGV[1]";
 113
 114     &runCommand($command);
 115 }
 116
 117 __END__
 118
 119 =head1 NAME
 120
 121 rsem-prepare-reference
 122
 123 =head1 SYNOPSIS
 124
 125 rsem-prepare-reference [options] reference_fasta_file(s) reference_name
 126
 127 =head1 ARGUMENTS
 128
 129 =over
 130
 131 =item B<reference_fasta_file(s)>
 132
 133 Either a comma-separated list of Multi-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.
 134
 135 =item B<reference name>
 136
 137 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).
 138
 139 =back
 140
 141 =head1 OPTIONS
 142
 143 =over
 144
 145 =item B<--gtf> <file>
 146
 147 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.
 148
 149 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 Multi-FASTA files is its transcript_id.
 150
 151 (Default: off)
 152
 153 =item B<--transcript-to-gene-map> <file>
 154
 155 Use information from <file> to map from transcript (isoform) ids to gene ids.
 156 Each line of <file> should be of the form:
 157
 158 gene_id transcript_id
 159
 160 with the two fields separated by a tab character.
 161
 162 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.
 163
 164 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.
 165
 166 (Default: off)
 167
 168 =item B<--allele-to-gene-map> <file>
 169
 170 Use information from <file> to provide gene_id and transcript_id information for each allele-specific transcript.
 171 Each line of <file> should be of the form:
 172
 173 gene_id transcript_id allele_id
 174
 175 with the fields separated by a tab character.
 176
 177 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 Multi-FASTA-formatted files.
 178
 179 (Default: off)
 180
 181 =item B<--no-polyA>
 182
 183 Do not add poly(A) tails to the end of reference isoforms. (Default: adding poly(A) tails to all transcripts)
 184
 185 =item B<--no-polyA-subset> <file>
 186
 187 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'.)
 188
 189 =item B<--polyA-length> <int>
 190
 191 The length of the poly(A) tails to be added. (Default: 125)
 192
 193 =item B<--bowtie>
 194
 195 Build Bowtie indices. (Default: off)
 196
 197 =item B<--bowtie-path> <path>
 198
 199 The path to the Bowtie executables. (Default: the path to Bowtie executables is assumed to be in the user's PATH environment variable)
 200
 201 =item B<--bowtie2>
 202
 203 Build Bowtie 2 indices. (Default: off)
 204
 205 =item B<--bowtie2-path>
 206
 207 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)
 208
 209 =item B<-q/--quiet>
 210
 211 Suppress the output of logging information. (Default: off)
 212
 213 =item B<-h/--help>
 214
 215 Show help information.
 216
 217 =back
 218
 219 =head1 DESCRIPTION
 220
 221 This program extracts/preprocesses the reference sequences for RSEM. It can optionally build Bowtie indices (with '--bowtie' option) and/or Bowtie 2 indices (with '--bowtie2' option) using their default parameters. If an alternative aligner is to be used, indices for that particular aligner can be built from either 'reference_name.idx.fa' or 'reference_name.n2g.idx.fa' (see OUTPUT for details). This program is used in conjunction with the 'rsem-calculate-expression' program.
 222
 223 =head1 OUTPUT
 224
 225 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', 'reference_name.n2g.idx.fa' and optional Bowtie/Bowtie 2 index files.
 226
 227 'reference_name.grp', 'reference_name.ti', 'reference_name.seq', and 'reference_name.chrlist' are used by RSEM internally.
 228
 229 B<'reference_name.transcripts.fa'> contains the extracted reference transcripts in Multi-FASTA format. Poly(A) tails are not added and it may contain lower case bases in its sequences if the corresponding genomic regions are soft-masked.
 230
 231 B<'reference_name.idx.fa' and 'reference_name.n2g.idx.fa'> are used by aligners to build their own indices. Poly(A) tails are added unless '--no-polyA' is set. In addition, all sequence bases are converted into upper case. The only difference between 'reference_name.idx.fa' and 'reference_name.n2g.idx.fa' is that 'reference_name.n2g.idx.fa' in addition converts all 'N' characters to 'G' characters. This conversion is in particular desired for aligners (e.g. Bowtie) that do not allow reads to overlap with 'N' characters in the reference sequences. Otherwise, 'reference_name.idx.fa' should be used to build the aligner's index files. RSEM uses 'reference_name.idx.fa' to build Bowtie 2 indices and 'reference_name.n2g.idx.fa' to build Bowtie indices. For visualizing the transcript-coordinate-based BAM files generated by RSEM in IGV, 'reference_name.idx.fa' should be imported as a "genome" (see Visualization section in README.md for details).
 232
 233 =head1 EXAMPLES
 234
 235 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'.
 236
 237 There are two ways to write the command:
 238
 239  rsem-prepare-reference --gtf mm9.gtf \
 240                         --transcript-to-gene-map knownIsoforms.txt \
 241                         --bowtie \
 242                         --bowtie-path /sw/bowtie \
 243                         /data/mm9/chr1.fa,/data/mm9/chr2.fa,...,/data/mm9/chrM.fa \
 244                         /ref/mouse_125
 245
 246 OR
 247
 248  rsem-prepare-reference --gtf mm9.gtf \
 249                         --transcript-to-gene-map knownIsoforms.txt \
 250                         --bowtie \
 251                         --bowtie-path /sw/bowtie \
 252                         /data/mm9 \
 253                         /ref/mouse_125
 254
 255 2) Suppose we also want to build Bowtie 2 indices in the above example and Bowtie 2 executables are found in '/sw/bowtie2', the command will be:
 256
 257  rsem-prepare-reference --gtf mm9.gtf \
 258                         --transcript-to-gene-map knownIsoforms.txt \
 259                         --bowtie \
 260                         --bowtie-path /sw/bowtie \
 261                         --bowtie2 \
 262                         --bowtie2-path /sw/bowtie2 \
 263                         /data/mm9 \
 264                         /ref/mouse_125
 265
 266 3) 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/Bowtie 2 indices, and will use an alternative aligner to align reads against either 'mouse_0.idx.fa' or 'mouse_0.idx.n2g.fa':
 267
 268  rsem-prepare-reference --transcript-to-gene-map mapping.txt \
 269                         --no-polyA \
 270                         mm9.fasta \
 271                         mouse_0
 272
 273 =cut