基因预测与比对

常用的原核生物基因预测软件有GeneMarks、Glimmer3、RAST

使用GeneMarks 进行原核生物基因的预测

1. GeneMarks 的下载和安装

wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_aQKBZ/genemark_suite_linux_64.tar.gz  
tar -zxvf genemark_suite_linux_64.tar.gz -C ~/tools/  
##添加环境变量  
vim ~/.bashrc  
export PATH=/home/wei/tools/gmsuite/:$PATH  
source ~/.bashrc  
## 添加密钥  
gzip -dv gm_key_64.gz  
mv gm_key_64 ~/.gm_key

2. GeneMarks 的使用

gmsn.pl --prok --format GFF --pdf --fnn --faa YP14.assembly.fasta

--output <string> 输出基因预测的结果文件，默认为输入的文件名
--prok 对原核生物序列进行基因预测
--format 输出的基因预测结果的格式文件，默认为LTS，还可以GFF，GFF3
--pdf 生成pdf 的图形结果
--fnn 生成预测基因的核酸序列
--faa 生成预测基因的氨基酸序列

输出的文件有.faa .fnn .gff .pdf 等文件，以方便下一步的比对工作

使用 Glimmer3进行原核生物基因的预测

Glimmer是用于寻找微生物DNA，特别是细菌、古菌和病毒中的基因。
其采用的方法为内插马尔科夫模型（interpolted Markov model，IMM）
来识别编码区域和非编码区域

1. glimmer3 的下载与安装

wget http://ccb.jhu.edu/software/glimmer/glimmer302b.tar.gz  
tar -zxvf glimmer302b.tar.gz -C ~/tools/  
cd ~/tools/glimmer3.02/src  
sudo make  
vim ~/.bashrc
export PATH=/home/wei/tools/glimmer3.02/bin/:$PATH
source ~/.bashrc

2. 使用方法

创建训练模型：
a. 用亲缘关系很近的物种的基因；
b. 用自身序列创建的orf数据；
c. 用基因组本身的已知信息；
下面用的自身数据作为训练集：

产生长orf 数据

long-orfs -n -t 1.15 YP14.seq run.longorfs  
-n 输出文件去除首行，只包含orf  
-t 熵距离得分阈值，小于阈值才被保留

提取数据集
extract -t YP14.seq run.longorfs > run.train
生成预测模型
build-icm -r run.icm < run.train

基因预测

glimmer3 -o50 -g110 -t30 YP14.seq run.icm run  
-o 最大重叠片段长度阈值，小于阈值保留  
-g 基因片段长度阈值，大于阈值则保留  
-t orf得分阈值，大于阈值保留 

根据预测结果提取序列
extract -t YP14.seq run.predict > predict.fasta

使用ncbi-blast+ 来进行序列的比对

1. ncbi-blast+ 的下载安装

下载地址: ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.8.0alpha/

2. 建立库文件

下载Swiss-Prot 的蛋白序列并构建Blast 数据库：
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz gzip -d uniprot_sprot.fasta.gz

建立数据库：
makeblastdb -in uniprot.fasta -dbtype prot -parse_seqids -hash_index -out uniprot

-in 输入数据库文件
-dbtype 蛋白库用prot,核酸用nucl
-out 所建数据库的名称
-parse_seqids      
-hash_index 创建哈希序列

将数据库的文件加入到特定位置

vim ~/.ncbirc  
[BLAST]
#BLAST数据库的存放位置
BLASTDB=/home/wei/database/
#BLASTDB_PROT_DATA_LOADER=//home/wei/database
#BLASTDB_NUCL_DATA_LOADER=//home/wei/database  

提供的脚本下载数据库

#展示可以下载的所有数据库  
update_blastdb.pl --showall    
#自动在后台下载数据库，然后自动解压  
nohup perl update_blastdb.pl --decompress nr >out.log 2>&1 &

查看数据库信息

blastdbcmd -db 16SMicrobial -dbtype nucl -info  
blastdbcmd -db 16SMicrobial -dbtype nucl -entry all | head  
blastdbcmd -db 16SMicrobial -dbtype nucl -entry NR_118889.1 | head

3. 数据的比对

blastp -query pritein.fasta -db uniprot -evalue 1e-5 -out blast.xml -outfmt 5 -num_alignments 10 -num_threads 4 [-remote]
参数说明：

