搜索
楼主: Jimmy

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

[复制链接]

0

主题

5

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-3-10 21:26:41 | 显示全部楼层
本帖最后由 chzhniu 于 2017-3-10 21:29 编辑

下载KEGG Orthology (KO) -Homo sapiens (human)对结果进行查看,结果如下(图1):
1.jpg
1
对以上结果进行解析:
1、 输出pathway及每个pathway下对应的基因号(图2);
[Perl] 纯文本查看 复制代码
perl -lane '{if (/^C/){/PATH:(hsa\d+)/;$path=$1;}elsif (/^D/){print "$path\t$F[1]\n";}}' hsa00001.keg >1.txt
2.jpg
2
2、 输出pathway及每个pathway下对应的基因名称(图3);
[Perl] 纯文本查看 复制代码
perl -lane '{if (/^C/){/PATH:(hsa\d+)/;$path=$1;}elsif (/^D/){print "$path\t$F[2]\n";}}' hsa00001.keg >2.txt
3.jpg
3
3、 输出pathway及每个pathway下对应基因的KO值(图4);
[Perl] 纯文本查看 复制代码
 perl -lane '{if (/^C/){/PATH:(hsa\d+)/;$path=$1;}elsif (/^D/){/(K\d{5})/; print "$path\t$1\n";}}' hsa00001.keg >KO.txt
