一行命令批量修改染色体和位置为RS号

鉴于太多人问我怎么批量根据chr:pos查找RS号,在这里出一个教程。

注意以下教程展示的是修改hg19基因组版本的RS号,如果你的数据是其他版本的,请修改为对应版本的数据。

假定数据是test.txt,内容如下:

现在希望根据第一列chr:pos找到对应的RS号,实现以下的效果:

则可以用dplyr::left_join参数,具体实现过程如下所示:

install.packages(dplyr)
library(dplyr)
tes = read.table("test.txt",header=T,check.names=F,sep="\t") #注意这里我设置的是制表符分隔符,如果你的文件不是制表符的话,需要修改成对应的分隔符
match = read.table("snp150_hg19.txt",header=T,check.names=F,sep="\t") 
need=dplyr::left_join(tes,match,by="chromosome:start") #如果snp150_hg19.txt文件中有对应的RS号,则比对到test.txt文件中,如果没有的话,就变为NA
write.table(need, file ="clean.txt", sep ="\t", row.names =FALSE, col.names =TRUE, quote =FALSE) #保存文件

上述命令需要用到snp150_hg19.txt文件,其内容如下所示:

snp150_hg19.txt文件可以从https://hgdownload.soe.ucsc.edu/downloads.html网站 下载,也可以通过后台回复RS号修改获取,下载后请解压后再使用;

这个文件很大,总共有234104111个位点,解压后5G左右,只能通过网盘传输;

另外,我也提供了hg38的基因组版本snp151_hg38.txt.gz,需要的同样后台回复RS号修改

由于本人只做人类,所以只提供了人类基因组版本的教程,做其他物种的小伙伴们需要自行从https://hgdownload.soe.ucsc.edu/downloads.html网站 下载。


致谢橙子牛奶糖(陈文燕),请用参考模版:We thank the blogger (orange_milk_sugar, Wenyan Chen) for XXX

感谢小可爱们多年来的陪伴, 我与你们一起成长~

原文地址:https://www.cnblogs.com/chenwenyan/p/15626207.html