yes, this is brilliant
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
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
Log in to answer this question.
Thanks very much
Sounds helpful.
I will try it tomorrow.