• 论坛首页
  • 我的丁香客
  • 找人
    查找好友
  • 更多
    丁香园
    丁香通
    丁香人才
    丁香会议
    丁香搜索
    丁香医生
    丁香无线
    丁香导航
    丁当铺
    文献求助
    医药数据库
    丁香诊所
    来问医生
登录 注册

生物信息

关注今日:15 | 主题:139794
论坛首页  >  生物信息学讨论版   >  Programming
  • 发帖
    每发1个新帖
    可以获得0.5个丁当奖励
  • 回帖

分享到:

  • 微信

    微信扫一扫

  • 微博
  • 丁香客
  • 复制网址

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

  • 只看楼主
  • 页码直达:
  • 直达末页
楼主 stringson
stringson
丁香园准中级站友

  • 148
    积分
  • 671
    得票
  • 460
    丁当
  • +1 积分
  • 1楼

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~

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

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

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

~~~喜欢我的分享,请用丁当鼓励我吧!~~~~

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

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

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



  • snp.csv(0.18k)
  • 邀请讨论
  • 不知道邀请谁?试试他们

    换一换
2020-09-03 11:36 浏览 : 2725 回复 : 0
  • 投票 4
  • 收藏 9
  • 打赏
  • 引用
  • 分享
    • 微信扫一扫

    • 新浪微博
    • 丁香客
    • 复制网址
  • 举报
    • 广告宣传推广
    • 政治敏感、违法虚假信息
    • 恶意灌水、重复发帖
    • 违规侵权、站友争执
    • 附件异常、链接失效
    • 其他
  • • 病史长达15年的老烂腿,多家医院诊疗无果,疼痛难忍……

关闭提示

需要2个丁当

丁香园旗下网站

  • 丁香园
  • 用药助手
  • 丁香通
  • 文献求助
  • 丁香人才
  • 丁香医生
  • 丁香导航
  • 丁香会议
  • 手机丁香园
  • 医药数据库

关于丁香园

  • 关于我们
  • 丁香园标志
  • 友情链接
  • 联系我们
  • 加盟丁香园
  • 版权声明
  • 资格证书

官方链接

  • 丁香志
  • 丁香园新浪微博
引用回复