Friday, January 25, 2013

quick ways to extract mRNA/transcript sequence for a gene

NOTE: please be aware that gtf_to_fasta tool in Tophat does NOT consider the strand info in input gff. 

1. Use "twoBitToFa -bed"
Download the latest twoBitToFa from Jim Kent's source, e.g. (link for linux version)

twoBitToFa -bed=input.bed /path/input.2bit inputbed.cDNA.fa

Here is the full usage:

twoBitToFa - Convert all or part of .2bit file to fasta
   twoBitToFa input.2bit output.fa
   -seq=name - restrict this to just one sequence
   -start=X  - start at given position in sequence (zero-based)
   -end=X - end at given position in sequence (non-inclusive)
   -seqList=file - file containing list of the desired sequence names
                    in the format seqSpec[:start-end], e.g. chr1 or chr1:0-189
                    where coordinates are half-open zero-based, i.e. [start,end)
   -noMask - convert sequence to all upper case
   -bpt=index.bpt - use bpt index instead of built in one
   -bed=input.bed - grab sequences specified by input.bed. Will exclude introns

Sequence and range may also be specified as part of the input
file name using the syntax:

To use the program, you need 2bit file for your genome. To convert fasta file to 2bit, use faToTwoBit ( It also requires your gene annotation in bed12 format.

2. Use "gffread -w" in Cufflinks utility (

gffread -w transcripts.fa -g /path/to/genome.fa transcripts.gtf

where the input annotation can be gtf or gff2/3.

Note that this tool is powerful in only extracting the CDS/protein etc. And to get fast, it's better to have genome.fa.fai index file in the same path.

Here is full usage:

gffread <input_gff> [-g <genomic_seqs_fasta> | <dir>][-s <seq_info.fsize>]
 [-o <outfile.gff>] [-t <tname>] [-r [[<strand>]<chr>:]<start>..<end> [-R]]
 [-CTVNJMKQAFGUBHZWTOLE] [-w <exons.fa>] [-x <cds.fa>] [-y <tr_cds.fa>]
 [-i <maxintron>]
 Filters and/or converts GFF3/GTF2 records.
 <input_gff> is a GFF file, use '-' if the GFF records will be given at stdin

  -g  full path to a multi-fasta file with the genomic sequences
      for all input mappings, OR a directory with single-fasta files
      (one per genomic sequence, with file names matching sequence names)
  -s  <seq_info.fsize> is a tab-delimited file providing this info
      for each of the mapped sequences:
      <seq-name> <seq-length> <seq-description>
      (useful for -A option with mRNA/EST/protein mappings)
  -i  discard transcripts having an intron larger than <maxintron>
  -r  only show transcripts overlapping coordinate range <start>..<end>
      (on chromosome/contig <chr>, strand <strand> if provided)
  -R  for -r option, discard all transcripts that are not fully
      contained within the given range
  -U  discard single-exon transcripts
  -C  coding only: discard mRNAs that have no CDS feature
  -F  full GFF attribute preservation (all attributes are shown)
  -G  only parse additional exon attributes from the first exon
      and move them to the mRNA level (useful for GTF input)
  -A  use the description field from <seq_info.fsize> and add it
      as the value for a 'descr' attribute to the GFF record

  -O  process also non-transcript GFF records (by default non-transcript
      records are ignored)
  -V  discard any mRNAs with CDS having in-frame stop codons
  -H  for -V option, check and adjust the starting CDS phase
      if the original phase leads to a translation with an
      in-frame stop codon
  -B  for -V option, single-exon transcripts are also checked on the
      opposite strand
  -N  discard multi-exon mRNAs that have any intron with a non-canonical
      splice site consensus (i.e. not GT-AG, GC-AG or AT-AC)
  -J  discard any mRNAs that either lack initial START codon
      or the terminal STOP codon, or have an in-frame stop codon
      (only print mRNAs with a fulll, valid CDS)

  -M/--merge : cluster the input transcripts into loci, collapsing matching
       transcripts (those with the same exact introns and fully contained)
  -d <dupinfo> : for -M option, write collapsing info to file <dupinfo>
  --cluster-only: same as --merge but without collapsing matching transcripts
  -K  for -M option: also collapse shorter, fully contained transcripts
      with fewer introns than the container
  -Q  for -M option, remove the containment restriction:
      (multi-exon transcripts will be collapsed if just their introns match,
      while single-exon transcripts can partially overlap (80%))

  -E  expose (warn about) duplicate transcript IDs and other potential
      problems with the given GFF/GTF records
  -Z  merge close exons into a single exon (for intron size<4)
  -w  write a fasta file with spliced exons for each GFF transcript
  -x  write a fasta file with spliced CDS for each GFF transcript
  -W  for -w and -x options, also write for each fasta record the exon
      coordinates projected onto the spliced sequence
  -y  write a protein fasta file with the translation of CDS for each record
  -L  Ensembl GTF to GFF3 conversion (implies -F; should be used with -m)
  -m  <chr_replace> is a reference (genomic) sequence replacement table with
      this format:
      <original_ref_ID> <new_ref_ID>
      GFF records on reference sequences that are not found among the
      <original_ref_ID> entries in this file will be filtered out
  -o  the "filtered" GFF records will be written to <outfile.gff>
      (use -o- for printing to stdout)
  -t  use <trackname> in the second column of each GFF output line
  -T  -o option will output GTF format instead of GFF3


  1. Chirag6:02 AM

    Hi Xianjun,

    For eg: how would you do this for Bed12, if you want only mRNA from the bed12. If you could illustrate an example, it would be nice.

    This is how i do:

    bed12ToBed6 -i File.bed | fastaFromBed -fi genome.fa -bed stdin -s -name -fo stdout -tab | groupBy -g 1 -o concat -c 2 | awk '{ print ">"$1"\n"$2 }'

  2. This comment has been removed by the author.

  3. Thanks for sharing these informations