目錄

寫在前面

功能注釋數(shù)據(jù)庫介紹

方法一: 以KEGG的注釋結(jié)果為主, 篩選出每個(gè)品種包含的特異通路及基因

方法二: 利用pfam的注釋結(jié)果, 過濾掉每個(gè)品種具有相同結(jié)構(gòu)域的基因, 對(duì)于剩下的基因, 利用其swissprot id 在DAVID數(shù)據(jù)庫中查找他們的功能信息, 單個(gè)進(jìn)行篩選

最后

寫在前面

事情是這樣的, 前段時(shí)間接到師兄布置的一個(gè)任務(wù), 需要從9個(gè)品種特異區(qū)段的總共4000多個(gè)基因中找到跟特異性狀相關(guān)的基因, 而且最好是單個(gè)找, 盡量不要漏掉每個(gè)信息, 然后我就拿到了這些基因的功能注釋文件, 先上圖吧

功能注釋文件
相信接觸過基因組注釋的同學(xué)對(duì)這個(gè)文件印象都很深刻, 這是折騰好多天整合多個(gè)數(shù)據(jù)庫得到的功能注釋結(jié)果, 這個(gè)結(jié)果中包含了GO, NCBI, Swissprot, KEGG, IPRSCAN和Pfam六個(gè)數(shù)據(jù)庫
這篇筆記沒有什么流程性可言, 只是我為了更快更好地完成任務(wù)想出來的小辦法, 如果哪位大神覺得有問題, 希望能給我指出, 大家共同進(jìn)步, 跪謝了~

功能注釋數(shù)據(jù)庫介紹

GO數(shù)據(jù)庫

Gene Ontology (GO)為了能夠使對(duì)各種數(shù)據(jù)庫中基因產(chǎn)物的功能描述相一致，發(fā)展了具有三級(jí)結(jié)構(gòu)的標(biāo)準(zhǔn)語言（ontologies）。根據(jù)基因產(chǎn)物的相關(guān)分子功能，生物學(xué)途徑，細(xì)胞學(xué)組成而給予定義，無物種相關(guān)性。目前，GO的定義法則已經(jīng)在多個(gè)合作的數(shù)據(jù)庫中使用，使得這些數(shù)據(jù)庫在查詢時(shí)具有很高的一致性
在基因表達(dá)分析中GO數(shù)據(jù)庫如何發(fā)揮它的作用呢？很簡單，將轉(zhuǎn)錄組數(shù)據(jù)引入GO注釋，可以揭示某一特定組具有相似表達(dá)模式基因的功能，因?yàn)楣脖磉_(dá)的基因可能編碼在同一個(gè)生物過程中出現(xiàn)的基因產(chǎn)物，或定位于同一個(gè)細(xì)胞組分。例如，某個(gè)未知的A基因和一些已經(jīng)被注釋到GO數(shù)據(jù)庫中某一生物學(xué)過程的B類基因共表達(dá)，那么這個(gè)未知的A基因很可能與B類基因一樣在同一個(gè)生物學(xué)過程中發(fā)揮功能
實(shí)用貼—注釋數(shù)據(jù)庫介紹之GO、KEGG數(shù)據(jù)庫

NCBI數(shù)據(jù)庫

ncbi就不過多介紹了
我感覺在注釋結(jié)果中, ncbi的結(jié)果是最多的, 同時(shí), 也可能是最不準(zhǔn)確的

Swissprot數(shù)據(jù)庫

Swissprot數(shù)據(jù)庫中的基因是經(jīng)過功能驗(yàn)證的, 我覺得參考價(jià)值應(yīng)該是最高的
正相反, uniprot數(shù)據(jù)庫中的另外一個(gè)TrEMBL數(shù)據(jù)庫包含的是沒有經(jīng)過功能驗(yàn)證的基因
每個(gè)序列條目由核心數(shù)據(jù)(Core Data)和注釋數(shù)據(jù)(Annotation)組成。核心數(shù)據(jù)包括序列、參考文獻(xiàn)和序列的生物來源，而注釋數(shù)據(jù)則描述了：①蛋白質(zhì)的功能；②蛋白質(zhì)的翻譯后加工修飾，如糖基化(Carbohydration)、磷酸化(Phosphorylation)、乙?；?Acetylation)、GPI錨定(GPI-anchor)等；③結(jié)構(gòu)域( Domains)和結(jié)合位點(diǎn)(Binding Sites)，如鈣結(jié)合區(qū)(Calcium Binding Regions)、ATP結(jié)合位點(diǎn)(ATP-Binding Sites)、鋅指結(jié)構(gòu)域(Zinc Fingers)、同源盒(Homeobox)、Kringle等；④二級(jí)結(jié)構(gòu)，四級(jí)結(jié)構(gòu)； ⑤和其他蛋白的序列相似性(Sequence Similarity)；⑥相關(guān)疾病(Associated Diseases)和序列變異( Variants)等 (來自付強(qiáng))

