lifting hg18 to hg38

Question

[hi.there]

Hi,

I am trying to align chromosome positions from genome reference hg18 to hg38 so that I can eventually update rsids. I have tried:

python liftmap.py -m plink.map -o ADNI_cluster_01_forward_757LONI

...where plink.map is the hg18 map file and ADNI_cluster_01_forward_757LONI is the prefix for the .bim, .fam and .bed file. The error is:

FileNotFoundError: [Errno 2] No such file or directory: 'ADNI_cluster_01_forward_757LONI.bed.unlifted'

I also have a chainfile for the conversion but am unsure of how to utilize it. Would anybody be able to help? It seems I am relegated to command line approaches as my version of r cannot import rtracklayer/bioconductor.

Thank you.

Answer 1

[LChart]

You can use the liftover tool from UCSC http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/ (whatever build you're running).

First map the .map file to .bed format, then lift it over, then push it back into a map format.

awk '{print $1"\t"$4-1"\t"$4"\t"$2"\t.\t."}' plink.map > plink_hg18.bed
liftOver plink_hg18.bed hg18ToHg38.over.chain.gz plink_hg38.bed plink_nomap
awk '{print $1"\t"$4"\t0\t"$3}' plink_hg38.bed > plink_hg38.map

[hi.there]

Thank you. Very appreciative.

[hi.there]

Hi LChart,

Very appreciative of your reply. I am having an issue though creating the other plink files. The ultimate goal is to update the rsids of the hg18 file as I have a file of rsids and chr position ends from UCSC in hg38. When I try to create the plink files with the new plink_hg38.bed and plink_hg38.map, I am getting an error with the command: system("./plink2 --file plink_hg38 --make-just-bim --out newBim")

The error is: Error: Empty .map file.

Would you know why this might be happening? Thanks.

[LChart]

Can you provide heads of these files?

[hi.there]

Hi LChart,

Thank you so much for the reply. Nothing reads out either in command line running "head plink_hg38.map" or in r "system('head plink_hg38.map')". Perhaps I am not understanding this process as well as I hoped I would. Do you have any suggestions?

[LChart]

If the heads are emtpy then something went wrong in either the extraction or liftover commands. Double-checking the contents of all intermediate files (and seeing when they become empty) will tell you where things got off track.

[hi.there]

Thanks LChart. My reply took a little extra time because of a hardware problem. plink_hg18.bed seems fine but it goes empty after the second line. There is nothing in plink_hg38.bed.

[LChart]

awk operates in a line-by-line fashion; so if plink_hg18.bed only has one line; this implies that the original .map file only has one line. Is this the case?

[hi.there]

Thank you LChart for your patience and help. I figured out where the problem was. The head of plink_hg18.bed looked as such: 1 696230 696231 rs12029736 . .

1 742428 742429 rs3094315 . .

1 743267 743268 rs3115860 . . ...thus prepending the string 'chr' in front of the first column fixed everything in the liftover. Thank you again for your help!