當(dāng)我們拿到某個研究的
GWAS summary data時头遭，因為大多是二代測序的緣故寓免，它的SNP ID可能不是正常的rsid计维，但是，我們下游有些分析卻要用到rsid蜈首，這時我們會想到去dbSNP庫去下載對應(yīng)版本參考基因組的rsID欠母。
當(dāng)下載完畢，并且整理成只剩下SNP CHR POS這三列時赏淌，又發(fā)現(xiàn)文件太大了，如果我們沒有基因型數(shù)據(jù)支持我們利用PLINK軟件注釋SNP ID的話六水，就會很難辦，因為用R或Python讀取大型文件又是個問題（
GRCh37大概15G）睛榄，這時候胯盯，當(dāng)然是分染色體去進(jìn)行注釋咯。
第一步
使用awk對SNP CHR POS三列文件的GRCh37 參考SNP ID分染色體保存
awk '{print > $1".txt"}' dbSNP_GRCh37
mkdir GRCh37 && mv *.txt GRCh37
第二步
要注釋的文件需包括CHR, POS, SNP三列
library(data.table)
setwd("~/new_run")
file = "reformatMETAL.gz"
df = fread(file)[, c(CHR, POS, SNP, A1, A2, freq, beta, SE, p, N)]
out = data.table()
for (chr in c(1:22, "X")){
setwd("~/dbSNP/GRCh37")
db = fread(paste0(chr, ".txt"), header=FALSE, col.names=c("CHR", "POS", "SNP"))
df_merge = merge(df, db, by=c("CHR", "POS"), all.x=TRUE, suffix=c("df", "db"))
df_merge$SNP = ifelse(is.na(df_merge$SNP.db), df_merge$SNP.df, df_merge$SNP.db)
out = rbind(out, df_merge)
}
fwrite(out, file="my_test.gz", row.names=FALSE, quote=FALSE, sep="\t")
就是這么簡單憎乙，散會２嫒ぁＥ⒈摺疗杉！
