This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Retrieve the sequence based on the start and end position in the cuffmerged.gtf

Hi guys,

Is there any one who knowns how to retrieve the gene sequence based on the staring and ending position in the cuffmerged.gtf file. Since there are some genes only tracking Ids and starting and ending positions available. I want to retrieve these sequences and annotate it. I will really appreciate for you guys help.

Thanks a lot

rna-seq

Thanks very much

Sounds helpful.

I will try it tomorrow.

2 answers

A one-line solution:

bedtools getfasta -fi genome.fa -bed cuffmerged.gtf -fo out.fa

Yes, the -bed parameter can actually take BED/GFF/VCF files. Full documentation here

Hi igor,

Thanks for your solution. However I want to extract the sequence corresponding to one cufflink tracking ID, the bedtools getfasta return several exon sequences for each tracking ID, Do you think there is a way to get around that? I also checked the cufflink website, there is a gffread utility.

which was designed to handle the cufflink output, however it extrac transcript sequences based on transcript ID in the cuffmerged.gtf not gene ID, Do you think there is a way to change it?

Thanks very much

First, you could build an indexed set of FASTA files for all of your chromosomes using samtools.

Second, you can use a conversion tool in the BEDOPS suite (for example) to generate a BED file containing coordinates for genes of interest, and run it against a BED-to-FASTA script to retrieve their sequences.

1)Say you are working withhg19. Create a work directory where your FASTA files and index files will go:

$ cd /foo/bar/baz
$ wget 'ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/*.fa.gz'
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

2)Use tools like BEDOPS gtf2bed to convert gene annotations to BED format.

For example:

$ gtf2bed < cuffmerged.gtf \
    | grep -w 'gene' - \
    | cut -f1-6 - \
    > genes.bed

We can convert genes to FASTA via a Perl script that callssamtoolsagainst the indexed reference sequence.

Here is an example of such a helper Perl script:

#!/usr/bin/env perl

use strict;
use warnings;

# fastaDir contains per-chromosome UCSC FASTA (.fa) files and samtools-indexed index (.fai) files
my $fastaDir = "/foo/bar/baz"; 

while (<STDIN>) {
    chomp;
    my ($chr, $start, $stop, $id, $score, $strand) = split("\t", $_);
    my $queryKey = "$chr:$start-$stop";
    my $queryResult = `samtools faidx $fastaDir/$chr.fa $queryKey`; chomp $queryResult;
    my @lines = split("\n", $queryResult);
    my @seqs = @lines[1..(scalar @lines - 1)];
    my $seq = join("", @seqs);
    if ($strand eq "-") { $seq = revdnacomp($seq); }
    my $header = ">".join(":",($chr, $start, $stop, $id, $score, $strand));
    print STDOUT $header."\n".uc($seq)."\n";
}

sub revdnacomp {
    my $dna = shift @_;
    my $revcomp = reverse($dna);
    $revcomp =~ tr/ACGTacgt/TGCAtgca/;
    return $revcomp;
}

Replace /foo/bar/baz with wherever you stored the FASTA sequence and index files.

You might run this script like so:

$ bed2fasta.pl < genes.bed > genes.fa

I like igor's solution, but wanted to mention that if you get a bed file you can pass it to pyfaidx without even having samtools installed:

$ pip install pyfaidx
$ faidx --bed genes.bed genome.fasta > genes.fasta

Clever alternative.

Log in to answer this question.