搜索
查看: 11152|回复: 34

生信编程直播第六题:下载最新版的KEGG信息,并且解析好

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-1-16 23:45:33 | 显示全部楼层 |阅读模式
很简单一件事,参见我的博客:http://www.bio-info-trainee.com/1188.html
下载得到文本文件,可以看到里面的结构层次非常清楚

C开头的就是kegg的pathway的ID所在行,D开头的就是属于它的kegg的所有的基因
A,B是kegg的分类,总共是6个大类,42个小类


需要写脚本变成:


代码如下:
[Perl] 纯文本查看 复制代码
perl -alne '{if(/^C/){/PATH:hsa(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' hsa00001.keg >kegg2gene.txt


上面得到的基因是ID,pathway也是ID,其实 你可以得到它们的ID与name的对应表格,请发挥自己的代码能力吧!

收集整理了最新的KEGG信息,就是什么基因对应什么通路,什么套路对应那些基因!
就可以去做富集分析了~
http://www.biotrainee.com/thread-409-1-1.html
就可以回答下面两个问题了:
人类有多少基因是有KEGG数据库注释信息的呢?http://www.biotrainee.com/thread-6-1-2.html
人类有多少基因是有GO数据库注释信息的呢?http://www.biotrainee.com/thread-249-1-3.html






上一篇:【R处理芯片数据】--总结
下一篇:缺失值和重复值的处理办法
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-3-6 12:11:18 | 显示全部楼层
本帖最后由 旭日早升 于 2017-3-6 16:19 编辑

201-python-无名

问题解析:问题说明:下载KEGG文件并且做解析。文件格式:keg文件详见楼主帖子http://www.bio-info-trainee.com/1188.html及下载地址http://www.genome.jp/kegg-bin/ge ... brite%2fhsa&length=。问题解析:keg文件是按一定结构(大类-小类-通路-基因)组织的,需要转换成更容易操作做的数据框。数据结构:{pathway1{gene1,gene2,...},pathway2{gene3,gene4,...}}。提高效率:暂时未知。
这一讲主要是keg数据结构的转换,结构转换在数据分析中挺常见的,我们需要解析成我们自己常用的结构。
我的尝试代码:
[Python] 纯文本查看 复制代码
###import module
import re, os
from collections import OrderedDict
###change directory if necessary
os.chdir("./")
###initialization variable
pathway = OrderedDict()
###read kegg file and transform to a dictory
with open("hsa00001.keg", 'rt') as f:
        for line in f:
                if not line.startswith(("C", "D")): # not pathway or gene line
                        continue
                if line.startswith("C"): # pathway line 
                        line = line.rstrip()
                        lst = line.split()
                        pathway_id = lst[1]
                        if pathway_id not in pathway:
                                pathway[pathway_id] = []
                        else:
                                print("Error: duplicated pathway ID")
                        continue 
                if line.startswith("D"): ## gene line
                        line = line.rstrip()
                        lst = line.split(";")
                        gene_id = lst[0].split()[1]
                        if len(lst[0].split()[2:]) > 1: # if gene symbol has more than one string
                                gene_symbol = " ".join(lst[0].split()[2:])
                        else:
                                gene_symbol = lst[0].split()[2]
                        if gene_id not in pathway[pathway_id]:
                                 pathway[pathway_id].append(str(gene_id)+"\t"+gene_symbol)
                        else:
                                print("Error: dup gene id in one pathway")
###output the result to a file
with open("KEGG2gene.txt", 'wt') as f:
        for path, genes in pathway.items():
                for gene in genes:
                        f.write(path+"\t"+gene+"\n")


输出结果:

可以看到第二列和第三列去重复后数目不同,有一些gene id不同但是symbol相同的情况