KEGG數(shù)據(jù)庫

KEGG數(shù)據(jù)庫是一個(gè)有助于了解高級(jí)功能和生物學(xué)系統(tǒng)的工具。實(shí)際上，KEGG通過對(duì)生物學(xué)過程進(jìn)行計(jì)算機(jī)化處理，構(gòu)建模塊并繪制圖表，從而對(duì)基因的功能進(jìn)行系統(tǒng)化的分析。KEGG也是將基因組中的一系列基因用一個(gè)細(xì)胞內(nèi)的分子相互作用網(wǎng)絡(luò)連接起來的過程，如一個(gè)通路或是一個(gè)復(fù)合物，通過他們來展現(xiàn)更高一級(jí)的生物學(xué)功能。目的是從基因組和分子水平上了解高一級(jí)層次的功能和與之作用的生物信息資源。KEGG數(shù)據(jù)庫中還包含疾病和藥物信息等

InterProScan數(shù)據(jù)庫

通過蛋白質(zhì)結(jié)構(gòu)域和功能位點(diǎn)數(shù)據(jù)庫預(yù)測(cè)蛋白質(zhì)功能。是EBI開發(fā)的一個(gè)集成了蛋白質(zhì)家族、結(jié)構(gòu)域和功能位點(diǎn)的非冗余數(shù)據(jù)庫。Interproscan整合了一些使用最普及的一些數(shù)據(jù)庫，并應(yīng)用于功能未知的蛋白進(jìn)行Interpro注釋和GO注釋 (來自陳連福的生信博客)

Pfam

Pfam數(shù)據(jù)庫是蛋白質(zhì)家族的大集合，每個(gè)蛋白質(zhì)家族由多序列比對(duì) 和隱馬爾可夫模型（HMM）表示
蛋白質(zhì)通常由一個(gè)或多個(gè)功能區(qū)組成，通常稱為結(jié)構(gòu)域。不同的結(jié)構(gòu)域組合產(chǎn)生了天然存在的各種蛋白質(zhì)。因此，鑒定蛋白質(zhì)內(nèi)發(fā)生的結(jié)構(gòu)域可以提供對(duì)其功能的見解

方法一: 以KEGG的注釋結(jié)果為主, 篩選出每個(gè)品種包含的特異通路及基因

文件準(zhǔn)備

https://www.genome.jp/kegg/catalog/org_list.html下載所有植物通路及包含基因的文件
目的區(qū)段功能注釋文件

基本思路

將植物通路構(gòu)建了一個(gè)字典dic1，通路名稱為key，包含的基因?yàn)関alue
9個(gè)文件，每個(gè)建一個(gè)字典dic2（以一個(gè)為例），geneID為key，KEGG注釋部分為value
取dic2 value中的KEGG_ID到dic1 value中遍歷，匹配成功后，將此 dic1 value對(duì)應(yīng)的key作為最終字典dic3的key，將基因ID+‘==>’+dic2 value作為value
這個(gè)dic3包含了9個(gè)文件中含有keggid能匹配到蘋果通路中的所有基因
從字典中選擇只包含唯一一個(gè)品種的通路及其value, 認(rèn)為此value包含的基因即為該品種的特異功能基因

腳本

這個(gè)腳本輸出部分寫的異常垃圾, 當(dāng)時(shí)只是覺得自己能看懂就行了

# -*- encoding:UTF-8 -*-
import re
import os
import pprint

def text_open(file_path):
    '''用于一個(gè)通路文件和9個(gè)特異基因文件的導(dǎo)入'''
    try:
        text_name = []
        for line in open(file_path):
            line = line.strip()
            text_name.append(line)
        return text_name
    except:
        print "The 'text_open' has a question!"