4.jpg
4
结果中比较奇怪的是有一段基因没有对应的pathway号(图5),
5.jpg
5
查看原文件后发现,有些没有标注[PATH,如Autophagy-yeastMitophagy-animal等(图6):
6.jpg
6
4、 如果不用单行代码,只需对如下代码做适当修改即可
[Perl] 纯文本查看 复制代码
while(<>){
       chomp;
       @tmp=split;
       if (/^C/){
              /PATH:(hsa\d+)/;[/align][align=left]              $path=$1;
       }elsif (/^D/){
              /(K\d{5})/;
              print "$path\t$1\n";
              #print "$path\t$tmp[1]";
              }
       }
结果如下(图7):
7.jpg
7
5、 对于没有pathway编号的通路,实际上前边有编号(图6):
6
所以,换中方式,取前边的编号及对应的基因:
代码如下:
[Perl] 纯文本查看 复制代码
while(<>){
       chomp;
       @tmp=split;
       if (/^C/){
              #/PATH:(hsa\d+)/;[/align][align=left]              #$path=$1;
              $path= $tmp[1]
       }elsif (/^D/){
              /(K\d{5})/;
              #print "$path\t$1\n";
              print "$path\t$1\n";
              }
       }
结果为(path_KO2.txt,图8):
8.jpg
8
发现,有些通路没有对应基因,如上图中显示的,04060就有一个基因没有KO编号,对应的基因为PLEKHO2(图9),
9.jpg
9
同时,还发现,有些通路没有对应的基因,如图9中的BR:hsa04052,所以以上代码也没有将这种没有对应基因的代谢通路取出来(8中没有04052这一通路),所以,以上代码有缺陷,还需优化(OMG……)。
6、 如果不考虑没有基因的通路,统计pathway数目及基因数目,结果如下(图10):
10.jpg
10
基因数目(不同KO数目),结果如下(图11):
11.jpg
11

回复 支持 反对

使用道具 举报

0

主题

4

帖子

41

积分

新手上路

Rank: 1

积分
41
发表于 2017-3-15 01:59:55 | 显示全部楼层
先照这panda姐的脚本写出来,
#!/usr/bin/perl
use strict;
my ($line,$path_id,@gene,$gene_id,$gene_name);
my (%path_count,%gene_count,$path_num,$gene_num,%all_path);
open IN,"hsa00001.keg";
open OUT1,">kegg_parse.txt";
open OUT2,">gene2all_path.txt";
while($line=<IN>){
        chomp($line);
        if($line=~s/^C\s+//){
                $path_id=(split /\s/,$line)[0];
        }elsif($line=~s/^D\s+//g){
                $path_count{$path_id}+=1;
                @gene = split /\s|;/,$line;
                $gene_id =$gene[0];
                $gene_name=$gene[1];
                $gene_count{$gene_id}+=1;
                push @{$all_path{$gene_id}},$path_id;
                print OUT1 $path_id,"\t",$gene_id,"\t",$gene_name,"\n";
        }
}
$path_num=keys %path_count;
$gene_num=keys %gene_count;
print "kegg_parse data stored at kegg_parse.txt\n";
print "Number of pathway : ",$path_num,"\n";
print "Number of gene_annotation:".$gene_num,"\n";
foreach my $gene_id(sort keys %all_path){
        print OUT2 $gene_id,"\t",join ("|",@{$all_path{$gene_id}}),"\n";
}
单行命令得到,pathway_ID \t gene_ID \t gene_name
perl -lane '{if($F[0] eq "C"){$path=$F[1]}elsif($F[0] eq "D"){$gene_id=$F[1];$gene_name=$F[2];};print "$path\t$gene_id\t$gene_name";}' hsa00001.keg | perl -p -e 's/;//'
单行命令得到 gene_ID pathway|pathway|pathway
perl -lane '%hash;{if($F[0] eq "C"){$path=$F[1]}elsif($F[0] eq "D"){$gene_ID=$F[1];$hash{$gene_ID} .= $path."|";}foreach $gene_ID (sort keys %hash ){print "$gene_ID\t$hash{$gene_ID}";}}' hsa00001.keg
shell命令注释的通路数
grep "^C\b" hsa00001.keg |wc -l
shell命令注释的基因数
grep "^D\b" hsa00001.keg |awk '{print $2}' | sort -u | wc -l
回复 支持 反对

使用道具 举报

0

主题

6

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2017-3-16 13:34:07 | 显示全部楼层
本帖最后由 highwind 于 2017-3-16 13:36 编辑

小板凳坐起,第5题比较花时间,晚点来做,先做第6题吧
Q:数据来源:http://www.genome.jp/kegg-bin/do ... rmat=htext&filedir=
A:问题拆解1(看看数据):

长得有点像html格式,不过处理起来看起头的字母就好

问题拆解2(想想算法)
逐行读入逐行记录
累加去重统计

问题拆解3(算法实现)
[Python] 纯文本查看 复制代码
import time
import re

starttime = time.clock()

def import_file(path):
    f = open(path,'r')
    return(f)

def close_file(file):
    file.close()
    
def cal(file):
    Cat,Catsub,Path,Gene = set(),set(),set(),set()
    Dl=list()
    for line in file:
        if(re.match(r'^A',line)):
            A = re.split(r'<b>|</b>',line)[1]
            Cat |= set([A])
        elif(re.match(r'^B  <b>',line)): #防止把B开头的空行统计了,不然下面的split会错
            B = re.split(r'<b>|</b>',line)[1]
            Catsub |= set([B])
        elif(re.match(r'^C',line)):
            C = re.split(r'    | \[',line)[1]
            Path |= set([C])
        elif(re.match(r'^D',line)):
            D = re.split(r'      |;',line)[1] #这里发现ID和gene名并不是一一对应的
            Gene |= set([D])
            Ds = re.split(r' ',D)
            Dl.append(Ds)          
        else:
            continue
    Catsub -= set(['Overview']) #这个到底算不算子类还真没搞懂?
    print("There are {} categories, {} sub-categories, {} pathways and {} annotated gene individually in KEGG"\
      .format(len(Cat),len(Catsub),len(Path),len(Gene)))
    return(Dl)

def find_dup(ls):#这个查找列表内重复元素的函数让我想了好半天,不知道有没有什么现成好用的模块?
    Dd = dict()
    j0 = 0 #方便后面反过来j0=1,j1=0去试试是不是有单ID对应多gene的情况【然而并没有】
    j1 = 1
    for i in ls:#循环构造一个嵌套字典,用key的唯一性和set的唯一性集合元素
        if(i[j1] in Dd.keys()): #这个判断模块也用到了go文件的编程中
            Dd[i[j1]] |= set([i[j0]])
        else:
            Dd[i[j1]]=set([i[j0]])
    print('Following gene have more than one ID')
    for i in Dd:
        if(len(Dd[i])>1):
            print('{:10s}:{}'.format(i,Dd[i]))
        else:
            continue

path = 'hsa00001.keg' 
file = import_file(path)
find_dup(cal(file))
close_file(file)

endtime=time.clock()
print("本次计算用时{:.3f}秒".format(endtime - starttime))  
计算结果:
There are 6 categories, 43 sub-categories, 473 pathways and 7249 annotated gene individually in KEGG
Following gene have more than one ID
putative  :{'100288562', '101930111', '101929627', '101929601'}
CRYAA     :{'102724652', '1409'}
ICOS      :{'29851', '102723996'}
olfactory :{'105369274', '102723532'}
本次计算用时0.776秒
[Python] 纯文本查看 复制代码
import time
import re

starttime = time.clock()

def import_file(path):
    f = open(path,'r')
    return(f)

def close_file(file):
    file.close()
    
def cal(file):
    gene = dict()
    for line in file:
        info = re.split(r'\s',line)
        if(info[0]=='9606'): #homosapine的taxid=9606的觉悟还是要有的!
            if(info[1] in gene.keys()):
                gene[info[1]] |= set([info[2]])
            else:
                gene[info[1]] = set([info[2]])
    print('{} gene have been annotated in GO'.format(len(gene)))

path = 'gene2go' 
file = import_file(path)
cal(file)
close_file(file)

endtime=time.clock()
print("本次计算用时{:.3f}秒".format(endtime - starttime)) 
计算结果:
18902 gene have been annotated in GO
本次计算用时15.853秒

问题拆解4(算法优化)
。。。

其他 :
理解了一下go文件的格式:ftp://ftp-trace.ncbi.nih.gov/gene/DATA/README


回复 支持 反对

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-3-27 15:19:39 | 显示全部楼层
本帖最后由 朝曦曦 于 2017-3-27 23:36 编辑

学员205 python
我想知道讲师的hsa00001.clean.keg是怎么下载的啊,我只看到了hsa00001.keg
代码主要是重复讲师的, 又参考了同学们的,最后终于跑出来啦!
[Python] 纯文本查看 复制代码
#file=sys.argv[1]
KEGG_file="hsa00001.keg"
def KEGG_deal(KEGG_file):
    hsaKEGG = 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'):
            mch = re.search('A<b>(.+)</b>',line)
            className = mch.group(1)
            hsaKEGG[className] = OrderedDict()
            
        elif line.startswith('B'):
            if line == 'B':
                continue
            else:
                mch = re.search('<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))
            pathway = pathID + ':' + pathName
            
            hsaKEGG[className][subClass][pathway] = [[],[]]
        
        elif line.startswith('D'):
            lst = line.split(';')
            geneInfo = lst[0].split('\t')
            mch = re.match('D\s+(\d+)\s(.+)', geneInfo[0])
            geneID = mch.group(1)
            gene = mch.group(2)
            
            hsaKEGG[className][subClass][pathway][0].append(gene)
            hsaKEGG[className][subClass][pathway][1].append(geneID)

with open('KEGG_hsa00001.txt', 'w') as f2:
    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])
                f2.write('\t'.join([ke, subk, ptwy, genes, geneIDs])+'\n')