看了东老师的视频,老师是按照另一种结构解析的。其中gene的name那里提醒了我们很多需要注意的细节,比如有些基因没有简写,是由几个单词外加一些注释组成,我但是只是把分号“;”分开,忽略了分号前还有注释。另外老师用了很多正则匹配,这也是我的一个短板,需要加强学习的内容。重新跑了老师的代码,和自己的结果做了比对。发现老师输出的gene name虽然已经考虑去掉重复,但是用set查看后还是存在重复,python里一个有意思的现象。
重跑老师代码:
[Python] 纯文本查看 复制代码
import os, re 
from collections import OrderedDict
##change directory 
os.chdir("./")
##read the keg file and transform each line 
hsaKEGG = OrderedDict()
with open("hsa00001.keg", 'rt') as f:
	for line in f:
		line = line.rstrip()
		if line.startswith("A"): # the Big class line
			mch = re.search("A<b>(.+)</b>",line) # match Big class name
			className = mch.group(1)
			hsaKEGG[className] = OrderedDict()
		elif line.startswith("B"): # the Small class line
			if line == "B":
				continue
			else:
				mch = re.search("<b>(.+)</b>", line) # match Small class name
				subClass = mch.group(1)
				hsaKEGG[className][subClass] = OrderedDict()
		elif line.startswith("C"): # Pathway line
			mch = re.search("(\d+)\s(.+)", line) # match pathwayID and pathwayName
			pathID = "hsa"+mch.group(1)
			pathName = re.sub("\s\[.+\]", "", mch.group(2)) # remove the annotation ID like "[PATH:hsa01200]"
			pathway = pathID + ":" + pathName
			hsaKEGG[className][subClass][pathway] = [[],[]]
		elif line.startswith("D"): # Gene line
			lst = line.split(";") # split gene and annotation  
			geneInfo = lst[0].split("\t") # split gene which has no abbreviation
			mch = re.match("D\s+(\d+)\s(.+)", geneInfo[0]) # match geneID and geneSymbol
			geneID = mch.group(1)
			gene = mch.group(2)
			hsaKEGG[className][subClass][pathway][0].append(gene)
			hsaKEGG[className][subClass][pathway][1].append(geneID)

## Write to a file
fh = open("hsa00001.cleaned.keg", 'wt')
for ke, val in hsaKEGG.items():
	for subk, subv in val.items():
		for ptwy, geneList in subv.items():
			genes = ";".join(geneList[0])
			geneIDs = ";".join(geneList[1])
			fh.write("\t".join([ke, subk, ptwy, genes, geneIDs])+"\n")

fh.close()

## count the number of pathways and genes
pthwyNum = 0
allGenes = []
with open("hsa00001.cleaned.keg", 'rt') as f:
	for line in f:
		line = line.rstrip()
		lst = line.split("\t")
		if len(lst) > 3:
			pthwyNum += 1
			geneList = lst[-2].split(";")
			allGenes = allGenes + [gene for gene in geneList if gene not in allGenes]

print("Number of non-empty pathways: %d" %pthwyNum)
print("Number of genes in all pathways: %d" %len(allGenes))
print("Number of genes in all pathways: %d" %len(set(allGenes))) # allGenes still has duplication






本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

1

主题

16

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-2-24 17:05:36 | 显示全部楼层
043-Perl+R+shell-涤生
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
print "gene_id\tpath_id\tsymblo\tpath_name\n";
while(<>){
        chomp;
        if(/^C/){
                /C\s+(\d+)\s+(\S+\s\S+)/;
                $path_name = $2;
                $path_id = $1;
        }
        if(/^D/){
                /D\s+(\d+)\s+(\w+)/;
                $gene_id = $1;
                $symblo = $2;
                print "$gene_id\t$path_id\t$symblo\t$path_name\n";
        }
}

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
#

print "gene_id\tpath_id\n";
while(<>){
        chomp;
        if(/^C/){
                /C\s+(\d+)\s+(\S+\s\S+)/;
        #       $path_name = $2;
                $path_id = $1;
        }
        if(/^D/){
                /D\s+(\d+)\s+(\w+)/;
                $gene_id = $1;
        }
        if(defined $gene_id){
        $hash{$gene_id} .= "$path_id|";
        }
}
foreach $id (keys %hash){
        print "$id\t$hash{$id}\n";
}

待完善

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-2-23 21:46:03 | 显示全部楼层
017中大普外 发表于 2017-2-23 21:40
主要是正则表达式和数据结构间嵌套,我还需要加成编程思维,数据嵌套容易搞混[mw_shl_code=python,true]imp ...

俺的主代码是抄讲师的。。。,封装成函数而已。。。
回复 支持 1 反对 0

使用道具 举报

0

主题

4

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-2-17 07:55:39 | 显示全部楼层
本帖最后由 huamixinhua 于 2017-2-17 07:59 编辑

