搜索
查看: 13124|回复: 34

生信编程直播第8题-几个ID转换咯

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-2-24 09:20:31 | 显示全部楼层 |阅读模式
请先学一点基础知识:ID转换大全(http://www.biotrainee.com/thread-862-1-1.html)

不管是什么ID转换,都是找到对应关系,然后match一下即可!
然后是练习题:

首先是affymetrix芯片的探针:运行里面的R代码,得到的变量my_probe就是我们需要转换的探针ID!

[AppleScript] 纯文本查看 复制代码
rm(list=ls())

library("hgu95av2.db")
ls('package:hgu95av2.db')
probe2entrezID=toTable(hgu95av2ENTREZID)
probe2symbol=toTable(hgu95av2SYMBOL)
probe2genename=toTable(hgu95av2GENENAME)

my_probe = sample(unique(mappedLkeys(hgu95av2ENTREZID)),30)

tmp1 = probe2symbol[match(my_probe,probe2symbol$probe_id),]
tmp2 = probe2entrezID[match(my_probe,probe2entrezID$probe_id),]
tmp3 = probe2genename[match(my_probe,probe2genename$probe_id),]

write.table(my_probe,'my_probe.txt',quote = F,col.names = F,row.names =F)
write.table(tmp1$symbol,'my_symbol.txt',quote = F,col.names = F,row.names =F)
write.table(tmp2$gene_id,'my_geneID.txt',quote = F,col.names = F,row.names =F)



如果用perl和python来写,就把对应关系保存到本地文件,或者去其它地方下载对应关系,然后用hash或者dictionary来转换!
ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/  下载对应关系
然后是illumina芯片的 探针:

[AppleScript] 纯文本查看 复制代码
library("illuminaHumanv4.db")
ls('package:illuminaHumanv4.db')

probe2entrezID=toTable(illuminaHumanv4ENTREZID)
probe2symbol=toTable(illuminaHumanv4SYMBOL)
probe2genename=toTable(illuminaHumanv4GENENAME)

my_probe = sample(unique(mappedLkeys(illuminaHumanv4ENTREZID)),30)

probe2symbol[match(my_probe,probe2symbol$probe_id),]
probe2entrezID[match(my_probe,probe2entrezID$probe_id),]
probe2genename[match(my_probe,probe2genename$probe_id),]


其实探针不管是什么平台,都是很简单的!
所有bioconductor支持的芯片平台对应关系:通过bioconductor包来获取所有的芯片探针与gene的对应关系
可以从NCBI的GPL平台下载:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6947
也可以直接加载对应的包!
或者直接去公司的主页下载manifest文件!

最后是基因的转换:
运行下面的R代码,得到的my_symbol_gene和my_entrez_gene就是我们需要转换的ID
[AppleScript] 纯文本查看 复制代码
library("illuminaHumanv4.db")
ls('package:illuminaHumanv4.db') 
my_entrez_gene = sample(unique(mappedRkeys(illuminaHumanv4ENTREZID)),30)
my_symbol_gene = sample(unique(mappedRkeys(illuminaHumanv4SYMBOL)),30) 

library("org.Hs.eg.db")
ls('package:org.Hs.eg.db') 

entrezID2symbol <- toTable(org.Hs.egSYMBOL)

entrezID2symbol[match(my_entrez_gene,entrezID2symbol$gene_id),]
entrezID2symbol[match(my_symbol_gene,entrezID2symbol$symbol),]


讲师可以不用直接拿R里面现成数据,可以秀一秀去其它地方获取对应关系文件!










上一篇:BS-seq 分析流程初学1
下一篇:根据基因名找找不同物种的特异基因!比如-人-老鼠
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

7

主题

25

帖子

678

积分

高级会员

Rank: 4

积分
678
发表于 2017-3-19 22:50:39 | 显示全部楼层
本帖最后由 xotong 于 2017-3-19 23:01 编辑

第八题ID转换作业:原文件:geneid2symbol.txt和geneid.list
解题思路:1.读入有两种ID对应关系的文件,并建立相应可以查询的字典;2.根据提供的其中一种ID,利用已建立好的字典查询出另外对应的ID。

[Python] 纯文本查看 复制代码
#用from ... import ...命令导入需要用到的外部模块中的单个字典对象
from collections import OrderedDict

#先建立一个空的字典
map = OrderedDict()

#用with open 打开有gene ID和genesymbol ID对应关系的文件,以“rt”模式进行操作,并赋值给函数f#(rt模式下,python在读取文本时会自动把windows下的换行符\r\n转换为python内部换行符\n)
with open('C:/Users/Administrator/Desktop/geneid2symbol.txt','rt') as f:
    for line in f:
        line = line.rstrip()
        a = line.split()
#建立字典map[geneid] = genesymbolid
        map[a[0]] = a[1]

#打开gene ID文件,根据上一步建立好的字典查找对应的genesymbol ID
with open('C:/Users/Administrator/Desktop/geneid.list','rt') as f:
    for line in f:
        line = line.rstrip()
#用print函数将查询到的结果在屏幕中打印出来
        print(line+' '+map[line])


输出:
3101 HK3
11343 MGLL
9477 MED20
5551 PRF1
6693 SPN
6511 SLC1A6
1910 EDNRB
1056 CEL
8604 SLC25A12
835 CASP2
10 NAT2


刚开始学习python,根据群主提供的两个文件,主要学习了python中如何建立字典,并进行查询。

回复 支持 1 反对 0

使用道具 举报

6

主题

23

帖子

201

积分

中级会员

Rank: 3Rank: 3

积分
201
发表于 2017-3-8 21:27:06 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-3-8 21:28 编辑

file1
a        2
s        3
c        1
d        4        
file2
a        3
b        2
c        6
d        2
[Shell] 纯文本查看 复制代码
awk 'BEGIN{FS=OFS="\t"}NR==FNR{a[$1]=$2;next}(a[$1]){print a[$1], $0}' file1 file2

2        a        3
1        c        6
4        d        2
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-3-13 23:41:44 | 显示全部楼层
本帖最后由 旭日早升 于 2017-3-14 10:09 编辑

201-python-无名
这一次作业主题是ID转换,就像楼主说的,只要找到对应关系,那就很简单。所以我也奇怪,经过前几次作业这次的编程任务明显简单啊,仔细想想楼主还是涉及到很多点的。首先,ID转换在生信是非常普遍的,尤其是芯片,一般都是probeID,可以对应到gene-symbol,symbol有可以对应到ensemblID或者entrezID等。另外楼主还给我们传达了获取对应关系的几种方式:gene对应关系可以在ncbi平台获取,芯片对应关系可以在GPL各平台、bioconductor或者各公司主页等。可以说给大家指明了一个从获取数据到处理数据的方向,其中我觉得获取数据是尤为重要。不过编程仍是我们的作业了。我就用python转换一个芯片的probe数据吧。
解题思路:1,获取probe和对应关系(楼主R代码);2、读入对应关系,用字典存储; 3、读入需要转换的probe,用list存储;4、用字典查询并输出结果。
我的代码:
[Python] 纯文本查看 复制代码
### Description: Cainiao8_v1, ID transform


##import module
import os
from collections import OrderedDict
##change directory if necessary
os.chdir("./")
##initialization variable
probe2symbol = OrderedDict()
my_probe = []
##read probe2symbol as a dictionary
with open("probe2symbol.txt",'rt') as f:
        for line in f:
                line = line.rstrip()
                lst = line.split()
                if lst[0] not in probe2symbol:
                        probe2symbol[lst[0]] = lst[1]
                else:
                        probe2symbol[lst[0]] += "\t" + lst[1]

##read my probe as a list
with open("my_probe.txt",'rt') as f:
        for line in f:
                line = line.rstrip()
                if line not in my_probe:
                        my_probe.append(line)

##write the matched gene
with open("matched-gene.txt",'wt') as f:
        f.write("probe"+"\t"+"gene-symbol"+"\n")
        for i in range(len(my_probe)):
                probe = my_probe
[i]                gene = probe2symbol[probe][/i]
[i]                f.write(probe+"\t"+gene+"\n")
结果展示:
probe   gene-symbol
40686_at        NMB
38807_at        HSPA4
38129_at        GK
35739_at        MTMR3
35288_at        ARHGEF9


今天晚了,明天看下老师视频。




这次python居然是楼主讲,记得楼主也是学习python不久,但是已经可以讲解python,着实比较厉害。哈哈哈,看到楼主视频就知道群主应该是学习perl高手,对于python的习惯也是不少吐槽。不过我个人因为没有学习过perl,所以我觉得python还是比较有规律的,不过看到楼主经常搞perl的一行代码,如果我之后有余力也想简单学习一下了。下面比较我和楼主的代码吧,因为差不多我就不贴了大家自己敲吧,哈哈。
与楼主代码比较:首先是match文件的读入,我还判断了这个probe是否已经在字典里,其实对于match文件probe一般都是唯一的,但是难免有bug,尤其是非官方的一些文本。第二个是my probe的读取,我是先把它读取到一个列表,而楼主是每读取一行做一次判断和输出,省去了一次循环,更加节省时间。最后感谢楼主花时间学习python和分享。
回复 支持 反对

使用道具 举报

0

主题

1

帖子

27

积分

新手上路

Rank: 1

积分
27
发表于 2017-3-14 09:34:59 | 显示全部楼层
本帖最后由 求知 于 2017-3-14 12:21 编辑

这个题目照猫画虎比较容易实现,就学了点R,perl和python还在学习计划中。就拿R实现一下吧。
解题思路:

1、R的基础知识还是要的,从bioconductor安装软件包。
2、然后按照群主的代码,一步一步实现就可以了(我是在R中建了一个R script)。
3、PS:其实以前我都是用在线工具实现各种ID的转化,包括GO ID,KEGG ID,Reactome ID等等,适合我这种比菜鸟还菜的人。也用R做过,用的是老Yu的clusterProfiler包中的[size=13.333333015441895px]bitr(geneID, fromType, toType, OrgDb, drop = TRUE),想用的自己查度娘。
[Perl] 纯文本查看 复制代码
rm(list=ls())
source("https://bioconductor.org/biocLite.R") #安装bioconductor的包
biocLite("hgu95av2.db")#安装bioconductor的包,安装了可以忽略
library("hgu95av2.db")#加载包
ls('package:hgu95av2.db')#查看包
probe2entrezID <- toTable(hgu95av2ENTREZID)#把probe 对应的Gene ID提出来
probe2symbol <- toTable(hgu95av2SYMBOL)#把probe 对应的Gene Symbol提出来
probe2genename <- toTable(hgu95av2GENENAME)#把probe 对应的Gene Name提出来
my_probe <-  sample(unique(mappedLkeys(hgu95av2ENTREZID)),30)#随机取30个probe
tmp1 = probe2symbol[match(my_probe,probe2symbol$probe_id),]#把30个probe对应的Gene Symbol找到
tmp2 = probe2entrezID[match(my_probe,probe2entrezID$probe_id),]#把30个probe对应的Gene ID找到
tmp3 = probe2genename[match(my_probe,probe2genename$probe_id),]#把30个probe对应的Gene Name找到
write.table(my_probe,'my_probe.txt',quote = F,col.names = F,row.names =F)#把probeID写成一个txt文件
write.table(tmp1$symbol,'my_symbol.txt',quote = F,col.names = F,row.names =F)#把Gene Symbol写成一个txt文件
write.table(tmp2$gene_id,'my_geneID.txt',quote = F,col.names = F,row.names =F)#把Gene ID写成一个txt文件
write.table(tmp3$gene_name,'my_geneName.txt',quote = F,col.names = F,row.names =F)#画蛇添足data <- merge(tmp1,tmp2,by="probe_id")#画蛇添足,把得到的结果合到一起
data <- merge(data,tmp3,by="probe_id")#画蛇添足,把得到的结果合到一起
head(data)#最后看一下结果


这是最后的结果:

   probe_id   symbol   gene_id                                                    gene_name
1    1134_at     TNK2    10188                                           tyrosine kinase, non-receptor, 2
2 1465_s_at   MYBL1      4603                                                  MYB proto-oncogene like 1
3    1696_at     POLB      5423                                                        polymerase (DNA) beta
4 2085_s_at  CTNNA1     1495                                                                  catenin alpha 1
5      235_at    S100B      6285                                            S100 calcium binding protein B
6  31996_at ARFGEF2   10564   ADP ribosylation factor guanine nucleotide exchange factor 2


帖子里老有[size=13.333333015441895px],怎样去掉,知道的告诉一下。
[size=13.333333015441895px]





回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-3-14 16:30:00 | 显示全部楼层
弱弱说一句,bioMart包也可以进行ID转换吧。
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-3-14 16:33:37 | 显示全部楼层
学习最快乐 发表于 2017-3-14 16:30
弱弱说一句,bioMart包也可以进行ID转换吧。

是不是因为不够标准,各种平台和芯片不一样的缘故吧
回复 支持 反对

使用道具 举报

0

主题

1

帖子

101

积分

注册会员

Rank: 2

积分
101
发表于 2017-3-14 20:51:23 | 显示全部楼层
本帖最后由 lifeutopia 于 2017-3-14 20:52 编辑

第八课作业
ID转换主要是probe_id,gene_id,gene_name, symbol之间的转换
首先是随机抽取30个人的基因芯片的probe_id,并将其转换成entrezID、genename和symbol。
编程思路:
1、通过ls()找到包里面和entrezID、genename和symbol相关的对象;
2、通过toTable将其格式化成data.frame的数据形式;
3、通过match进行匹配(match出来的结果是位置信息);
4、再通过数据框的操作将匹配的entrezID、genename和symbol提取出来;
5、将probe_id,gene_id,gene_name, symbol数据进行整合,通过merge循环;
6、输出结果到txt文档
相应的代码即大神的代码
[AppleScript] 纯文本查看 复制代码
#将probe_id对应的symbol、entrezID、genename找到
rm(list=ls())
library("hgu95av2.db")
ls("package:hgu95av2.db")
probe2entrezID=toTable(hgu95av2ENTREZID)
probe2symbol=toTable(hgu95av2SYMBOL)
probe2genename=toTable(hgu95av2GENENAME)

my_probe = sample(unique(mappedLkeys(hgu95av2ENTREZID)),30)

tmp1 = probe2symbol[match(my_probe,probe2symbol$probe_id),]
tmp2 = probe2entrezID[match(my_probe,probe2entrezID$probe_id),]
tmp3 = probe2genename[match(my_probe,probe2genename$probe_id),]

multimerge = function(dat=list(),...){
    if(length(dat)<2)return(as.data.frame(dat))
    mergedat = dat[[1]]
    dat[[1]] = NULL
    for(i in dat){
        mergedat<-merge(mergedat,i,...)
    }
    return(mergedat)
}
tmp4 = multimerge(list(tmp1,tmp2,tmp3))

write.table(my_probe,'my_probe.txt',quote = F,col.names = F,row.names =F)
write.table(tmp1$symbol,'my_symbol.txt',quote = F,col.names = F,row.names =F)
write.table(tmp2$gene_id,'my_geneID.txt',quote = F,col.names = F,row.names =F)
write.table(tmp3$gene_name,'gene_name.txt',quote = F,col.names = F,row.names =F)write.table(tmp4,'ID_4.txt',quote = F,col.names = F,row.names =F)

当然也可以将symbol、entrezID、genename转换成其它的类型,这里似乎只能先转换成probe_id,再转换成其它的,直接互转还需要借助其它的包,没有接触,不甚明白。
仅以此纪念生信练习第一帖。


回复 支持 反对

使用道具 举报

0

主题

5

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-3-14 22:14:54 | 显示全部楼层
本帖最后由 王土土哈 于 2017-4-30 14:32 编辑

322 + python + shell
主要思想就是寻找对应关系
1 shell
(1) join
[Bash shell] 纯文本查看 复制代码
join -1 1 <(sort -k1,1 geneid.list) -2 1 <(sort -k1,1 geneid2symbol.txt) -a 1 -o 1.1 -o 2.2 -e "NULL" >tmp.txt

(2) 用awk
[Bash shell] 纯文本查看 复制代码
cat geneid.list|while read line
do
    awk -v id=$line '{if($1==id)print $1"\t"$2}' geneid2symbol.txt 
done>tmp

(3) 用grep
[Bash shell] 纯文本查看 复制代码
grep -w -f geneid.list geneid2symbol.txt >>tmp

2 python
[Python] 纯文本查看 复制代码
import os
from collections import OrderedDict

os.chdir("E:\\work\\biopython\\biotrainee\\lesson8")
gene_info = OrderedDict()

with open("geneid2symbol.txt",'rt') as f1:
    for line1 in f1:
        lst = line1.rstrip().split(' ')
        gene_info[lst[0]] = lst[1]
 
with open('gene2symbol.txt','wt') as f2:
    print('gene\tsymbol',file=f2)      
    with open("geneid.list",'rt') as f3:
        for line2 in f3:
            id = line2.rstrip()
            if id in gene_info.keys():
                print(id,gene_info[id],sep='\t',file=f2)


目前先解答到这,接下来想通过Python整合所有的基因信息,实现无论提交什么基因信息,都能相对应地输出其他的信息

回复 支持 反对

使用道具 举报

0

主题

16

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-3-15 10:36:31 | 显示全部楼层
本帖最后由 azazelcc 于 2017-3-15 10:41 编辑

生信作业第八节
307-Python-Richard
本来是想学Python的,结果发现R都要忘光了,所以先用R提交作业。虽然寥寥两行,但是也用了大半个晚上。基础差没得说。
为了熟悉一下各种数据结构,画蛇添足的用了数据框。编程思路:
1,加载现成的芯片对应关系
2,随机生成或者读取现有的目标芯片号
3,用现有的对应关系和需要转换的号码进行比对,得到位置关系
4,用位置信息输出总关系表格,得到最终结果
[] 纯文本查看 复制代码
library("hgu95av2.db") #载入关系包

ls('package:hgu95av2.db') #查看包内容

my_probe <- sample(unique(mappedLkeys(hgu95av2ENTREZID)),30) #随机抽样30个作为测试对象

probe2entrezID <- toTable(hgu95av2ENTREZID) 
probe2genename <- toTable(hgu95av2GENENAME)
probe2symbol <- toTable(hgu95av2SYMBOL) #将对应关系提取到列表

df = data.frame(probe2entrezID,probe2genename[["gene_name"]], probe2symbol[["symbol"]]) #综合几个列表作为查询数据框,方便后期更改

result = df[match(my_probe, df$probe_id),] #将测试对象与查询表格匹配,并输出结果到result


write.table(result,'my_result.txt',quote = F,col.names = F,row.names =F) #将结果输出为文件格式


问题1,本来做完就做完了,结果今早再次运行一直报错,提示toTable 不能打开”character“,奇怪了,明明什么都没变。不知道谁能解答一下。

问题2,不知道为什么match输出结果全是匹配到的位置,难度不应该是一个逻辑结果吗?
回复 支持 反对

使用道具 举报

2

主题

4

帖子

63

积分

注册会员

Rank: 2

积分
63
发表于 2017-3-15 21:01:28 | 显示全部楼层
Shell baker
用之前做基因组注释的数据来演示,字段包括old_id, new_id, offical gene symbol, go id, go term。oldId2newId.txt
[Bash shell] 纯文本查看 复制代码
old_id new_id
evm.TU.scaffold337.1    AMOCG03370010516
evm.TU.scaffold337.29   AMOCG03370051016
evm.TU.scaffold337.22   AMOCG03370020000
evm.TU.scaffold337.33   AMOCG03370091516
evm.TU.scaffold225.49   AMOCG02250181616
evm.TU.scaffold225.48   AMOCG02250171616
evm.TU.scaffold225.10   AMOCG02250051016
evm.TU.scaffold225.51   AMOCG02250190216
evm.TU.scaffold225.25   AMOCG02250130315
evm.TU.scaffold225.46   AMOCG02250160216
evm.TU.scaffold225.14   AMOCG02250080711


oldId2genename.txt
[Bash shell] 纯文本查看 复制代码
evm.TU.scaffold10.110   SAV1
evm.TU.scaffold10.111   NIN
evm.TU.scaffold10.137   TAF5L
evm.TU.scaffold10.15    IMPAD1
evm.TU.scaffold10.150   DCTD
evm.TU.scaffold10.18    DARS2
evm.TU.scaffold10.32    DDX46
evm.TU.scaffold10.38    GLE1
evm.TU.scaffold10.58    FBXO5
同样还有id2go的对应文件
把这些信息按照整合成一个大的table并去掉old_id的信息。
[Bash shell] 纯文本查看 复制代码
join <(sort oldId2newId.txt) <(sort oldId2genename.txt)|join - <(sort oldId2go.txt)|cut -f2- > newId_genename_goid_goterm.txt

结果:
[Bash shell] 纯文本查看 复制代码
AMOCG00010031316        CERS1   GO:0000139      Golgi membrane
AMOCG00010031316        CERS1   GO:0003674      molecular_function
AMOCG00010031316        CERS1   GO:0005576      extracellular region
AMOCG00010031316        CERS1   GO:0005783      endoplasmic reticulum
AMOCG00010031316        CERS1   GO:0005789      endoplasmic reticulum membrane
AMOCG00010031316        CERS1   GO:0005794      Golgi apparatus
AMOCG00010031316        CERS1   GO:0006629      lipid metabolic process
AMOCG00010031316        CERS1   GO:0008083      growth factor activity
AMOCG00010031316        CERS1   GO:0016020      membrane
AMOCG00010031316        CERS1   GO:0016021      integral component of membrane
AMOCG00010031316        CERS1   GO:0030148      sphingolipid biosynthetic process
AMOCG00010031316        CERS1   GO:0035690      cellular response to drug
AMOCG00010031316        CERS1   GO:0036146      cellular response to mycotoxin
AMOCG00010031316        CERS1   GO:0043231      intracellular membrane-bounded organe
AMOCG00010031316        CERS1   GO:0046513      ceramide biosynthetic process
AMOCG00010031316        CERS1   GO:0050291      sphingosine N-acyltransferase activit
AMOCG00010031316        CERS1   GO:0051974      negative regulation of telomerase act
AMOCG00010031316        CERS1   GO:0071492      cellular response to UV-A
AMOCG00010031316        CERS1   GO:0072721      cellular response to dithiothreitol
AMOCG00010041616        COPE    GO:0005198      structural molecule activity
AMOCG00010041616        COPE    GO:0005515      protein binding
AMOCG00010041616        COPE    GO:0005622      intracellular
AMOCG00010041616        COPE    GO:0006888      ER to Golgi vesicle-mediated transpor
AMOCG00010041616        COPE    GO:0006890      retrograde vesicle-mediated transport
AMOCG00010041616        COPE    GO:0006891      intra-Golgi vesicle-mediated transpor
AMOCG00010041616        COPE    GO:0030126      COPI vesicle coat

使用grep 就可以match 要找的基因,mylist 包含关注的geneid。
[Bash shell] 纯文本查看 复制代码
grep -Fwf mylist newId_genename_goid_goterm.txt





回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:30 , Processed in 0.048943 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.