之前我们爬取过NCBI基因信息，DailyMed药物信息，今天来爬取 PubChem中的药物，包括药物名、各种ID、2级结构图片等信息。
点击查看之前的消息：
使用R语言爬取DailyMed药物信息
R语言批量爬取NCBI基因注释数据
假如我现在有如下药物清单（大约5000个），想批量获取其在pubchem中的信息。
包括 分子式、分子量结构等信息：
pubchem中每个药物有很多别名，但药物的ID是唯一的，也就是说一个ID对应多个药物名，为了准确的获得我们药物清单中的药物信息，首先得获得该药物的ID，pubchem药物ID-drug name对应信息文件可以在 ftp://ftp.ncbi.nlm.nih.gov/pubchem/Compound/Extras/ 中获得:
下载 文件
CID-Synonym-filtered，这个文件只包含两列，第一列是CID，第二列是药物名。
现在将我们的药物清单与该文件的药物名进行匹配，shell代码及注释如下:
inputfile="drug.txt"
outfile="./pubchem_CID.txt"
tmpfile="./tmp.txt"
rm "$tmpfile"
:>"${outfile}"
i=0
awk 'BEGIN{FS="\t"}{print $1}' "$inputfile" | while read line
do
# 读入每行第一列的药名：
drug_name[${i}]=`echo ${line} `
# 跳过第一行 colname ，然后进行 grep 药物名。
if(($i>0));then
echo "正在处理药物："$i" "${drug_name[$i]}
j=0
# 将 这个药 所有匹配到的行拆分保存到两个数组：drug_CID，drug_name
# 如果完全匹配则跳过，否则选择最短那个 grep项
# 注意只在前20000万匹配项中找，否则有的有几十万的太久了。
grep -i "${drug_name[$i]}" ./Extra/CID-Synonym-filtered -m 10000 >"$tmpfile"
# 如果文件为空，则只输出药物名：
if [[ ! -s "$tmpfile" ]];then
echo -e "${drug_name[$i]}\tNA\tNA\t文件为空，没有匹配到"
echo -e "${drug_name[$i]}\tNA\tNA\t没有匹配到" >>"${outfile}"
fi
# 如果文件不为空：
if [[ -s "$tmpfile" ]];then
j=0
nrow_tmp=`wc -l "$tmpfile"|awk '{print $1}'`
# 注意这里 的IFS 重新定义分割符 和 -r选项保证读入的内容是原始的内容，意味着反斜杠转义的行为不会发生
OLD_IFS=$IFS #保存原始值
cat "$tmpfile" | while IFS=$'\t' read -r col1 col2
do
# 存储原始 CID 和 name:
drug_CID[$j]="${col1}"
drug_name[$j]="${col2}"
IFS=$OLD_IFS #还原IFS的原始值
# 药物名转为小写进行比较：
drug_tmp_name["$j"]=`echo "${col2}" | awk '{print tolower($0)}'`
A_tmp_name=`echo "${drug_name[$i]}" | awk '{print tolower($0)}'`
# 如果完全匹配，然后赋值：
if [ "${A_tmp_name}" == "${drug_tmp_name[$j]}" ];then
echo -e "${drug_name[$i]}\t${drug_CID[$j]}\t${drug_name[$j]}\t完全匹配到"
echo -e "${drug_name[$i]}\t${drug_CID[$j]}\t${drug_name[$j]}\t完全匹配" >>"${outfile}"
break 1
fi
# 如果没有完全匹配到的，循环结束后选择最短的一个 drug_name 项，因为最短的最相似：
let j=j+1
if [ "${nrow_tmp}" -eq "${j}" ] ;then
OLD_IFS=$IFS #保存原始值
IFS=$'\t' #改变IFS的值
# 找出最小的那个数组元素的序号：
min_order=`awk -v drug_name="${drug_tmp_name[*]}" 'BEGIN{ split(drug_name,a,"\t"); OFS="\t" }END{min=length(a[1]);min_order=1;
for(i=1;i<=length(a);i++){
if( length(a[i]) <= min ){min=length(a[i]);min_order=i}
#print i,a[i],length(a[i])
}
print min_order
} '`
IFS=$OLD_IFS #还原IFS的原始值
let min_order_shell=$min_order-1
echo -e "${drug_name[$i]}\t${drug_CID[$min_order_shell]}\t${drug_name[$min_order_shell]}\t没有完全匹配，选第"$min_order"个"
echo -e "${drug_name[$i]}\t${drug_CID[$min_order_shell]}\t${drug_name[$min_order_shell]}\t选第"$min_order"个" >>"${outfile}"
fi
done
# 删除数组：
unset drug_CID
unset drug_name
unset drug_tmp_name
unset A_tmp_name
rm "$tmpfile"
fi
fi
let i=i+1
done
结果分为三类，第一类是药物名完全匹配的，第二类是pubchem中没有的，第三类是与pubchem药物名部分匹配的，对于第三类结果，选择相似度最高的药物名最为最终结果，部分结果如下：
然后根据CID爬取每个药物的信息，R代码及注释如下：
rm(list=ls())
library(data.table)
library(openxlsx)
library(stringr)
library(stringi)
library(VennDiagram)
library(RCurl)
library(stringr)
library(XML)
library(downloader)
# 读取包含 CID 的并且手动检查修改过的药物信息文件：
A_drug <- fread("./pubchem_id.txt",header = T,stringsAsFactors = F)
A_drug <- A_drug [,c("A_name","pubchem_CID","pubchem_name")]
A_drug $pubchem_url <- paste("https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/",A_drug $pubchem_CID,"/XML/?response_type=save&response_basename=CID_",A_drug $pubchem_CID,sep="")
A_drug $img_2D_url <- paste("https://pubchem.ncbi.nlm.nih.gov/image/imagefly.cgi?cid=",A_drug $pubchem_CID,"&width=500&height=500",sep="")
for(i in 1:nrow(A_drug )){
# 如果药物CID为NA则跳过循环：
if(is.na(A_drug [i,pubchem_CID])){
print(paste("处理第",i,"个药,CID为：",A_drug [i,pubchem_CID],",已经跳过"))
next
}
# 定义XML文件路径，和图片路径：
drug_XML_file <- "./pubchem_tmp.xml"
Img_2D_file <- paste("./Img_2D/",A_drug [i,pubchem_CID],".png",sep="")
cat(paste("处理第",i,"个药,CID为：",A_drug [i,pubchem_CID],"\t"))
# XML文件存在则直接删除：
if(file.exists(drug_XML_file)){
file.remove(drug_XML_file)
cat("XML文件存在，已经删除!")
# 下载每个药物的 XML 文件：
download.file(A_drug [i,pubchem_url], destfile = drug_XML_file, quiet = TRUE)
cat(ifelse(file.exists(drug_XML_file),"下载成功!","下载失败!"))
}else if(!file.exists(drug_XML_file)){
cat("XML文件不存在")
# 下载每个药物的 XML 文件：
download.file(A_drug [i,pubchem_url], destfile = drug_XML_file, quiet = TRUE)
cat(ifelse(file.exists(drug_XML_file),"下载成功!","下载失败!"))
}
# 如果文件存在则进行下面任务：
if(file.exists(drug_XML_file)){
# 由于在根节点包含 namespace 信息，而有 namespace 无法使用 Xpath提取信息，所以去掉 namespace 信息：
edit_file <- readLines(drug_XML_file, encoding = "UTF-8")
# 由于部分字符编码不一致，去除这些不一致字符如：℃ 等
enc<-lapply(edit_file,function(x){stri_enc_mark(x)})
enc<-as.character(enc)
if(length(which(enc=="UTF-8"))>0){
remove_Non_ASCII<-function(x){
# 找出哪个不是ASCII编码的,注意这里用 str_split 会出错，特殊符号偶尔分不开来。
x_split <- strsplit(x,"")[[1]]
x_All_ASCII <- paste( x_split[which(stri_enc_isascii(x_split))],sep=">",collapse = "" )
x_All_ASCII
}
edit_file[which(enc=="UTF-8")] <- as.character(lapply(edit_file[which(enc=="UTF-8")],remove_Non_ASCII))
cat(paste("去掉非ASCII字符",length(which(enc=="UTF-8")),"项\t"))
}
del_rows <- which(str_detect(edit_file[1:20],"^ xs:schemaLocation|^ xmlns:xs=|^ xmlns="))
edit_file<-edit_file[-del_rows]
writeLines(edit_file,drug_XML_file,useBytes=T)
cat(paste("删除了第",paste(del_rows,collapse = ","),"行\t "))
# 获得网页内容,如果包含空格，则使用 asText = TRUE
doc <- xmlParse(drug_XML_file, useInternal = TRUE,encoding="ASCII")
# 获取 drug title:
Drug_tile <- xmlValue(doc[['/Record/Section/Section/Information[preceding-sibling::TOCHeading[text()="Record Title" and position()=1]]/StringValue']])
# 获得 Canonical SMILES
Canonical_SMILES<- xmlValue(doc[['/Record//Section/Information[preceding-sibling::TOCHeading[text()="Canonical SMILES" and position()=1]]/StringValue']])
# 注意这两种方法结果一样：
# xmltop = xmlRoot(doc)
# getNodeSet(xmltop, '//Section/Information[preceding-sibling::TOCHeading[text()="Canonical SMILES" and position()=1]]/StringValue')
# 获得 CAS：
CAS <- xmlValue(doc[['/Record//Section/Section[preceding-sibling::TOCHeading[text()="Other Identifiers" and position()=1]]/Information/StringValue']])
# 获得 Molecular Weight：
Mol_Weight_value <- xmlValue(doc[['/Record//Section/Information/Table/Row/Cell[preceding-sibling::Cell/StringValue[text()="Molecular Weight" and position()=1]]/NumValue']])
Mol_Weight_unit <- xmlValue(doc[['/Record//Section/Information/Table/Row/Cell[preceding-sibling::Cell/StringValue[text()="Molecular Weight" and position()=1]]/ValueUnit']])
Mol_Weight <- paste(Mol_Weight_value,Mol_Weight_unit,sep=" ")
# 获得 Molecular Formula:
Mol_Formula <- xmlValue(doc[['/Record//Section/Information[preceding-sibling::TOCHeading[text()="Molecular Formula" and position()=1]]/StringValue']])
# 保存到数据框：
A_drug [i,"Drug_tile"] <- Drug_tile
A_drug [i,"Canonical_SMILES"] <- Canonical_SMILES
A_drug [i,"CAS"] <- CAS
A_drug [i,"Mol_Weight"] <- Mol_Weight
A_drug [i,"Mol_Formula"] <- Mol_Formula
# 下载2D结构图片，必须以二进制模式下载文件：
# download(A_drug [i,img_2D_url],Img_2D_file,quiet = TRUE,mode = "wb")
download.file(A_drug [i,img_2D_url],Img_2D_file,mode="wb",quiet = TRUE)
cat(ifelse(file.exists(Img_2D_file),"数据提取成功，图片下载成功!","图片下载失败!"))
cat("\n")
}
}
结果如下：
爬的过程中遇到不少问题，感觉用python爬数据可能更方便一些。但..谁让我爱R(不会)呢...
更多原创精彩视频敬请关注生信杂谈:
网友评论