天天看点

2021.01.05丨根据基因名称拼接表达量与相关注释

这一步是在进行最后的数据汇总工作中用到的,将基因的count与FPKM值和基因注释的结果组合在一起,得到一个完整的数据。方便客户进行后续研究。算法与之前那篇基因ID匹配注释文本一文相似,用了两个for循环嵌套进行比对,O=n²,在此也希望能够抛砖引玉,得到大神指点。

输入文件:

anno.DEG.txt

2021.01.05丨根据基因名称拼接表达量与相关注释

all.anno.xls #这里用的Editplus打开

2021.01.05丨根据基因名称拼接表达量与相关注释

本来之前我对all.anno.xls的geneID已经处理过了,但是正好遇到ftp出问题,无法下载最新文件,就将就前两天的结果进行处理

#读取文件列
anno_file = open('C:/Users/bbplayer/Downloads/all.anno.xls','r')
diffgene_file = open('C:/Users/bbplayer/Downloads/anno.DEG.txt','r')
genome_line = anno_file.readlines()
diffgene_list = diffgene_file.readlines()
#保留首行
title = diffgene_list[0]
#设置输出文件名
newfile_name = 'anno.DEG2'
desktop_path = 'C:/Users/bbplayer/Downloads/'
file_path = desktop_path+newfile_name+'.txt'
file = open(file_path,'w') #打开文件名
file.write(title) #输入首行
for count_fpkm in diffgene_list:
    count_fpkm = count_fpkm.replace("\n","") #替换换行符
    #print(count_fpkm) 
    gene_ID = count_fpkm.split('\t') #根据分隔符进行分段
    #print(gene_ID[0])
    for line in genome_line:
        str_line = str(line)
        anno_ID = str_line.split(':')
        anno_line = str_line.split('\t',1) #根据分隔符进行分成2段
        #print(anno_line)
        if  anno_ID[0] == gene_ID[0]:
            num = len(anno_line) #注释为空时num=1
            if num !=2:
                continue
            else:
                #print(num)
                file.write(count_fpkm + "\t" + anno_line[1])
        else:
            continue
#关闭文件
file.close()
diffgene_file.close()
anno_file.close()
           

tips:里面注释掉的print()用来测试输出文本,最简单直白的测试方式。

结果展示(只截了中间结合起来的部分,即表达下调down和GO注释连接处):

2021.01.05丨根据基因名称拼接表达量与相关注释

其实还可以这样,如果都用excel打开,然后对geneID统一升序或者降序,是否能直接对应上geneID,直接把注释文件粘贴到FPKM值后面。当然,这种办法一个是手动不方便,另一个是数据量大,中间计算复杂,不一定能保证完全匹配。

最近在看BLAST的算法,觉得人家好厉害,可以想到那么低计算度的方式。之后也要多研究研究算法了。

继续阅读