天天看點

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的算法,覺得人家好厲害,可以想到那麼低計算度的方式。之後也要多研究研究算法了。

繼續閱讀