liftover from Crab Eating Macau to hg_19

Question

[nir.galun]

Hi, how can I use rtracklayer and liftover function to liftover from Macaca Fascicularis genome to hg_19? This is an example of how my data frame looks like in R.

chr1 213 214 0.87738095

chr1 225 226 0.82777778

chr1 528 529 0.49857599

chr1 563 564 0.88650609

chr1 574 575 0.91061254

I already downloaded the chain file "macFas5ToHg19.over.chain" from UCSC. thanks.

here is the download link: ftp://hgdownload.cse.ucsc.edu/goldenPath/macFas5/liftOver/macFas5ToHg19.over.chain.gz

[h.mon]

Did you read Changing genomic coordinate systems with rtracklayer::liftOver? Did you try to follow its instructions?

[nir.galun]

Yes, but it didn't seem to work.

[h.mon]

Explain "it didn't seem to work". What command did you use? Did you get any error messages?

[nir.galun]

Hi,

This is the command I got stack In. 
> path = system.file(package="liftOver", "extdata", "macFas5ToHg19.over.chain.gz")
> ch = import.chain(path)
Error in connectionForResource(resource(x), open = open) : 
  path cannot be an empty string
> ch
Error: object 'ch' not found

Thanks.

[WouterDeCoster]

I added code markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:

Check the command where you assing to path. The issue is there.

[Michael]

Assuming that you downloaded the chain file, you have to specify the right path. The path is empty because it’s not found where system.file expects it.

[nir.galun]

Hi, after I ran this code:

library(gwascat)
cur = makeCurrentGwascat() 
data(cur)
cur
        library(rtracklayer)
        library(gwascat)
        library(IRanges)
        chain <- import.chain("C:/Users/nirgalun/Downloads/macFas5ToHg19.over.chain")
        str(chain[[1]])
        chain
        seqlevelsStyle(cur) = "UCSC" 
        cur19 = liftOver(cur,chain)
        class(cur19)
        cur19 = unlist(cur19)
        genome(cur19) = "hg19"
        cur19 = new("gwaswloc", cur19)
        cur19
        length(cur)-length(cur19)

    setdiff(mcols(cur)$SNPS, mcols(cur19)$SNPS)

I got:

1] "rs9275390"   "rs7584330"   "rs2242652"   "rs163182"    "rs1800562"   "rs855791"    "rs4783244"  
   [8] "rs3771180"   "rs2381416"   "rs4845783"   "rs3129943"   "rs9500927"   "rs2069408"   "rs10202497" 
  [15] "rs10745527"  "rs12344583"  "rs9349379"   "rs183211"    "rs6977820"   "rs10235789"  "rs3739070"  
  [22] "rs4900109"   "rs1024161"   "rs4947296"   "rs6457617"   "rs4950806"   "rs11675251"  "rs11842874

as a finale resolet. How can I convert that to genome positions? thanks.

[Michael]

It looks as if you are running code that you cannot fully understand. The last command gives you the difference between original snp set and hg19 snap set. Like the ones that couldn’t be converted.

Answer 1

[Michael]

how can I use rtracklayer and liftover function to liftover from Macaca Fascicularis genome to hg_19?

Even though this file exists, may I question to which extent it makes sense to lift-over between two different primate species? This certainly depends on your study and the scientific question, but it might mean that many loci will simply be dropped.

How can I convert them to genome positions?

First your example used human SNPs, therefore the outcome is not representative for the case where you use Macaque data. You already have them. Your lifted data is in cur19 after running your code.

class(cur19)

[1] "gwaswloc" attr(,"package") [1] "gwascat"

This tells you that your data now is in a object of Class gwaswloc which extends GRanges.

Therefore, see ?'gwaswloc-class' and ?GRanges for a list of accessor methods to retrieve genomic coordinates.

Some of the most basic examples:

start(cur19)
end(cur19)
seqnames(cur19)
data.frame(cur19)

The data.frame call is used to export the content into a TSV file.