103-python-如梦
第5题看了老师的视频,要补充的知识点很多,最近时间比较紧,跟的有点吃力,等学好知识点后会再补上。
这是第6题,下了Arabidopsis thaliana拟南芥 最新kegg,代码有点冗杂

[Python] 纯文本查看 复制代码
import os,re

f=open(r'C:\Users\zxx\Desktop\ath00001.keg','r')
for line in f:
    sc=re.search("^C.*\[PATH:(.*\d+)\]",line)   #获取kegg基因,每次遇到C开头就更新jy
    if sc:
        jy=sc.group(1)
    sc1=re.search("^D\s*([\dA-Za-z]+) ",line)    #获取每行通路ID(以D开头)
    if sc1:
        print(jy+' '+sc1.group(1))               #将基因与对应的通路ID拼接输出





本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-2-23 21:40:14 | 显示全部楼层
主要是正则表达式和数据结构间嵌套,我还需要加成编程思维,数据嵌套容易搞混
[Python] 纯文本查看 复制代码
import os
import re
from collections import OrderedDict
import sys
os.chdir('F:/生信菜鸟团讲课/生信菜鸟团作业/6')
#file=sys.argv[1]
def KEGG_deal(kegg_file):
    kegg_dic=OrderedDict()
    with open(kegg_file) as f:
        for line in f:
            line=line.rstrip()
            if line.startswith('+D'):
                continue
            if  line.startswith("#"):
                continue
            if line.startswith(";"):
                continue
            if line.startswith('!'):
                continue
            if line.startswith("%"):
                continue
            if line.startswith("A"):                
                distribution=re.search("A<b>(.*)</b>",line).group(1)
                kegg_dic[distribution]=OrderedDict()
            elif line.startswith('B'):
                if line=="B":
                    continue
                else:
                    biology=re.search("<b>(.*)</b>",line).group(1)
                    kegg_dic[distribution][biology]=OrderedDict()#d这个地方使用多维字典,                   
            elif line.startswith('C'):
                mch=re.search("(\d+)\s(.*)",line)#匹配前面的数字
                pathwayID="hsa"+mch.group(1)#第一个括号内的内容                
                pathwayname=re.sub('\s\[.*\]','',mch.group(2))#提取pathway的功能
                pathway=pathwayID+':'+pathwayname  
                kegg_dic[distribution][biology][pathway]=[[],[]]#继续嵌套字典和列表,总觉得有可能搞混,可能我的编程思维不够强
                #{e:{f:{d:{a:b}}}}=[[x],[y]]这种结构了
                
            elif line.startswith('D'):
                line=line.split(';')
                geneid=re.search("D\s+(\d+)\s(.*)",line[0]).group(1)
                genesymbol=re.search("D\s+(\d+)\s(.*)",line[0]).group(2)
                
                
                kegg_dic[distribution][biology][pathway][0].append(geneid)
                kegg_dic[distribution][biology][pathway][1].append(genesymbol)
            
    with open('kegg.txt','w') as f2:
        for k,it in kegg_dic.items():#A层        
            for x,y in it.items():#B层
                for x1,genelist in y.items():#C.D层
                    geneidx=";".join(genelist[0])
                    genenamex=';'.join(genelist[1])
                    f2.write("\t".join([k,x,x1,geneidx,genenamex])+'\n')                    
            
KEGG_deal("hsa00001.keg")        


下图是我做的结果

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-2-25 12:01:02 | 显示全部楼层
本帖最后由 xxnku 于 2017-2-26 12:12 编辑

[Python] 纯文本查看 复制代码
import re
import sys
from collections import OrderedDict
args = sys.argv
hsaKEGG = OrderedDict()