def pathway_dic(pathway_text):
    '''通路文件生成字典，通路名稱為key，通路包含的基因?yàn)関alue'''
    try:
        dic = {}
        partern1 = re.compile('^C')
        partern2 = re.compile('^D')
        for line in pathway_text:
            if partern1.search(line):
                key_name =line
                dic[key_name] = []
                continue
            if partern2.search(line):
                dic[key_name].append(line)
        return dic
    except:
        print "The 'pathway_dic' has a question!"

def geneid2keggid_dic(pear_text):
    '''特異基因信息生成字典1，基因id為key，kegg id為value'''
    try:
        dic = {}
        partern1 = re.compile('K\d{5}')
        for line in pear_text:
            line_table = line.split('\t')
            if partern1.search(line_table[8]):
                dic[line_table[0]] = partern1.search(line_table[8]).group()
        return dic
    except:
        print "The 'geneid2keggid_dic' has a question!"

def geneid2keggline_dic(pear_text):
    '''特異基因信息生成字典2，基因id為key，kegg注釋信息為value'''
    try:
        dic = {}
        partern1 = re.compile('K\d{5}')
        for line in pear_text:
            line_table = line.split('\t')
            if partern1.search(line_table[8]):
                dic[line_table[0]] = line_table[8]
        return dic
    except:
        print "The 'geneid2keggline_dic' has q question!"

def nine_list(osdir):
    '''將9個(gè)文件的列表寫入一個(gè)更大列表中'''
    try:
        specific_genefile_list = []
        big_table = []
        partern1 = re.compile('^P\w+_.*')
        # partern1 = re.compile('^sun.*')
        for lists in os.listdir(osdir):
            if partern1.search(lists):
                specific_genefile_list.append(lists)
        for file in specific_genefile_list:
            big_table.append(text_open(osdir+'/'+file))
        return big_table
    except:
        print "The 'nine_list' has a question!"

def our_assume1(big_dic):
    try:
        new_dic = {}
        for big_dic_key in big_dic.keys():
            if big_dic[big_dic_key] != []:
                n = 1
                small_table = big_dic[big_dic_key][-1].split('==>')
                partern2 = re.compile('^P[a-zA-Z]+')
                word = partern2.search(small_table[0]).group()
                partern1 = re.compile(word)
                for big_dic_value in big_dic[big_dic_key]:
                    if not partern1.search(big_dic_value):
                        n +=1
                if n == 1:
                    new_dic[big_dic_key] = big_dic[big_dic_key]
        return new_dic
    except:
        print "The 'our_assume1' has a question!"

def our_assume2and3(big_dic):
    try:
        new_dic1 = {}
        new_dic2 = {}
        new_dic3 = {}
        new_dic4 = {}
        for big_dic_key in big_dic.keys():
            if big_dic[big_dic_key] != []:
                words = []
                for big_dic_value in big_dic[big_dic_key]:
                    small_table = big_dic_value.split('==>')
                    partern2 = re.compile('^P[a-zA-Z]+')
                    words.append(partern2.search(small_table[0]).group())
                if 'Pbr' in words and 'Ppc' in words and 'Puc' in words and 'Pco' in words and 'Psi' in words and  'Ppw' not in words and 'Puw' not in words and 'Ppy' not in words and 'Pbe' not in words:
                    new_dic1[big_dic_key] = big_dic[big_dic_key]
                if 'Pbr' not in words and 'Ppc' not in words and 'Puc' not in words and 'Pco' not in words and 'Psi' not in words and  'Ppw' in words and 'Puw' in words and 'Ppy' in words and 'Pbe' in words:
                    new_dic2[big_dic_key] = big_dic[big_dic_key]
                if 'Pbr' in words and 'Ppc' in words and 'Puc' in words and 'Ppw' in words and 'Psi' in words and 'Puw' in words and 'Pbe' in words and 'Pco' not in words and 'Ppy' not in words:
                    new_dic3[big_dic_key] = big_dic[big_dic_key]
                if 'Pbr'not in words and 'Ppc'not in words and 'Puc' not in words and 'Ppw'not in words and 'Psi'not in words and 'Puw'not in words and 'Pbe' not in words and 'Pco' in words and 'Ppy' in words:
                    new_dic4[big_dic_key] = big_dic[big_dic_key]
        return new_dic1, new_dic2, new_dic3, new_dic4
    except:
        print "The 'our_assume2' has a question!"