KEGG_deal("hsa00001.keg")

然后结果长这样

Screen Shot 2017-03-27 at 5.35.22 PM.png
回复 支持 反对

使用道具 举报

2

主题

9

帖子

235

积分

中级会员

Rank: 3Rank: 3

积分
235
发表于 2017-3-29 22:22:41 | 显示全部楼层
本帖最后由 wangfangce 于 2017-3-29 22:24 编辑

206-Perl
看了panda姐的视频,感谢分享的一维hash数组,非常实用的技巧。
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
open IN,"/home/wangfangce/biosoft/kegg/hsa00001.keg";
open OUTPUT,">kegg_prase.txt";
open OUTPUT2,">gene2all_path.txt";
while($line=<IN>){
chomp($line);
if($line=~s/^C\s+//){
$path_id=(split /\s/,$line)[0];
}elsif($line=~s/^D\s+//g){
$path_count{$path_id}+=1;
@gene=split /\s|;/,$line;
$gene_id=$gene[0];
$gene_name=$gene[1];
$gene_count{$gene_id}+=1;
push @{$all_path{$gene_id}},$path_id;
print OUTPUT $path_id,"\t",$gene_id,"\t",$gene_name,"\n";
}
}
$path_num=keys %path_count;
$gene_num=keys %gene_count;
print "kegg_parse data stored at kegg_prase.txt\n";
print "Number of pathway: ",$path_num,"\n";
print "Number of gene_annotation: ",$gene_num,"\n";

foreach my $gene_id(sort keys %all_path){
 print OUTPUT2 $gene_id,"\t",join("|",@{$all_path{$gene_id}}),"\n";
}
回复 支持 反对

使用道具 举报

0

主题

6

帖子

291

积分

