zoukankan      html  css  js  c++  java
  • 寻找与疾病相关的SNP位点——R语言从SNPedia批量提取搜索数据

       SNP是单核苷酸多态性,人的基因是相似的,有些位点上存在差异,这种某个位点的核苷酸差异就做单核苷酸多态性,它影响着生物的性状,影响着对某些疾病的易感性。SNPedia是一个SNP调査百科,它引用各种已经发布的文章,或者数据库信息对SNP位点进行描述,共享着人类基因组变异的信息。我们可以搜索某个SNP位点来寻找与之相关的信息,也可以根据相关疾病,症状来寻找相关的SNP。

    初次使用SNPedia

      SNPedia主页网址为http://snpedia.com/index.php/SNPedia,比如我想查找与crouzon综合症相关的SNP,只需要在SNPedia中搜索crouzon syndrome,即会出现许多相关的SNP搜索结果

      如果这时候我想看每个SNP的相关信息,我就要每个链接分别点进去

      后来发现我们只需要提取里面的部分信息,Orientation,Stabilized,Reference,Chromosome,Position,Gene,还有clinvar表格信息,这时候我们就可以从网页中利用RCurl包,XML包,正则表达是把所需要的内容提取出来,有效抓取有用信息。

    知识准备

    RCurl包和XML包

       在前一篇博文R语言从小木虫网页批量提取考研调剂信息 http://www.cnblogs.com/ywliao/p/6420501.html中已经提过,这里再提一个XML包中之前没有介绍的函数。
       readHTMLTable(doc) #doc 是XML或者HTML格式文本,可以是文件名,也可以是刚刚parse的html对象,该函数返回XML或HTML中的表格

    正则表达式

    这里阐述基本的正则表达式使用 **
      [ ]中括号,匹配中括号里面的任意字符,例如[a]匹配"a"
      [a-z]表示匹配a到z任意字母,[A-Z]匹配大写A到Z,[0-9]匹配0-9任意数字
      [ ]*中括号加*表示匹配任意次,[ ]+表示匹配至少一次,例如[a-zA-z,;: ]+表示匹配小写和大写字母,;:和空格至少一次
      [ a|b ] 匹配a或者b
      直接输入字符,实现精确定位。比如"apple[a-zA-z,;: ]+",定位到apple开头的后面匹配小写和大写字母,;:和空格至少一次的内容
      [u4E00-u9FA5]匹配汉字
    **

    R语言gregexpr函数

      使用方法:```gregexpr(pattern,istring, fixed = FALSE) #pattern就是要匹配正则表达是,istring是待匹配的字符串矢量,比如c("abc","cdf"),fixed,
    如果设置为true,默认pattern是真正的字符串,不会作为其它使用,相当于转义,
    函数返回列表,包括每个字符串的匹配长度和是否匹配)

    #实例
     这里直接上代码,代码里面有着详细解释,许多函数以后可以直接复制使用,或者放进一个自己做的R包
    

    !/usr/bin/env Rscript

    download <- function(strURL){

    输入网址返回html树格式文件

    strURL:网页链接地址 return: html树文件

    h <- basicTextGatherer()# 查看服务器返回的头信息
    txt <- getURL(strURL, headerfunction = h$update,.encoding="gbk") ## 字符串形式
    htmlParse(txt,asText=T,encoding="gbk") #选择gbk进行网页的解析
    }

    getinf <- function(strURL){

    主要提取网页信息函数

    strURL:网页链接网址 return:包括所要的所有信息的data.frame

    doc<- download(strURL)

    写如标题

    info<- data.frame("Title"=strsplit(xmlValue(getNodeSet(doc,'//title')[[1]])," -")[[1]][1]) #"rs... - SNPedia"进行split

    写入"Geno Mag Summary "table

    GMS_table <- readHTMLTable(doc)
    GMS_index <- 0
    if (length(GMS_table)>2){
    for (p in 1:6){
    if (length(GMS_table[[p]])3){
    GMS_index <- p
    }
    }
    }
    if (GMS_index!=0){
    for (i in 1:length(GMS_table[[GMS_index]])){
    tmp <- ""
    for (t in 1:nrow(GMS_table[[GMS_index]][i])){
    if(tmp
    ""){tmp <-as.vector(GMS_table[[GMS_index]][i][t,1])}else{
    tmp <- paste(tmp,as.vector(GMS_table[[GMS_index]][i][t,1]),sep=";")
    }
    }
    if (i1){info$Geno <-tmp}
    else if (i
    2){info$Mag <-tmp}
    else if (i==3){info$Summary <- tmp
    tmp <- ""
    }
    }
    }else{
    info$Geno <-" "
    info$Mag <-" "
    info$Summary <- " "
    }

    写入剩下table信息

    mes <- getNodeSet(doc,'//td')
    mes2 <- list()
    for (c in mes){
    d <- xmlValue(c)
    if (d""){
    }else{
    mes2=c(mes2,d)
    }
    }
    tmp <- greg_return_string("Make[-A-Za-z0-9_.%;(), ]+",mes2)
    if (length(tmp)
    2){info$"Make"=paste(strsplit(tmp[[1]]," ")[[1]][2],strsplit(tmp[[2]]," ")[[1]][2],sep=";")}else{info$"Make"=" "}
    for (i in (1:length(pattlistMainTable))){
    tmp <- greg_return_index(pattlistMainTable[[i]],mes2)
    if (i1 && length(tmp)1){info$"Orientation"=strsplit(mes2[[tmp+1]]," ")[[1]]}else if (i1 && length(tmp)!=1){info$"Orientation"=" "}
    else if (i
    2 && length(tmp)1){info$"Stabilized"=strsplit(mes2[[tmp+1]]," ")[[1]]}else if (i2 && length(tmp)!=1){info$"Stabilized"=" "}
    else if (i3 && length(tmp)1){info$"Reference"=strsplit(mes2[[tmp+1]]," ")[[1]]}else if (i3 && length(tmp)!=1){info$"Reference"=" "}
    else if (i
    4 && length(tmp)1){info$"Chromosome"=strsplit(mes2[[tmp+1]]," ")[[1]]}else if (i4 && length(tmp)!=1){info$"Chromosome"=" "}
    else if (i5 &&length(tmp)1){info$"Position"=strsplit(mes2[[tmp+1]]," ")[[1]]}else if (i5 && length(tmp)!=1){info$"Position"=" "}
    else if (i
    6&&length(tmp)1){info$"Gene"=strsplit(mes2[[tmp+1]]," ")[[1]]}else if (i6 && length(tmp)!=1){info$"Gene"=" "}
    }

    写入clivar

    mes <- getNodeSet(doc,'//tr')
    mes2 <- list()
    for (c in mes){
    d <- xmlValue(c)
    if (d""){
    }else{
    mes2=c(mes2,d)
    }
    }
    for (i in (1:length(pattlistClinvar))){
    tmp <- greg_return_string(pattlistClinvar[i],mes2)
    if (length(tmp)!=0){tmp <- tmp[[1]]}
    if (i
    1 && length(tmp)!=0){info$"Risk"=strsplit(tmp," ")[[1]][3]}else if (i1 && length(tmp)0){info$"Risk"=" "}
    else if (i2 && length(tmp)!=0){info$"Alt"=strsplit(tmp," ")[[1]][3]}else if (i2 && length(tmp)0){info$"Alt"=" "}
    else if (i
    3 && length(tmp)!=0){info$"ReferenceBase"=strsplit(tmp," ")[[1]][3]}else if (i3&& length(tmp)0){info$"ReferenceBase"=" "}
    else if (i4 && length(tmp)!=0){info$"Significance"=strsplit(tmp," ")[[1]][2]}else if (i4 && length(tmp)0){info$"Significance"=" "}
    else if (i
    5&& length(tmp)!=0){info$"Disease "=strsplit(tmp," ")[[1]][3]}else if (i5 && length(tmp)0){info$"Disease "=" "}
    else if (i6 && length(tmp)!=0){info$"CLNDBN"=strsplit(tmp," ")[[1]][3]}else if (i6 && length(tmp)0){info$"CLNDBN"=" "}
    else if (i
    7 && length(tmp)!=0){info$"Reversed"=strsplit(tmp," ")[[1]][3]}else if (i7 && length(tmp)0){info$"Reversed"=" "}
    else if (i8 && length(tmp)!=0){info$"HGVS"=strsplit(tmp," ")[[1]][3]}else if (i8 && length(tmp)0){info$"HGVS"=" "}
    else if (i
    9 && length(tmp)!=0){info$"CLNSRC"=strsplit(tmp," ")[[1]][3]}else if (i9 && length(tmp)0){info$"CLNSRC"=" "}
    else if (i10 && length(tmp)!=0){info$"CLNACC "=strsplit(tmp," ")[[1]][3]}else if (i10 && length(tmp)==0){info$"CLNACC "=" "}

    }

    info
    }

    greg_return_string <- function(pattern,stringlist){

    greg_return_stirng 指定匹配全部字符串列表,返回匹配的字符串

    pattern:匹配模式,比如"abc[a-z]*" stringlist:字符串列表,list("abc","abcde","cdfe") return : 列表里字符串匹配结果,"abc""abcde"

    findlist <- gregexpr(pattern,stringlist)
    needlist <- list()
    for (i in which(unlist(findlist)>0)){
    preadress <- substr(stringlist[i],findlist[[i]],findlist[[i]]+attr(findlist[[i]],'match.length')-1)
    needlist<- c(needlist,list(preadress))
    }
    return(needlist)
    }
    greg_return_index <- function(pattern,stringlist){

    greg_return_stirng 指定匹配全部字符串列表,返回存在匹配的字符串列表index

    pattern:匹配模式 stringlst:待匹配字符串列表 return:存在返回匹配的字符串在列表中的index

    findlist <- gregexpr(pattern,stringlist)
    needlist <- list()
    which(unlist(findlist)>0)
    }

    extradress <- function(strURL){

    将strURL网页里面我们所需要链接提取出来并加工

    strURL:网页链接网址 return:网址列表,包括所有提取加工后的网址链接

    pattern <- "/index.php/Rs[0-9]+"
    prefix <- "https://snpedia.com" #网址改为https起始
    links <- getHTMLLinks(download(strURL)) # getHTMLLinks不能解析https的网址,因此先用download解析网址
    needlinks <- gregexpr(pattern,links)
    needlinkslist <- list()
    for (i in which(unlist(needlinks)>0)){
    preadress <- substr(links[i],needlinks[[i]],needlinks[[i]]+attr(needlinks[[i]],'match.length')-1)
    needlinkslist<- c(needlinkslist,list(preadress))
    adresses <- lapply(needlinkslist,function(x)paste(prefix,x,sep=""))

    }
    adresses
    }

    greg <- function(pattern,istring){

    greg函数查看单个字符串istring,并且返回匹配的部分,不匹配返回空

    gregout <- gregexpr(pattern,istring)
    substr(istring,gregout[[1]],gregout[[1]]+attr(gregout[[1]],'match.length')-1)
    }

    library(RCurl)
    library(XML)

    自定义部分

    strURL <- "https://snpedia.com/index.php?title=Special%3ASearch&profile=default&fulltext=Search&search=Congenital+adrenal+hyperplasia" #snpedia网址都已改为https开头
    output <- "ouput.txt"

    message(paste("[prog]",strURL,output,sep=" "))

    strURLs <- extradress(strURL)
    pattlistMainTable <- list("Orientation$","Stabilized$","Reference$","Chromosome$","Position$","Gene$")

    此匹配模式列表用于返回该字符串所在index,而对应的值是index是该index+1

    pattlistClinvar <- list("Risk [-A-Za-z0-9_.%;(), ]+","Alt [-A-Za-z0-9_.%;(), ]+",
    "Reference [-A-Za-z0-9_.%;() ]+","Significance [A-Za-z ]+","Disease [A-Za-z ]+",
    "CLNDBN [-A-Za-z0-9_.% ]+","Reversed [0-9]+", "HGVS [-A-Za-z0-9_.%:> ]+","CLNSRC [-A-Za-z0-9_.% ]+","CLNACC [-A-Za-z0-9_.%, ]+")

    此匹配模式列表用于返回相应clinvar

    inf <- NULL
    for ( strURL in strURLs){
    dat <- getinf(strURL)
    if (length(inf) == 0){
    inf <- dat
    }else{
    inf <- rbind(inf,dat)
    }
    }

    write.table(inf, file = output, row.names = F, col.names=T,quote = F, sep=" ") # tab 分隔的文件
    message("完成!")

    
    结果可以直接打开,也可以用excel的自文本打开,方便查看
    ![](http://images2015.cnblogs.com/blog/1093203/201703/1093203-20170308182910844-1533827742.png)
  • 相关阅读:
    单链表的实现C语言版
    顺序表的基本方法实现C语言版
    算法
    Redis
    Linux安装python3.6.1
    Markdown 基本使用手册
    设计Django个人博客网站
    RabbitMQ消息队列
    堡垒机 Paramiko 模块
    进程、线程、协程总结
  • 原文地址:https://www.cnblogs.com/ywliao/p/6517063.html
Copyright © 2011-2022 走看看