if __name__ == '__main__':
    '''problems 不能找到其他物種沒有定義的通路；沒有KEGG注釋的基因沒辦法參與定位；沒有定位到的基因可能位于非特異區(qū)域'''

    big_table = nine_list('/Users/suncongrui/home/03_workshop/specific_gene/target_gene')
    pathway_text = text_open('/Users/suncongrui/home/03_workshop/specific_gene/plant/ko00001.keg')
    pathway_dic1 = pathway_dic(pathway_text)
    big_dic = {}
    for def2_key in pathway_dic1.keys():
        big_dic[def2_key] = []
    for specise in big_table:
        geneid2keggline_dic1 = geneid2keggline_dic(specise)
        geneid2keggid_dic1 = geneid2keggid_dic(specise)
        for geneid in geneid2keggid_dic1.keys():
            partern1 = re.compile(geneid2keggid_dic1[geneid])
            for pathway_dir_key in pathway_dic1.keys():
                for pathway_dir_value in pathway_dic1[pathway_dir_key]:
                    if partern1.search(pathway_dir_value):
                        big_dic[pathway_dir_key].append(geneid+'==>'+geneid2keggline_dic1[geneid])

    sun = our_assume1(big_dic)
    cong1, cong2, rui1, rui2 = our_assume2and3(big_dic)
    pp = pprint.PrettyPrinter(indent=4)
    # pp.pprint(sun)
    # pp.pprint(cong1)
    # pp.pprint(cong2)
    # pp.pprint(rui1)
    pp.pprint(rui2)

部分結(jié)果展示

輸出異常垃圾...

垃圾輸出

此方法的問題

9個(gè)品種的很多基因沒有KEGG注釋結(jié)果
有的基因在很多通路中起作用, 但可能有主次之分, 這樣的基因也過濾掉了
沒有定位到的基因可能位于非特異區(qū)域

方法二: 利用pfam的注釋結(jié)果, 過濾掉每個(gè)品種具有相同結(jié)構(gòu)域的基因, 對(duì)于剩下的基因, 利用其swissprot id 在DAVID數(shù)據(jù)庫中查找他們的功能信息, 單個(gè)進(jìn)行篩選

文件準(zhǔn)備

目的區(qū)段功能注釋文件
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz

腳本

腳本一: 過濾掉相同結(jié)構(gòu)域的基因

# -*- encoding:UTF-8 -*-
import re
import os

def text_open(file_path):
    '''用于一個(gè)通路文件和9個(gè)特異基因文件的導(dǎo)入'''
    try:
        text_name = []
        for line in open(file_path):
            line = line.strip()
            text_name.append(line)
        return text_name
    except:
        print "The 'text_open' has a question!"

def nine_list(osdir):
    '''將9個(gè)文件的id-pfam列寫入一個(gè)更大列表中'''
    try:
        specific_genefile_list = []
        big_table = []
        partern1 = re.compile('^P\w+_.*')
        for lists in os.listdir(osdir):
            if partern1.search(lists):
                specific_genefile_list.append(lists)
        for file in specific_genefile_list:
            pfam = []
            text_name = text_open(osdir+'/'+file)
            for line in text_name:
                line_table = line.split('\t')
                pfam.append(line_table[0]+'\t'+line_table[-1])
            big_table.append(pfam)
        return big_table
    except:
        print "The 'nine_list' has a question!"

if __name__ == '__main__':
    big_table = nine_list('/Users/suncongrui/home/03_workshop/specific_gene/target_gene')
    table = []
    for i in range(9):
        big_table_copy = big_table[:]
        table1 = []
        table2 = []
        del big_table_copy[i]
        for small_table_copy in big_table_copy:
            for line in small_table_copy:
                line_table1 = line.split('\t')
                table1.append(line_table1[1])
        for line in big_table[i]:
            line_table2 = line.split('\t')
            if line_table2[1] in table1:
                table2.append(line_table2[0]+'\t'+line_table2[1]+'\t'+'NO')
            else:
                table2.append(line_table2[0]+'\t'+line_table2[1]+'\t'+'YES')
        table.append(table2)
    for b in table:
        for c in b:
            print c

