Extracting multiple regions from MAF file

[Urjaswita Yadav]

Dear UCSC support team,

I have ChiP-Seq peak regions in rat and I am interested in extracting the corresponding multiple genome alignment blocks from UCSC MAF files. Could you suggest an easy tool to do that? Basically I want to take a bed file with rat coordinates as input and extract all the corresponding regions wherever there is one with all other species in the 99 vertebrate alignment MAF files.

I tried mafExtractor from mafTools, which works fine for a single start-end coordinate, but there are a few complications. 

1. It does not take multiple maf files as input (using *.maf), so I need to know which file will have my sequence. Or maybe there is another way?

2. It takes only one set of start-end pair as input, but I have ~2000.

3. In the peak file output by Macs, I think the coordinates always correspond to + strand, but in maf files regions may be reverse complemented also. So I need to convert the coordinates to + strand and then input in mafExtractor, but I don't know which of my peaks may correspond to - strand alignments.

I am trying to program this, but it seems like a basic problem and I wonder if there is an easy tool already available? Thanks a lot for your help!

[Brian Lee]

Dear Urjaswita Yadav,

Thank you for using the UCSC Genome Browser and your question about extracting multiple regions from MAF files.

You may be interested in the utility called mafFetch. You can find a precompiled version here: http://hgdownload.soe.ucsc.edu/admin/exe/

Here is an example of the usage statement:

$ mafFetch
mafFetch - get overlapping records from an MAF using an index table
usage:
mafFetch db table overBed mafOut

Select MAF records overlapping records in the BED using the
the database table to lookup the file and record offset.
Only the first 3 columns are required in the bed.

For example, if you had a region around the gene SOD1 on chr21 of interest you could use an input file called overBed like the following:

chr21 33031935 33041243

Then you could use the command as follows:

mafFetch hg19 multiz100way overBed mafOut

(see the "Using the MySQL server with our utilities" section about enabling the tools with installed MySQL to access our public server: http://genome.ucsc.edu/goldenPath/help/mysql.html)

The results would be a file called mafOut.

You may be interested to learn you could further turn these results into a visualizable file called "bigMaf". Following the steps on our bigMaf help page ( http://genome.ucsc.edu/goldenPath/help/bigMaf.html ) you could create a file that when hosted at a server that supports byte-ranges can be visualized.

mafToBigMaf hg19 mafOut stdout | sort -k1,1 -k2,2n > bigMaf.txt

bedToBigBed -type=bed3+1 -as=bigMaf.as -tab bigMaf.txt hg19.chrom.sizes bigMaf.bb

Then use as a custom track:
track type=bigMaf name="bigMaf Example without summary or frames files (see help page)" bigDataUrl=http://URL/to/bigMaf.bb

Thank you again for your inquiry and using the UCSC Genome Browser. If you have any further questions, please reply to gen...@soe.ucsc.edu. All messages sent to that address are archived on a publicly-accessible forum. If your question includes sensitive data, you may send it instead to genom...@soe.ucsc.edu.

All the best,

Brian Lee
UCSC Genomics Institute

--

---
You received this message because you are subscribed to the Google Groups "UCSC Genome Browser Public Support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to genome+un...@soe.ucsc.edu.
To post to this group, send email to gen...@soe.ucsc.edu.
Visit this group at https://groups.google.com/a/soe.ucsc.edu/group/genome/.
To view this discussion on the web visit https://groups.google.com/a/soe.ucsc.edu/d/msgid/genome/CALwtSZKQQ5DcxmgVteFgTfdYyLmwfXkToot1qnCrRX_%2BjBuJyg%40mail.gmail.com.
For more options, visit https://groups.google.com/a/soe.ucsc.edu/d/optout.