-query 输入要比对的文件路径
-db 格式化后的数据库名称
-evalue 设定输出结果中的e-value阈值
-out 输出文件名
-num_alignments 输出比对上的序列的最大值条目数
-num_threads 线程数
-num_descriptions 序列描述信息，一般跟tabular格式连用
-remote 在联网状态下会利用NCBI的服务器来进行比对
-outfmt      
 0 = pairwise,
 1 = query-anchored showing identities,
 2 = query-anchored no identities,
 3 = flat query-anchored, show identities,
 4 = flat query-anchored, no identities,
 5 = XML Blast output,输出的信息最全
 6 = tabular,表格的格式
 7 = tabular with comment lines,
 8 = Text ASN.1,
 9 = Binary ASN.1
 10 = Comma-separated values

使用diamond 来进行快速比对

1. diamond 的基本使用

diamond 的比对速度要比blast快的多

wget http://github.com/bbuchfink/diamond/releases/download/\  
v0.9.22/diamond-linux64.tar.gz  
tar xzf diamond-linux64.tar.gz
#使用方法与blast类似，但只能比对蛋白数据  
#加入到环境变量  
vim ~/.bashrc
export PATH=/home/wei/tools/diamond:$PATH  
diamond help #查看参数信息  
diamond makedb --in uniprot_sprot.fasta --db uniprot  
diamond blastp --query YP14.faa --db uniprot --out   \
11.xml --outfmt 5  --evalue 1e-5     

比对信息的提取脚本

比对出来的信息可能需要进一步提取，如果需要具体的描述信息和identity信息
我们需要将格式5 与格式6 的数据综合一下
shell 脚本为：

#一共设置了4个命令行参数    
$1 -query $2 -db $3 -out $4.xml -outfmt 5 -evalue 1e-5 -num_alignments 1 -num_threads 4;
$1 -query $2 -db $3 -out $4.tab -outfmt 6 -evalue 1e-5 -num_alignments 1 -num_threads 4;  

下面是提取数据用到的python 脚本

"""
将数据存储为文本文件格式
"""
from lxml import etree  
import sys

a = []
#第一个命令行参数为blast 6 输出的文件
f = open(sys.argv[1])
for i in f:
    j = i.split("\t")[2]
    a.append(j) 

all_info_list = []
header = ['Iteration_query_def','Hit_def','Hit_accession','Hsp_evalue','Hsp_score','identity']
all_info_list.append(header)
def get_info(html):  
   #先锁定每个节点
    infos = html.xpath('//Iteration')
    m = 0
    for info in infos:
        #过滤没有匹配到的hits
        try:
            Hit_def = info.xpath('Iteration_hits/Hit/Hit_def/text()')[0]
        except IndexError:
            #pass continue and break
            continue
        Iteration_query_def = info.xpath('Iteration_query-def/text()')[0]
        Hit_accession = info.xpath('Iteration_hits/Hit/Hit_accession/text()')[0]
        Hsp_evalue = info.xpath('Iteration_hits/Hit/Hit_hsps/Hsp/Hsp_evalue/text()')[2]
        Hsp_score = info.xpath('Iteration_hits/Hit/Hit_hsps/Hsp/Hsp_score/text()')[0]
        info_list = [Iteration_query_def,Hit_def,Hit_accession,Hsp_evalue,Hsp_score,a[m]]
        all_info_list.append(info_list)
        m += 1
    print (all_info_list)
if __name__ == '__main__':
    #第二个命令行参数为 blast 5输出的文件
    html = etree.parse(sys.argv[2])
    get_info(html)
    #第三个命令行参数为合并后的文件输出
    f = open(sys.argv[3],'w')
    for i in all_info_list:
        for j in i:
            f.write(j) 
            f.write("\t")
        f.write("\n")
    f.close()
"""
将数据存储为excel 表格形式
"""
from lxml import etree 
import xlwt
import sys
a = []
f = open(sys.argv[1])
for i in f:
    j = i.split("\t")[2]
    a.append(j)    

all_info_list = []
def get_info(html):
    #先锁定每个节点
    infos = html.xpath('//Iteration')  
    m = 0 
    for info in infos:
        #过滤没有匹配到的hits
        try:
            Hit_def = info.xpath('Iteration_hits/Hit/Hit_def/text()')[0]
        except IndexError:
            #pass continue and break
            continue
        Hit_accession = info.xpath('Iteration_hits/Hit/Hit_accession/text()')[0]
        Hsp_evalue = info.xpath('Iteration_hits/Hit/Hit_hsps/Hsp/Hsp_evalue/text()')[0]
        Iteration_query_def = info.xpath('Iteration_query-def/text()')[0]
        Hsp_score = info.xpath('Iteration_hits/Hit/Hit_hsps/Hsp/Hsp_score/text()')[0]
        info_list = [Iteration_query_def,Hit_def,Hit_accession,Hsp_evalue,Hsp_score,a[m]]
        all_info_list.append(info_list)   
        m += 1
    print (all_info_list)
