目录

使用说明

完整julia程序

使用说明

本程序为julia程序，用于计算蛋白序列间的一致性（identity）。

运行本程序前请确保已经安装好了julia，如未安装，前往：

https://julialang.org/downloads

确保已安装运行本程序所需要的包，如未安装，则再命令行下输入：

julia

进入julia环境，接着依次输入：

import Pkg
Pkg.add("BioSequences")
Pkg.add("BioAlignments")
Pkg.add("ProgressMeter")
Pkg.add("Bio")

等待安装完成。

运行本程序前需要将待计算一致性（identity）的所有蛋白序列放在同一fasta文件中。

用文本编辑器打开protein_aligment.jl程序，修改倒数第3行的fasta文件路径。

重新打开命令行，进入protein_aligment.jl所在目录，输入：

julia protein_alignment.jl

等待一致性（identity）计算完成，生成一致性（identity）矩阵alignmentscores。

可以使用python直接读取该矩阵，如：

import numpy as np

alignment_scores = np.loadtxt(‘alignmentscores‘)

完整julia程序

# protein_alignment.jl

using BioSequences
using BioAlignments
using DelimitedFiles
using LinearAlgebra
using ProgressMeter

function file_to_array(file)
    sequences = []
    reader = FASTA.Reader(open(file, "r"))
    for record in reader
        seq = FASTA.sequence(record)
        push!(sequences, seq)
    end
    return sequences
end


function calculate_perc_identity(sequence1, sequence2)
    scoremodel = AffineGapScoreModel(BLOSUM62, gap_open=-10, gap_extend=-1)
    res = pairalign(LocalAlignment(), sequence1, sequence2, scoremodel);
    aln = alignment(res)

    return count_matches(aln) / min(length(sequence1), length(sequence2))
end


function calculate_score_matrix(file)
    sequence_list = file_to_array(file)
    kernel_matrix = zeros(length(sequence_list), length(sequence_list))
    p = Progress(Int64(round((length(sequence_list)^2)/2, digits=0)))

    for i in 1:length(sequence_list), j in i:length(sequence_list)
        kernel_matrix[i,j] = calculate_perc_identity(sequence_list[i], sequence_list[j])
        next!(p)
    end

    kernel_matrix = kernel_matrix + kernel_matrix‘
    kernel_matrix = kernel_matrix - Diagonal(kernel_matrix)/2
    return kernel_matrix
end


seqfile = "protein_sequences.fasta"  # fasta文件路径
m = calculate_score_matrix(seqfile)
writedlm("alignmentscores.txt", m)  # 输出路径

Julia实现基于双序列比对（Pairwise sequence Alignment）计算蛋白序列间一致性（identity）