with open(args[1]) as f:
    for line in f:
        line = line.rstrip()
        if line.startswith('A'):
            mch = re.search('A<b>(.+)</b>',line)
            className = mch.group(1)
            hsaKEGG[className] = OrderedDict()

        elif line.startswith('B'):
            if not line == 'B':
                mch = re.search('B\s+<b>(.+)</b',line)
                subclass = mch.group(1)
                hsaKEGG[className][subclass] = OrderedDict()

        elif line.startswith('C'):
            mch = re.search('(\d+)\s(.+)',line)
            pathID = 'hsa' +  mch.group(1)
            pathName = re.sub('\s\[.+\]','',mch.group(2))
            #这个替换不太懂,有的C开头的行后面没有[PATH:******]   
            pathway = pathID + ':' + pathName
            hsaKEGG[className][subclass][pathway] = []  

        elif line.startswith('D'):
            lst = line.split(';')
            mch = re.search('D\s+(\d+)\s(.+)',lst[0])
            geneID = mch.group(1)
            geneName = mch.group(2)
            gene = geneID + ':' + geneName               
            hsaKEGG[className][subclass][pathway].append(gene)

with open('hsa00001_cleaned.keg','w') as f:
    for k,v in hsaKEGG.items():
        for subk,subv in v.items():
            for ptwy,genelist in subv.items():
                genes = ','.join(genelist)
                f.write('\t'.join([k,subk,ptwy,genes]) + '\n')


pthwyNum = 0
allGenes = []
with open('hsa00001_cleaned.keg') as f:
    for line in f:
        line = line.rstrip()
        lst = line.split('\t')
        if len(lst) > 3:
            pthwyNum += 1
            genelist = lst[-1].split(',')
            for gene in genelist:
                geneID = gene.split(':')[0]
                if geneID not in allGenes:
                    allGenes.append(geneID)
print ("Number of non-empty patways: %d" %pthwyNum)
print ("Number of genes in all pathways: %d" %len(allGenes))


谢谢Python-东
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-2-25 16:06:09 | 显示全部楼层
034-perl-shell

输入文件:hsa00001.keg
输出文件格式 : 第一列 大类;第二列 小类;第三列 通路id以及名称;第四列 通路下基因id以及symbol;第五列通路下K号以及K号对应的名称

perl代码如下:

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $in = "hsa00001.keg";
my $out = "reslut.txt";
open (my $fh, $in) or die "Cannot open $in";
open (my $fh_out, ">$out") or die "Cannot write to $out";
my ($A,$B,$C);
my (@gene_all,@k_id_all);

while(<$fh>){
	chomp;
	next unless($_ =~ /^\w+/);
	if ($_ =~ /^A<b>(\S+)<\/b>/){
		$A = $1;
		next; 
	}
	next if ($_ =~ /^B$/);
	if ($_ =~ /^B\s+<b>(\S+)<\/b>/){
		$B = $1;
		next;
	}
	if ($_ =~ /^C\s+(\d+)\s+(.*?)\s+\[.*\]/){
		$C = $1.":".$2;
		next unless (@gene_all);
		next unless (@k_id_all);
		print $fh_out "$A\t$B\t$C\t";
		my $tmp = join "\|", @gene_all;
		print $fh_out $tmp."\t";
		@gene_all =();
		my $tmp2 = join "\|", @k_id_all;
		print $fh_out $tmp2."\n";
		@k_id_all = ();
		next;
	}
	if ($_ =~ /^D/){
		my @array = split /\t/, $_;
		my $gene = $1 if ($array[0] =~ /.*?(\d+\s+\w+)/);
		my $K_id = $1 if ($array[1] =~ /(.*?);/);
		push @gene_all, $gene;
		push @k_id_all, $K_id;
	}

}



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

2

主题

34

帖子

773

积分

高级会员

Rank: 4

积分
773
发表于 2017-2-26 01:03:52 | 显示全部楼层
本帖最后由 x2yline 于 2017-5-28 12:38 编辑

077-R-python-shell x2yline

数据下载时间:2017-02-26 00:22

python代码如下:
[Python] 纯文本查看 复制代码
import os
os.chdir(r'E:\r\biotrainee_demo\class6')
#数据下载时间:2017-02-26 00:22
#数据下载地址:[url=http://www.genome.jp/kegg-bin/download_htext?htext=hsa00001&format=htext&filedir=]http://www.genome.jp/kegg-bin/do ... rmat=htext&filedir=[/url]
def analysis_kegg(file_path):
    kegg_dict={}
    with open(file_path,'r') as f:
        for line in f:
            if line.startswith('C'):
                try:
                    path_num = line.split()[1]