中级会员

Rank: 3Rank: 3

积分
291
发表于 2017-4-15 10:27:29 | 显示全部楼层
思路:将所需信息用正则提取,分层次存储到字典中,迭代写入新的文件。
脚本:
[Python] 纯文本查看 复制代码
import re
from collections import OrderedDict

hsaKEGG = OrderedDict()
with open('E:\Project\python_work\work6\hsa00001.keg','r') as f:
    for line in f:
        line = line.strip()
        if line.startswith('A'):
            mch = re.search(r'A<b>(.+)</b>', line)
            classname = mch.group(1)
            hsaKEGG[classname] = OrderedDict()
        elif line.startswith('B'):
            if line == 'B':
                continue
            elif line.startswith('B'):
                mch = re.search(r'<b>(.+)</b>', line)
                subclass = mch.group(1)
                hsaKEGG[classname][subclass] = OrderedDict()
        elif line.startswith('C'):
            mch = re.search(r'C\s+(\d+)\s(.+)', line)
            pathID = 'hsa' + mch.group(1)
            pathname = re.sub(r'\s\[.+\]', '', mch.group(2))
            pathway = pathID + ':' + pathname
            print pathway
            hsaKEGG[classname][subclass][pathway] = [[], []]
        elif line.startswith('D'):
            line = line.split(';')
            geneinfo = line[0].split('\t')[0]
            mch = re.search(r'D\s+(\d+)\s(.+)', geneinfo)
            geneID = mch.group(1)
            hsaKEGG[classname][subclass][pathway][0].append(geneID)
            genesymbol = mch.group(2)
            hsaKEGG[classname][subclass][pathway][1].append(genesymbol)

with open('E:\Project\python_work\work6\hsa00001.cleaned.keg','w') as f_out:
    for cn,val in hsaKEGG.items():
        for sn,val1 in val.items():
            for pw,val2 in val1.items():
                geneIDs = ';'.join(val2[0])
                genesymbols = ';'.join(val2[1])
                f_out.write('\t'.join([cn, sn, pw, geneIDs, genesymbols]) + '\n')

小结:目前已基本掌握整理统计文本文件的思路,还需多加练习,做到熟练运用。
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-4-18 22:46:18 | 显示全部楼层
本帖最后由 sxuan 于 2017-4-19 10:24 编辑

111111111111
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-4-18 22:47:15 | 显示全部楼层
本帖最后由 sxuan 于 2017-4-19 10:52 编辑

2222222222222222222
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-4-19 10:59:38 | 显示全部楼层
本帖最后由 sxuan 于 2017-4-19 11:01 编辑

15行python代码解决 23333333333333333333333333   
[Python] 纯文本查看 复制代码
import re

with open(r'E:\迅雷下载\hsa00001.keg', 'r') as f:
    # filter and keep the line startwith 'C' OR 'D'
    l = [line for line in f.readlines() if line.startswith('C')
         or line.startswith('D')]

# get the boundry index of every class 
C_index = [idx for idx, value in enumerate(l) if value.startswith('C')]

# match the ID of each element in l
l_ID=[re.search(r'[0-9]+', element).group() for element in l]
# print(l_ID)

# get every pathway ID 
pathway_ID = [l_ID[index] for index in C_index]

gene_ID = []
# for every for cycle make the genes in each class to be a sublist of gene_ID 
for i in range(len(C_index)):[i]
    try:
        gene_ID.append(l_ID[(C_index[i]+1): C_index[i+1]])
    except:
        gene_ID.append(l_ID[C_index[i]+1:])

keggDict = dict(zip(pathway_ID, gene_ID))

# obtain the number of genes of every class
gene_number =[len(l) for l in gene_ID]
kegg_numberDict=dict(zip(pathway_ID,gene_number))
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-4-21 20:11:40 | 显示全部楼层
朝曦曦 发表于 2017-3-27 15:19
学员205 python
我想知道讲师的hsa00001.clean.keg是怎么下载的啊,我只看到了hsa00001.keg
代码主要是重复 ...

你这with open 这个嵌套的if 条件写的估计老师看了想打人,前面这几个continue的没必要全部写出来,你要想写出来可以写成一句,中间用or连接,优雅的写法是直接把你后面想要的几个写出来,例如startwith B,C,D的这几句写出来,最后写个else: continue。就ok了,要优雅~
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-20 21:10 , Processed in 0.038502 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.