This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Quick Search For Overlaps With A Mask File

Hi there,

I am interested in computing some statistics for a subset of the genome. I have compiled a list of 28M sites and I am now checking if this site passed the 1000 genomes filters. In order to do this, I am using samtools to access the mask file (in fasta format) for the particular position using -

samtools faidx maskfile chr_pos

However, for 28M sites, this is very slow. Does anyone have any suggestions to speed up the search? I am fairly good at Perl programming but I dont think using the standard search methods that I am familiar with will significantly speed things up. Any suggestions would be very helpful?

Also as step2, I need to lookup the results from a table for which I am currently using perl hashes but again too slow for such large amount of data. I would be most grateful if you have any suggestions to more the speed.

-Diviya

search
• 3,279 views • link
written by diviya.smith • diviya.smith

What are the attributes for these 28M sites ? How you are matching with genomes ? Are you trying to match based on some interval range etc?

• 271 views • link
written by always_learning always_learning

3 answers

you could always try using ANNOVAR filtering option, which would give you 1000 genomes' population frequencies, if by "checking if this site passed the 1000 genomes filters" you mean that you want to want to know if the 1000 genomes project detected those sites as variants:

annotate_variation.pl -filter -dbtype 1000g2012apr_eur -buildver hg19 YOUR_FILE humandb/
• 0 views • link
written by Jorge Amigo Jorge Amigo

Pardon my ignorance, But Can you please explain that how this is going help above posed question ? This is purely for my knowledge only !

• 0 views • link
written by always_learning always_learning

I think this is an attempt to be more specific about "this site passed the 1000 genomes filters".

• 0 views • link
written by Michael Dondrup Michael Dondrup

thanks Michael, you got it. I was just suggesting to filter (or to highlight) 1000 genomes' variants, if the original 28M sites are to be queried for known variation.

• 0 views • link
written by Jorge Amigo Jorge Amigo

You should calculate range overlaps between query regions and 1kg regions to approach this in the most efficient way. Therefore, convert both files to ranges either in GFF or BED format and use bedtools or Bioconductor/IRanges to calculate overlaps. See What is the quickest algorithm for range overlap?

To do this you need to face the following challenges:

  • Is the mask file the correct and optimal way of getting all the "site passed the 1000 genomes filters" and what does that mean?
  • How is the fact encoded in the mask file and how can this be converted to BED or GFF? If you show us where you got that file from would certainly help.
• 0 views • link
written by Michael Dondrup Michael Dondrup

Sorry, perhaps I was not very clear. I am interested in extracting all the CG sites across the genome (~28M) and checking if these sites passed the 1000 genomes filters (information for which is included in a mask file in fasta format). I am currently using samtools to pull out the sites from the fasta file but the search is very slow. Any suggestions on how I can improve the search time?

• 0 views • link
written by diviya.smith • diviya.smith

this should be a comment, not an answer. anyway, I would go for Michael's suggestion: getting 1000genomes' masked regions in bed format and intersecting them with your sites.

• 0 views • link
written by Jorge Amigo Jorge Amigo

Thanks all for the suggestions. This is very helpful.

• 0 views • link
written by diviya.smith • diviya.smith

Log in to answer this question.

More posts like this