한타바이러스의 시퀀스를 받아서 MSA를 해보자 (후편)

https://koreanraichu.tistory.com/824

한타바이러스의 시퀀스를 받아서 MSA를 해보자 (전편)

우리 여기서 MSA 돌릴 FASTA파일 저장했죠? 후편에서는 그걸로 분석을 돌릴건데… 아… 이거 트리땜에 개노가다했음… 그것도 트리 얘기하면서 얘기해드리겠음.

---

# 바이러스 DNA가 세그먼트별로 섞여있습니다. (S, M, L)
# 이거 분류 안하면 MSA 뻑나요. 
s_segments = [rec for rec in vir_sequence if 1600 <= len(rec.seq) <= 2000]
SeqIO.write(s_segments, "hantavirus_segmemt_s.fasta", "fasta")

# 세그먼트 ID를 이름으로 바꾸는 절차라고 보시면 됩니다. 
for rec in s_segments:
    rec.id = f"{rec.id}_{len(rec.seq)}"

print("Completed.")

바이러스 DNA를 세그먼트별로 나눠야 한다. 왜죠? 하나의 바이러스가 여러 개의 파일(segment)로 나뉘어 있는 유전체니까요. 저게 스몰 미디엄 라지가 길이별로 다 다른데 안 나누고 묶어서 MSA 돌리잖아요? 그러면 컴퓨터가 장난 똥때리나 마!!! 하면서 갭을 긴 놈한테 맞춰서 결과가 개발살납니다.

 

print('MSA start... ')

# MSA 분석 시-작
try: 
    result = subprocess.run([muscle_exe, "-align", "hantavirus_segmemt_s.fasta", "-output", "Hantavirus_alignment_s.afa"], check=True, capture_output=True, text=True)
    print("Completed. ")
except subprocess.CalledProcessError as e: 
    print(f"MSA failed: {e}")

2차 똥타임을 가질 수 있는 절호의 찬스. 이게 개수가 많아지면 MSA 하는데도 시간이 엄청 걸립니다. 얘는 좀 덜하긴 함… 인플루엔자는 하나 돌리는데 2~3분 걸리더라… 100갠가 갖고왔는데…

 

# 다 됐으면 몇 개만 확인해보자. 
print("====== MSA Result ======")
alignment = AlignIO.read("Hantavirus_alignment_s.afa", "fasta") # FASTA 니네 확장자가 몇개냐... 

for record in alignment:
    print(f"{record.id[:15]:<15} : {record.seq[:50]}")

결과를 확인할 수 있는 코드이다. 보고 싶은 시퀀스 길이를 조절하고 싶으시다고요? record.seq[:50]를 조절하십쇼.

 

# 코어 시퀀스는 어디? (바이러스라고 앞뒤 안가리고 다 변형하는거 아님)
def calculate_conservation(alignment):
    length = alignment.get_alignment_length()
    scores = []
    for i in range(length):
        column = alignment[:, i]
        most_common = max(column, key=column.count)
        score = column.count(most_common) / len(column)
        scores.append(score)
    return scores

scores = calculate_conservation(alignment)

print(f"해당 구간의 평균 보존율: {np.mean(scores)*100:.2f}%")

# 변이 핫스팟 찾기 (엔트로피 분석)
variation_scores = [1 - s for s in scores] # 위에서 계산한 그거 

# 핫스팟 시각화
plt.figure(figsize=(15, 5))
plt.plot(variation_scores, color='red', alpha=0.7)
plt.fill_between(range(len(variation_scores)), variation_scores, color='red', alpha=0.2)

plt.title("Viral Variation Hotspots (Potential Immune Evasion Sites)")
plt.xlabel("Sequence Position")
plt.ylabel("Mutation Frequency (Entropy)")
plt.show()

이건 또 뭐임? 바이러스가 변이율이 개삽쩌는건 다들 아시죠? 근데 얘네들도 앞뒤 안가리고 일단 변이부터 하고 보지는 않는다. 덮어놓고 변이했다가 게놈에 문제생겨서 증식을 못하거나 껍질이 제대로 안 만들어지거나 하면 바이러스 입장에서는 ㄹㅇ 초유의 사태이기 떄문이다. 그래서 얘들도 건드려도 되는 부위와 건드리면 안 되는 부위가 있거든요? 거기를 엔트로피 분석을 통해서 보는거다.

 

