1. Use "twoBitToFa -bed"
Download the latest twoBitToFa from Jim Kent's source, e.g. http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa (link for linux version)
Example:
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
usage:
twoBitToFa input.2bit output.fa
options:
-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:
/path/input.2bit:name
or
/path/input.2bit:name
or
/path/input.2bit:name:start-end
To use the program, you need 2bit file for your genome. To convert fasta file to 2bit, use faToTwoBit (http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit). It also requires your gene annotation in bed12 format.
2. Use "gffread -w" in Cufflinks utility (http://cufflinks.cbcb.umd.edu/gff.html)
Example:
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:
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
Options:
-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
Here is full usage:
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
Options:
-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
Hi Xianjun,
ReplyDeleteFor 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 }'
This comment has been removed by the author.
ReplyDeleteThanks for sharing these informations
ReplyDelete