【一文就够】批量查询SNP的基因组位置、附近序列等信息（hg38/grch38）

2年前，我连续发了2个帖子，分别解决，

1、在“只知道SNP号（rsID）时，如何批量查询SNP的基因组位置信息”的问题（【一文就够】批量查询SNP的基因组位置等信息）

2、在“只知道SNP的基因组位置，如何查询SNP附近序列”的问题（【一文就够】批量下载SNP附近序列）

正好有坛友问我，如何把这2步合并成1步，即“只知道SNP号（rsID），如何查询SNP附近序列”。

在工作之余，我就顺手完成了。考虑到SNP注释信息的更新，我也一并进行了更新。

这些更新包括：

1、SNP的位置信息，从GRCh37/hg19，更新到GRCh38/hg38。

2、dbSNP数据库的更新，新增了许多SNP。

3、相应地，有些二等位SNP，更新为了多等位SNP。

以下是具体代码（已经过测试）。示例的输入文件“snp.csv”，见附件。格式如下图，

```r

# BiocManager::install("biomaRt", version = "3.11")

# BiocManager::install("BSgenome", version = "3.11")

# BiocManager::install("BSgenome.Hsapiens.NCBI.GRCh38", version = "3.11")

library(biomaRt)

library(BSgenome)

library(BSgenome.Hsapiens.NCBI.GRCh38)

library(dplyr)

df <- read.csv("snp.csv", header=T,stringsAsFactors = F)

# 选择b38版本基因组，SNP数据

mart.snp <- useMart(host=" http://aug2017.archive.ensembl.org ", # Ensembl90

          biomart="ENSEMBL_MART_SNP",

          dataset="hsapiens_snp")

getAnno <- function(rs = "rs3043732", mart = mart.snp) {

 results <- getBM(attributes = c( "chr_name","refsnp_id","chrom_start","chrom_end","allele"),

          filters  = "snp_filter", values = rs, mart = mart)

 return(results)

}

want <- df[,1]

snp <- getAnno(want) %>% 

 select(rsID = refsnp_id,

     chr = chr_name,

     pos = chrom_start,

     allele = allele) %>% 

 group_by(rsID, pos) %>%

 filter(row_number() == 1) %>%

 ungroup() %>%

 filter(!grepl('HSC|HG', chr))

# save(snp, file="snp.rda")

# load("snp.rda")

df1 <- snp

# Four cols with head:

# rsID chr pos allele

# rs778661560 chr12 33538213 [A/C]

offset <- 100

seq <- as.data.frame(paste0(getSeq(Hsapiens,df1$chr,start=(df1$pos-offset),end=(df1$pos-1)),

              "[",

              df1$allele,

              "]",

              getSeq(Hsapiens,df1$chr,start=(df1$pos+1),end=(df1$pos+offset))))

seq$rsID <- df1$rsID

seq$allele <- df1$allele

seq$chr <- df1$chr

seq$pos <- df1$pos

colnames(seq) <- c("seq","rsID","allele","chr","pos")

seq <- seq[,c("rsID","chr","pos","allele","seq")]

write.csv(seq,"snp-nearby100bp-Ensembl90.csv",

     row.names=F)

```

运行结果如下图，

保存的csv文件，使用excel打开，即可查看。

enjoy~

~~~~~~~~~~~~~~~~~~~~~~~~~~~~

另外，我分享了一个目录：版内《生信问题及解决方案》目录汇总

旨在 完善知识体系、方便查询所需代码、减少搜索的时间；

~~~喜欢我的分享，请用丁当鼓励我吧！~~~~

如你想我发布其他的生信经验分享，可以私信留言给我，

我会不定期挑选一些发布在坛子里

~~~~~~~~~~~~~~~~~~~~~~~~~~~~