if __name__ == '__main__':
    html = etree.parse(sys.argv[2])
    get_info(html)
    header = ['Iteration_query_def','Hit_def','Hit_accession','Hsp_evalue','Hsp_score','identity']
    book = xlwt.Workbook(encoding="utf-8")
    sheet = book.add_sheet('Sheet1')
    for h in range(len(header)):
        sheet.write(0,h,header[h])#写入表头
    i = 1
    for list in all_info_list:
        j = 0
        for data in list:
            sheet.write(i,j,data)
            j+= 1
        i += 1
book.save(sys.argv[3])

除此之外实验室进行的大量分菌测序工作需要16s比对
拿到的结果为上海生工的不标准序列文件，合并脚本如下面：

import sys
import os
#遍历指定的目录，显示目录下的所有文件名
seq_list = []
file_list = []
def eachFile(filepath):
    pathDir = os.listdir(filepath)
    for allDir in pathDir:
        child = os.path.join('%s%s' %(filepath,allDir))
        file_list.append(child)
    return file_list
def readFile(filename):
    line_num = 0
    fopen = open(filename,'r')
    for line in fopen:
        seqname = filename.split('/')[-1]
        if '^^' in line:
            line = line.replace('^^','>' + seqname.replace('拼接结果',''))
        line_num += 1
	#跳过第一行打印输出
        if (line_num != 1):
            seq_list.append(line.strip())
            print (line.strip())
    fopen.close()     
if __name__ == '__main__':  
    #文件的打开路径
    for i in eachFile(sys.argv[1]):
        readFile(i)  
    #设置文件的保存路径
    f = open(sys.argv[2],'w')
    for i in seq_list:
        f.write(i + '\n')
    f.close()

COG 的数据库注释需要写脚本来解决

from lxml import etree 
import sys
"""
将whog 文件的数据提取出来存储为字典格式
"""
f = open("sys.argv[1]")
dict = {}
for line in f:
    if line[0] == "\n":
        continue
    if line[0] == "_":
        continue
    if line[0] == "[":
        key = line.strip()
        dict[key] = []
    else:
        dict[key].append(line.strip())
        
"""
blast输出 6 格式的文件identity的提取
"""
a = []
f = open(sys.argv[2])
for i in f:
    j = i.split("\t")[2]
    a.append(j)
    
"""
blast输出 5 格式的文件详细信息的提取
"""    
all_info_list = []
header = ['Iteration_query_def','Hit_def','Hit_accession','Hsp_evalue',
          'Hsp_score','identity','Function_class','COG_num','Functional_description']
all_info_list.append(header)
def get_info(html):  
   #先锁定每个节点
    infos = html.xpath('//Iteration')
    m = 0
    for info in infos:
        #过滤没有匹配到的hits
        try:
            Hit_def = info.xpath('Iteration_hits/Hit/Hit_def/text()')[0]
        except IndexError:
            #pass continue and break
            continue
        Iteration_query_def = info.xpath('Iteration_query-def/text()')[0]
        Hit_accession = info.xpath('Iteration_hits/Hit/Hit_accession/text()')[0]
        for key,values in dict.items():
            for i in values:
                if  Hit_accession in i:
                    n= key.split(" ") 
                    Function_class = n[0].replace("[",'').replace("]",'')
                    COG_num = n[1]
                    Functional_description = ' '.join(n[2:])
        Hsp_evalue = info.xpath('Iteration_hits/Hit/Hit_hsps/Hsp/Hsp_evalue/text()')[0]
        Hsp_score = info.xpath('Iteration_hits/Hit/Hit_hsps/Hsp/Hsp_score/text()')[0]       
        info_list = [Iteration_query_def,Hit_def,Hit_accession,Hsp_evalue,Hsp_score,a[m],
                     Function_class,COG_num,Functional_description]
        all_info_list.append(info_list)
        m += 1
    print (all_info_list)
if __name__ == '__main__':
    #第二个命令行参数为 blast 5输出的文件
    html = etree.parse(sys.argv[3])
    get_info(html)
    #第三个命令行参数为合并后的文件输出
    f = open(sys.argv[4],'w')
    for i in all_info_list:
        for j in i:
            f.write(j) 
            f.write("\t")
        f.write("\n")
    f.close()

GeneMarks、glimmer3 与blast+ 的使用以及比对信息的提取