#                    if len(path_num) > 10:
#                        print(line)
                    kegg_dict[path_num]= []
                except:
                    print('Warning: '+line)
            elif line.startswith('D'):
                try:
                    gene_name = line[:line.find(';')].split()[-1]
                    kegg_dict[path_num].append(gene_name)
                except:
                    print('Warning: '+line)
    non_empty_path = []
    all_gene = []
    for keys in kegg_dict.keys():
        if kegg_dict[keys]:
           non_empty_path.append(keys)
           all_gene += kegg_dict[keys]
    print('The number of pathways in kegg is {}\nThe number of gene in kegg is {}'.format(len(set(non_empty_path)),len(set(all_gene))))

file_path = 'hsa00001.keg'
analysis_kegg(file_path)

#结果The number of pathways in kegg is 322
#The number of gene in kegg is 7217





回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-2-26 10:00:55 | 显示全部楼层
本帖最后由 chapman 于 2017-3-8 09:46 编辑

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

f = open('hsa00001.keg','r')
kegg = OrderedDict()
for line in f:
    line = line.rstrip()
    if line.startswith('A'):
        mch = re.search(r'<b>(.+)</b>',line)
        class_name = mch.group(1)
        kegg[class_name] = OrderedDict()

    elif line.startswith('B'):
        if line == "B":
            continue
        else:
            mch = re.search(r'<b>(.+)</b>',line)
            sub_name = mch.group(1)
            kegg[class_name][sub_name] = OrderedDict()

    elif line.startswith('C'):
        mch = re.search('(\d+)\s(.+)',line)
        pathID = "has"+mch.group(1)
        pathname = re.sub(r'\[.+\]',"",mch.group(2))
        pathway = pathID + ":" +pathname
        kegg[class_name][sub_name][pathway] =[[],[]]

    elif line.startswith('D'):
        lst = line.split(';')
        geneinfo = lst[0].split('\t')
        mch = re.search(r'D\s+(\d+)\s(.+)',geneinfo[0])
        geneID = mch.group(1)
        gene = mch.group(2)

        kegg[class_name][sub_name][pathway][0].append(gene)
        kegg[class_name][sub_name][pathway][1].append(geneID)


fo = open('kegg_inf','w')
for ke,val in kegg.items():
    for subk , subv in val.items():
        for ptwy,genelist in subv.items():
            genes = ";".join(genelist[0])
            geneIDs = ";".join(genelist[1])
            fo.write('\t'.join([ke , subk,ptwy,genes,geneIDs])+"\n")

f.close()


pathwy_num = 0
allgene =[]
with open('kegg_inf','r') as f:
    for line in f:
        line = line.rstrip()
        lst = line.split('\t')
        if len(lst) > 3:
            pathwy_num += 1
            genelist = lst[-2].split(';')
            allgene = allgene + [gene for gene in genelist if gene not in allgene]
print('number of non-empty pathways:%d'% pathwy_num)
print('number of gene in all pathways:%d'% len(allgene)) 
            
            
        
        
    



number of non-empty pathways:317
number of gene in all pathways:7227


回复 支持 反对

使用道具 举报

6

主题

23

帖子

201

积分

中级会员

Rank: 3Rank: 3

积分
201
发表于 2017-2-26 20:31:57 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-2-26 20:33 编辑

123-python-qin
根据要求将GO数据整理了一下,发现只是储备太少了,基础知识用的时候就全乱掉了,还是要多联系哈
[AppleScript] 纯文本查看 复制代码
#!/bin/python
from collections import OrderedDict
import re
import os
hsaGO = OrderedDict()
with open('gene2go') as f:
        for line in f:
                line = line.rstrip()
                lst = line.split('\t')
                if line.startswith('#'):
                        continue
                if lst[7] == 'Component':
                        GO_ID = lst[2]
                        GO_term = lst[5]
                        pathway = GO_ID+':'+GO_term
                        GeneID = lst[1]
                        if pathway not in hsaGO:
                                hsaGO[pathway] = []
                        hsaGO[pathway].append(GeneID)
fn = open('gene2go.cleaned','wt')                        
for pathway ,value in hsaGO.items():
        value = list(set(value))
        geneIDs = ','.join(value)
        #print pathway,geneIDs
        fn.write('\t'.join([pathway,geneIDs])+'\n')
fn.close()
结果如下,只有Component部分。



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-17 07:12 , Processed in 0.043671 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.