Use the JC69 model to calculate the substitution rate at all first, second and third codon positions in the two aligned sequences in SSA3. Assume the first position is 1 so the codon positions will be: 123123123123 etc in the sequences. #Here is the code for reading in Seub_SSA3 and Spar_SSA3 with open("SSA3.fasta") as f: seqA = "" seqB = "" current_seq = "" seq_header = "" for line in f: line = line.strip() if line.startswith(">"): if seq_header == ">Seub_SSA3": seqA = current_seq elif seq_header == ">Spar_SSA3": seqB = current_seq current_seq = "" seq_header = line else: current_seq += line if seq_header == ">Seub_SSA3": seqA = current_seq elif seq_header == ">Spar_SSA3": seqB = current_seq with open("SeqA.fasta", "w") as f: f.write(seqA) with open("SeqB.fasta", "w") as f: f.write(seqB) print("Length of seqA:", len(seqA)) print("Lenght of seqB :", len(seqB)) # JC69 model def jc69(seq1, seq2): # Convert sequences to uppercase seq1 = seq1.upper() seq2 = seq2.upper() # Check that sequences have the same length if len(seq1) != len(seq2): raise ValueError("Sequences have different lengths") # Calculate number of differences diff_count = sum([seq1[i] != seq2[i] for i in range(len(seq1))]) # Calculate proportion of differences p = diff_count / len(seq1) # Calculate substitution rate K = -(3/4)*np.log(1 - p*(4/3)) return K
Use the JC69 model to calculate the substitution rate at all first, second and third codon positions in the two aligned sequences in SSA3. Assume the first position is 1 so the codon positions will be: 123123123123 etc in the sequences.
#Here is the code for reading in Seub_SSA3 and Spar_SSA3
with open("SSA3.fasta") as f:
seqA = ""
seqB = ""
current_seq = ""
seq_header = ""
for line in f:
line = line.strip()
if line.startswith(">"):
if seq_header == ">Seub_SSA3":
seqA = current_seq
elif seq_header == ">Spar_SSA3":
seqB = current_seq
current_seq = ""
seq_header = line
else:
current_seq += line
if seq_header == ">Seub_SSA3":
seqA = current_seq
elif seq_header == ">Spar_SSA3":
seqB = current_seq
with open("SeqA.fasta", "w") as f:
f.write(seqA)
with open("SeqB.fasta", "w") as f:
f.write(seqB)
print("Length of seqA:", len(seqA))
print("Lenght of seqB :", len(seqB))
# JC69 model
def jc69(seq1, seq2):
# Convert sequences to uppercase
seq1 = seq1.upper()
seq2 = seq2.upper()
# Check that sequences have the same length
if len(seq1) != len(seq2):
raise ValueError("Sequences have different lengths")
# Calculate number of differences
diff_count = sum([seq1[i] != seq2[i] for i in range(len(seq1))])
# Calculate proportion of differences
p = diff_count / len(seq1)
# Calculate substitution rate
K = -(3/4)*np.log(1 - p*(4/3))
return K
Trending now
This is a popular solution!
Step by step
Solved in 2 steps
