Given a specific DNA polymer coming from hg19 (i.e. AATAATAATCA..), I need to localize it inside hg19 and retrieve all the annotations in the UCSC database.
I tried to search these polymers by BLAST or BLAT with no responses, but I'm sure that hg19 contains them since they are extracted from the FASTA sequence.
I could download the entire USCS mysql database, localize all the positions of the input sequence and query the database by start and end positions of annotations.
Is there a tool that can localize such short sequences, one by one, and retrieve all correlated annotations?
probably BLAT is not giving you an ouput because it is not a very unique sequence and exists all over the genome. What is the length of the sequence you are querying?
As you sure that sequence isn't repeat-masked ?
Here is a script that will output all the genomic locations of any DNA sequence you input.
It requires you to get a rather large fasta file for the hg19 genome
you can get it from here:
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
decompress it on your work folder and save the script then just save the script as a perl file script.pl
the command to find your sequence becomes:
perl script.pl
ex:
perl script.pl hg19.fa AATAATAATCA
I hope this helps
Duarte
########################### script.pl ##############################
use strict;
use Bio::SeqIO;
my $usage = "perl dnamotif.pl ";
my $fasta_filename = shift(@ARGV) or die("Usage: $usage $!");
my $motif = shift(@ARGV) or die("Usage: $usage $!");
my $motif_length = length ($motif);
my $fasta_parser = Bio::SeqIO->new(-file => $fasta_filename, -format => 'Fasta');
while(my $seq_obj = $fasta_parser->next_seq())
{
printf("Searching sequence '%s'...", $seq_obj->id);
if((my $pos = index($seq_obj->seq(), $motif)) != -1)
{
printf("motif found at position %d ... %d!\n", $pos + 1, $pos+1+$motif_length);
}
else
{
printf("motif not found.\n");
}
}
BLAT can be set to search with higher sensitivity at the expense of speed. If you are using the version of BLAT embedded in some Genome Browser website, you may not have access to the proper settings. There are several parameters that may affect its sensitivity and the ability to locate short sequences:
-tileSize=n
you want this to be less than the length of the sequence you are looking for (L), if it is more than half of the desired length, you may need to lower -stepSize=S to at most L-n
-ooc=N.ooc
I think you may want to avoid using this. It speeds up the search, but ignores overrepresented sequence patterns
-minMatch=mm
often this will be set to 2 for DNA for a better speed/sensitivity ratio. If n>L/2, however, this definitely needs to be set to 1, otherwise the hits will be missed.
Whether you use BLAT, or some other tool (grep, VMATCH, PRIMEX) you end up with a list of positions. One way to get the annotations around these positions is to use the Table Browser. Alternatively, R with biomaRt package can provide functions to extract annotations at certain positions.
Thank you so much.
Duarte: i'm just using the unmasked file you posted. Sequences' lengths range from 6 to 100 nucleotides, some of these occur only once in each chromosome, others may occur thousands of times, but i'm sure that every sequence is present in all chromosomes. I'm also working on new indexing data structures (such as suffix arrays or FM-index) better performed to extract these salients words and their locations but i need to understand if they have a biological relevance.
after you get the genomic positions of these sequences I agree with @Matej Lexa you can use the table browser to instersect these positions of features of interest.
I tend to just download all the annotations I am interested in bed format and then use tools like bedtools to intersect the files and find what these overlaps with those coordinates...
Good luck
- Badges
- Science topic
- Similar topics
- Philosophy
- Esthetics
- Interpretation