腳本二: 因?yàn)閟wissprot注釋結(jié)果沒有id, 所有用下載的數(shù)據(jù)庫匹配對(duì)應(yīng)基因的id

# -*- encoding:UTF-8 -*-
import re
import os

def text_open(file_path):
    '''用于一個(gè)通路文件和9個(gè)特異基因文件的導(dǎo)入'''
    try:
        text_name = []
        for line in open(file_path):
            line = line.strip()
            text_name.append(line)
        return text_name
    except:
        print "The 'text_open' has a question!"

def nine_list(osdir):
    '''將9個(gè)文件的id-swissprot-id列寫入一個(gè)更大列表中'''
    try:
        specific_genefile_list = []
        big_table = []
        partern1 = re.compile('^P\w+_.*')
        for lists in os.listdir(osdir):
            if partern1.search(lists):
                specific_genefile_list.append(lists)
        for file in specific_genefile_list:
            text_name = text_open(osdir+'/'+file)
            for line in text_name:
                line_table = line.split('\t')
                big_table.append(line_table[0]+'\t'+line_table[7])
        return big_table
    except:
        print "The 'nine_list' has a question!"

if __name__ == '__main__':
    big_table = nine_list('/home/Dai_XG/software/user/74/target_gene')
    swissprot_text = text_open('/home/Dai_XG/software/user/74/uniprot_sprot.head')
    table = []
    for element1 in big_table:
        line_table = element1.split('\t')
        n = 0
        for element2 in swissprot_text:
            if line_table[1] == 'NA':
                continue
            if line_table[1] in element2:
                split = element2.split('|')
                print line_table[0]+'\t'+split[1]+'\t'+line_table[1]
                n += 1
        if n == 0:
            print line_table[0]+'\t'+'NO'+'\t'+line_table[1]

查找部分核心展示

這里要感謝一下chuxinxzz同學(xué)的幫助, 效率提高了不少

將運(yùn)行腳本的到的結(jié)果id寫入單獨(dú)的txt文件中
提交到DAVID數(shù)據(jù)庫中, 可以得到比較全面的功能列表
DAVID
然后挨著往后讀, 找到疑似結(jié)果就記錄一下, 這比單個(gè)找出來去每個(gè)數(shù)據(jù)庫搜索方便了不少

最后

夏天很熱 -_-

熱死我了

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

「生信」基因組功能注釋后的基因篩選

「生信」基因組功能注釋后的基因篩選

寫在前面

功能注釋數(shù)據(jù)庫介紹

GO數(shù)據(jù)庫

NCBI數(shù)據(jù)庫

Swissprot數(shù)據(jù)庫

KEGG數(shù)據(jù)庫

InterProScan數(shù)據(jù)庫

Pfam

方法一: 以KEGG的注釋結(jié)果為主, 篩選出每個(gè)品種包含的特異通路及基因

文件準(zhǔn)備

基本思路

腳本

部分結(jié)果展示

此方法的問題

方法二: 利用pfam的注釋結(jié)果, 過濾掉每個(gè)品種具有相同結(jié)構(gòu)域的基因, 對(duì)于剩下的基因, 利用其swissprot id 在DAVID數(shù)據(jù)庫中查找他們的功能信息, 單個(gè)進(jìn)行篩選

文件準(zhǔn)備

腳本

查找部分核心展示

最后

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

「生信」基因組功能注釋后的基因篩選

寫在前面

功能注釋數(shù)據(jù)庫介紹

GO數(shù)據(jù)庫

NCBI數(shù)據(jù)庫

Swissprot數(shù)據(jù)庫

KEGG數(shù)據(jù)庫

InterProScan數(shù)據(jù)庫

Pfam

方法一: 以KEGG的注釋結(jié)果為主, 篩選出每個(gè)品種包含的特異通路及基因

文件準(zhǔn)備

基本思路

腳本

部分結(jié)果展示

此方法的問題

方法二: 利用pfam的注釋結(jié)果, 過濾掉每個(gè)品種具有相同結(jié)構(gòu)域的基因, 對(duì)于剩下的基因, 利用其swissprot id 在DAVID數(shù)據(jù)庫中查找他們的功能信息, 單個(gè)進(jìn)行篩選

文件準(zhǔn)備

腳本

查找部分核心展示

最后

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av