Question: liftover from Crab Eating Macau to hg_19
0
gravatar for nir.galun
6 months ago by
nir.galun0
nir.galun0 wrote:

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

alignment R • 336 views
ADD COMMENTlink modified 6 months ago by Michael Dondrup45k • written 6 months ago by nir.galun0

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

ADD REPLYlink written 6 months ago by h.mon24k

Yes, but it didn't seem to work.

ADD REPLYlink written 6 months ago by nir.galun0

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

ADD REPLYlink written 6 months ago by h.mon24k

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.

ADD REPLYlink modified 6 months ago by WouterDeCoster37k • written 6 months ago by nir.galun0

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:

101010 Button

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

ADD REPLYlink written 6 months ago by WouterDeCoster37k

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.

ADD REPLYlink written 6 months ago by Michael Dondrup45k

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.

ADD REPLYlink modified 6 months ago • written 6 months ago by nir.galun0
1

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.

ADD REPLYlink modified 6 months ago • written 6 months ago by Michael Dondrup45k
0
gravatar for Michael Dondrup
6 months ago by
Bergen, Norway
Michael Dondrup45k wrote:

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.

ADD COMMENTlink modified 6 months ago • written 6 months ago by Michael Dondrup45k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1374 users visited in the last hour