...빗이여? 근데 엔트로피는 열역학 아니예요? 그게 왜 갑자기 여기서 나옴? 뭐 바이러스가 변이할때마다 엔트로피 증가함? 아니 그거 아님. 섀넌 엔트로피, 혹은 정보 엔트로피라고 하는건데 확률변수의 불확실성(uncertainty)을 정량화하는 척도를 의미한다. 근데 엔트로피랑 저거랑 뭔 상관이냐고? 변이가 많다는 건 염기가 다르다는거잖아요. 그럼 저기에 A가 오는지 T가 오는지 G가 오는지 C가 오는지 몰라요. 이거는 슈뢰딩거도 박스 까봐야 알아. 그럼 막대기가 낮은 부분은요? 상대적으로 변이가 덜 되는거다. 어쩌면 막대가 상대적으로 낮은 부분이 변이하면 바이러스 입장에서 개망하는 부분일수도 있지... 

 

# Phylogenic tree

# 1. 거리 계산
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(alignment)

# 2. Phylogenic tree 생성(NJ)
constructor = DistanceTreeConstructor(calculator, 'nj')
tree = constructor.build_tree(alignment)

# 이름이 모예요 
clean_name_dict = {}
for record in vir_sequence:
    # "Hantaan virus" 같은 종 이름만 추출
    species_name = record.description.split("segment")[0].split(",")[0].strip()
    # 너무 긴 학명 단어는 제거 (가독성 향상)
    species_name = species_name.replace("orthohantavirus", "").strip()
    clean_name_dict[record.id] = species_name

for leaf in tree.get_terminals():
    found = False
    for original_id, clean_name in clean_name_dict.items():
        if original_id.split('.')[0] in leaf.name.replace('_', '.'):
            leaf.name = str(clean_name) # 확실하게 문자열로 변환해서 주입
            found = True
            break
    if not found:
        leaf.name = leaf.name.split('_')[0]

# Inner 노드 소멸
for clade in tree.find_clades():
    if not clade.is_terminal():
        clade.name = None

# Phylogenic tree 시각화
fig = plt.figure(figsize=(26, 25), dpi=100)
axes = fig.add_subplot(1, 1, 1)
Phylo.draw(tree, axes=axes, do_show=False, show_confidence=False, label_func=lambda x: str(x.name) if x.is_terminal() else "")

# 오케이 트리 봐봐 
plt.rc('font', size=10) # 내부 글꼴 사이즈
plt.rc('axes', titlesize=20) # 제모옥은 이 크기로 하겠습니다 

plt.title("Hantavirus S-Segment Phylogenetic Tree") # 근데 이제 이걸 곁들인
plt.xlabel("Genetic Distance (Substitution per site)", fontsize=12)
plt.ylabel("Viral Strains", fontsize=12)
plt.tight_layout()
plt.savefig("Hantavirus_Final_Tree.png", dpi=300, bbox_inches='tight')
plt.show()

이게 왜 시간을 이렇게 잡아먹었냐... 일단 처음에 그려진 트리가 어떻게 생겼는지 보여드림.

 

이게 Phylogenic tree다. 근데 안이쁘죠? 저기 이너 어쩌고 보이는것도 그렇고, 폰트는 무슨 현미경 써야 보일 것 같고, 결정적으로 그래서 쟤가 무슨 바이러스인지 모르잖음. 저거 저렇게 주면 연구자들도 NCBI 가서 저놈이 뭔 바이러스인지 찾아야 한다 이거거든요. inner는 바로 빠졌는데 트리 크기하고 이름 출력하는것때문에 ㄹㅇ 애먹었음.

 

이름 붙여놨더니 2차로 노안 간접체험 하게 생겼네… 근데 이거는 글꼴 크기 조정하면 해결은 됩니다. 아무튼 이게 NCBI에서 한타바이러스 게놈 긁어와서 S segment만 걸러서 Phylogenic tree 그린 결과임. 저기서 뭉쳐있는 놈들끼리는 친척인거고 거리가 좀 떨어져 있으면 사돈의 8촌정도 되는? 그런 애들이라고 생각하시면 됩니다. hierarchical clustering 하면 저게 90도 꺾인 그림이 나오는데 그게 덴드로그램임. 근데 좀 아쉬운건 이름 출력이 저게 최선이라 저게 그래서 한탄이여 스울이여 푸말라여는 안나옴. 쟤네 진짜 쪽수 많습니다.