DNA to Protein seq

DNA seq to Protein seq

from Bio.Seq import Seq
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            length = 3 * ((len(seq)-frame) // 3)  # Adjust length for frame
            for pro in nuc[frame:frame+length].translate(table=trans_table).split("*"):
                if len(pro) >= min_protein_length:
                    answer.append(pro)
    return answer
    
# FASTA 파일에서 서열 읽기 및 ORF 번역
all_proteins = []
with open("./fasta/Hymenobacter rubidus.fasta", "r") as file:
    for record in SeqIO.parse(file, "fasta"):
        orfs = find_orfs_with_trans(record.seq, "Standard", min_protein_length=100)
        all_proteins.extend(orfs)
        
# 결과 저장
with open("protein/Hymenobacter_rubius_translated_proteins.fasta", "w") as output_file:
    for i, protein in enumerate(all_proteins):
        record = SeqRecord(Seq(str(protein)), id=f"protein_{i}", description="")
        SeqIO.write(record, output_file, "fasta")

 

 

번역된 단백질 서열에서 recA 서열 찾기

from Bio import pairwise2
from Bio.pairwise2 import format_alignment

# RecA 단백질의 아미노산 서열
recA_sequence = "MSKDATKEISAPTDAKERSKAIETAMSQIEKAFGKGSIMKLGAESKLDVQVVSTGSLSLDLALGVGGIPRGRITEIYGPESGGKTTLALAIVAQAQKAGGTCAFIDAEHALDPVYARALGVNTDELLVSQPDNGEQALEIMELLVRSGAIDVVVVDSVAALTPRAEIEGDMGDSLPGLQARLMSQALRKLTAILSKTGTAAIFINQVREKIGVMYGNPETTTGGRALKFYASVRLDVRKIGQPTKVGNDAVANTVKIKTVKNKVAAPFKEVELALVYGKGFDQLSDLVGLAADMDIIKKAGSFYSYGDERIGQGKEKTIAYIAERPEMEQEIRDRVMAAIRAGNAGEAPALAPAPAAPEAAEA"  # RecA 아미노산 서열 (실제 서열로 교체 필요)

# 유사도 임계값 설정
similarity_threshold = 0.8  # 예: 80% 유사도


fasta_file_path = "./protein/protein.fasta"  # fasta 파일의 경로를 지정하세요
translated_proteins = []

with open(fasta_file_path, "r") as fasta_file:
    current_sequence = ""
    for line in fasta_file:
        line = line.strip()  # 줄 끝의 공백 및 개행 문자를 제거합니다.
        if line.startswith(">"):
            # 헤더 줄이면 현재 서열을 리스트에 추가하고 새로운 서열 준비
            if current_sequence:
                translated_proteins.append(current_sequence)
            current_sequence = ""
        else:
            # 서열 데이터를 현재 서열에 추가
            current_sequence += line

# 마지막 서열을 추가 (파일의 마지막 줄에 헤더가 없을 경우)
if current_sequence:
    translated_proteins.append(current_sequence)

# translated_proteins 리스트에는 모든 번역된 단백질 서열이 포함됨.


# RecA와 유사한 서열 찾기
similar_sequences = []
for protein in translated_proteins:
    # 최적의 정렬과 점수 계산
    alignments = pairwise2.align.globalxx(recA_sequence, protein)
    alignment_result = pairwise2.format_alignment(*alignments[0])
    score_line = alignment_result.split('\n')[-2]  # 마지막 줄에서 score를 찾습니다.
    score = score_line.split("=")[1]  # "=" 문자를 기준으로 score를 추출합니다.
    score = score.strip()  # 문자열 양쪽의 공백을 제거합니다.

    
    # 유사도 계산
    similarity = float(score) / len(recA_sequence)
    
    # 유사도 임계값 이상인 경우 저장
    if similarity >= similarity_threshold:
        similar_sequences.append((protein, similarity))

# 유사한 서열 출력
print("Sequences similar to RecA:")
for seq, sim in similar_sequences:
    print(f"Sequence: {seq}")
    print(f"Similarity: {sim}")

• 출력된 recA 유사 서열 중에 similarity가 가장 높은 서열을 recA 서열로 채택
• 사용된 recA 서열은 BAA21330.1 RecA [Deinococcus radiodurans]의 서열

BlastP (Protein Blast) in linux

# Blast 설치
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.15.0+-x64-linux.tar.gz
tar -xvzf ncbi-blast-2.15.0+-x64-linux.tar.gz
export PATH=$PATH:/home/seona/DNAtoPROTEIN/ncbi-blast-2.15.0+/bin
source ~/.bashrc

# blastp로 두 개의 단백질 서열 비교
blastp -query ./protein/M29_recA.fasta -subject ./protein/recA.fasta