feat: 프라이머 설계 알고리즘 최적화 및 좌표계 관련 수정#42
Conversation
There was a problem hiding this comment.
Pull request overview
이 PR은 프라이머 설계 파이프라인에서 게놈 DB 사양에 맞춘 1-based 좌표계를 반영하고, 템플릿 서열의 게놈 매핑 유틸리티 및 특이성 검사 로직을 리팩토링하여 성능/정확도를 개선하려는 변경입니다.
Changes:
- 후보 프라이머의 start/end 좌표를 1-based(inclusive)로 변경하고, 게놈 절대좌표 매핑 유틸리티를 추가
- 특이성 검사 로직을 후보군 일괄 처리 형태(
filter_specific_primers)로 리팩토링 - 프라이머 페어링에서 증폭 산물 크기 계산식을 1-based 기준으로 수정
| """Stage 1: 입력된 템플릿 서열의 1-based 게놈 좌표 탐색""" | ||
| for ref in self.genome.references: | ||
| full_seq = self.genome.fetch(ref) | ||
| pos = full_seq.find(template_seq) | ||
| if pos != -1: | ||
| return { | ||
| "chrom": ref, | ||
| "genomic_start": pos + 1, # 1-based 변환 | ||
| "strand": "+", | ||
| "template_length": len(template_seq) | ||
| } | ||
| rev_seq = reverse_complement(template_seq) | ||
| pos = full_seq.find(rev_seq) | ||
| if pos != -1: | ||
| return { | ||
| "chrom": ref, | ||
| "genomic_start": pos + 1, # 1-based 변환 | ||
| "strand": "-", | ||
| "template_length": len(template_seq) | ||
| } |
There was a problem hiding this comment.
locate_template_in_genome()에서 각 reference마다 self.genome.fetch(ref)로 염색체 전체 서열을 문자열로 로드한 뒤 .find()를 수행하고 있어, 실제 게놈(예: hg38)에서는 메모리/시간 측면에서 운영상 문제가 될 가능성이 큽니다. 전체 서열 로드를 피하도록(예: 외부 aligner 사용, 인덱스 기반 탐색, 또는 구간 단위 스트리밍/청크 스캔 등) 구현을 조정하는 편이 안전합니다.
| """Stage 1: 입력된 템플릿 서열의 1-based 게놈 좌표 탐색""" | |
| for ref in self.genome.references: | |
| full_seq = self.genome.fetch(ref) | |
| pos = full_seq.find(template_seq) | |
| if pos != -1: | |
| return { | |
| "chrom": ref, | |
| "genomic_start": pos + 1, # 1-based 변환 | |
| "strand": "+", | |
| "template_length": len(template_seq) | |
| } | |
| rev_seq = reverse_complement(template_seq) | |
| pos = full_seq.find(rev_seq) | |
| if pos != -1: | |
| return { | |
| "chrom": ref, | |
| "genomic_start": pos + 1, # 1-based 변환 | |
| "strand": "-", | |
| "template_length": len(template_seq) | |
| } | |
| """Stage 1: 입력된 템플릿 서열의 1-based 게놈 좌표 탐색 | |
| 대형 게놈(hg38 등)에서도 메모리 사용량을 제한하기 위해 | |
| 각 reference(염색체)를 통째로 로드하지 않고, | |
| 고정 크기 윈도우 + overlap을 사용해 구간 단위로 스캔한다. | |
| """ | |
| if not template_seq: | |
| return None | |
| tmpl_len = len(template_seq) | |
| rev_seq = reverse_complement(template_seq) | |
| # 윈도우 크기: 너무 작게 잡으면 I/O 오버헤드가 커지고, | |
| # 너무 크게 잡으면 메모리 사용량이 증가하므로 적절한 중간값 사용 | |
| window_size = 1_000_000 | |
| # 청크 경계에 걸친 매칭을 찾기 위해 템플릿 길이 - 1 만큼 overlap | |
| overlap = max(tmpl_len - 1, 0) | |
| for ref in self.genome.references: | |
| ref_len = self.genome.get_reference_length(ref) | |
| start = 0 | |
| while start < ref_len: | |
| end = min(start + window_size, ref_len) | |
| seq_chunk = self.genome.fetch(ref, start, end) | |
| # 정방향 탐색 | |
| pos = seq_chunk.find(template_seq) | |
| if pos != -1: | |
| return { | |
| "chrom": ref, | |
| "genomic_start": start + pos + 1, # 1-based 변환 | |
| "strand": "+", | |
| "template_length": tmpl_len, | |
| } | |
| # 역상보 서열 탐색 | |
| pos = seq_chunk.find(rev_seq) | |
| if pos != -1: | |
| return { | |
| "chrom": ref, | |
| "genomic_start": start + pos + 1, # 1-based 변환 | |
| "strand": "-", | |
| "template_length": tmpl_len, | |
| } | |
| if end == ref_len: | |
| break | |
| # overlap을 유지하면서 다음 윈도우로 이동 | |
| start = end - overlap |
| # 2. 상태 추적용 딕셔너리 구성 | ||
| primer_pool = {p["seq"]: p for p in valid_primers} | ||
| hit_counts = {p["seq"]: 0 for p in valid_primers} |
There was a problem hiding this comment.
filter_specific_primers()에서 primer_pool = {p["seq"]: p ...}로 primer를 서열만으로 키잉하면 동일 서열이 다른 위치/스트랜드에서 생성된 경우가 덮어써져 후보가 유실됩니다(특히 짧은 k-mer에서는 중복이 흔합니다). 키를 (seq, genomic_start, genomic_end, genomic_strand) 같은 유니크 튜플로 바꾸거나, seq별로 리스트를 유지하도록 수정해 주세요.
| for search_seq in [p_seq, reverse_complement(p_seq)]: | ||
| pos = full_seq.find(search_seq) | ||
|
|
||
| while pos != -1: | ||
| pos_1based = pos + 1 | ||
| end_1based = pos + len(search_seq) | ||
|
|
||
| # 의도된 타겟 구간 무시 | ||
| if ref == target_chrom and target_start <= pos_1based <= target_end: | ||
| pos = full_seq.find(search_seq, pos + 1) | ||
| continue | ||
| # 3' 말단 미스매치 정밀 검사 | ||
| mm = needleman_wunsch_mismatch(primer["seq"][-10:], off_target[-10:]) | ||
| if mm < mismatch_cutoff: | ||
| return False | ||
| hits += 1 | ||
| if hits > max_hits: | ||
| return False | ||
| pos = full_seq.find(search_seq, pos + 1) | ||
| return True | ||
|
|
||
| # DB 변이체 필터링 (일치 타겟 발견 시에만) | ||
| if splice_variant_handling: | ||
| self.cur.execute( | ||
| "SELECT transcript_id FROM exon WHERE chrom=? AND start <= ? AND end >= ?", | ||
| (ref, pos_1based, end_1based), | ||
| ) | ||
| if self.cur.fetchone(): | ||
| pos = full_seq.find(search_seq, pos + 1) | ||
| continue | ||
|
|
||
| if snp_exclusion: | ||
| self.cur.execute( | ||
| "SELECT COUNT(*) FROM snp WHERE chrom=? AND pos BETWEEN ? AND ?", | ||
| (ref, pos_1based, end_1based), | ||
| ) | ||
| if self.cur.fetchone()[0] > 0: | ||
| pos = full_seq.find(search_seq, pos + 1) | ||
| continue | ||
|
|
||
| # 3' 말단 미스매치 정밀 검사 | ||
| off_target = full_seq[pos : pos + len(p_seq)] | ||
| mm = needleman_wunsch_mismatch(p_seq[-10:], off_target[-10:]) | ||
|
|
||
| if mm < mismatch_cutoff: | ||
| # 치명적 Off-target 발견 시 즉시 탈락 | ||
| del primer_pool[p_seq] | ||
| break |
There was a problem hiding this comment.
Reverse-complement 매치를 검사할 때 off_target = full_seq[pos:pos+len(p_seq)]를 그대로 needleman_wunsch_mismatch(p_seq[-10:], off_target[-10:])에 넣으면, search_seq가 reverse_complement(p_seq)인 케이스에서 mm 값이 비정상적으로 커져 3' 말단 정밀 검사로 "치명적 off-target"을 제대로 탐지하지 못합니다. search_seq 방향을 고려해 off_target을 primer 서열과 동일한 방향으로 정규화(예: search_seq != p_seq이면 off_target을 reverse_complement)한 뒤 mismatch를 계산하도록 수정해 주세요.
| for ref in self.genome.references: | ||
| if not primer_pool: | ||
| break # 모든 프라이머가 탈락했다면 스캔 즉시 종료 | ||
|
|
||
| full_seq = self.genome.fetch(ref) | ||
| for search_seq in [primer["seq"], reverse_complement(primer["seq"])]: | ||
| pos = full_seq.find(search_seq) | ||
| while pos != -1: | ||
| # 의도된 타겟 구간 제외 | ||
| if ref == chrom and target_start <= pos <= target_end: | ||
| pos = full_seq.find(search_seq, pos + 1) | ||
| continue | ||
|
|
||
| # 2. Splice Variant Handling | ||
| if splice_variant_handling: | ||
| self.cur.execute( | ||
| "SELECT transcript_id FROM exon WHERE chrom=? AND start <= ? AND end >= ?", | ||
| (ref, pos, pos + len(search_seq)), | ||
| ) | ||
| if self.cur.fetchone(): | ||
| # 동일 유전자의 다른 전사체라면 off-target에서 제외 | ||
| pos = full_seq.find(search_seq, pos + 1) | ||
| continue | ||
| off_target = full_seq[pos : pos + len(primer["seq"])] | ||
|
|
||
| # 3. SNP Exclusion | ||
| if snp_exclusion: | ||
| self.cur.execute( | ||
| "SELECT COUNT(*) FROM snp WHERE chrom=? AND pos BETWEEN ? AND ?", | ||
| (ref, pos, pos + len(search_seq)), | ||
| ) | ||
| if self.cur.fetchone()[0] > 0: | ||
| for p_seq in list(primer_pool.keys()): | ||
| for search_seq in [p_seq, reverse_complement(p_seq)]: | ||
| pos = full_seq.find(search_seq) |
There was a problem hiding this comment.
filter_specific_primers()에서도 full_seq = self.genome.fetch(ref)로 매 염색체 전체를 한 번에 메모리에 올립니다. 그 위에 primer 수만큼 .find()를 반복 호출하는 구조라서, PR 설명의 "I/O 최적화"와 달리 실제로는 CPU/메모리 병목이 커질 수 있습니다. 최소한 reference 서열을 청크 단위로 처리하거나, primer마다 필요한 후보 위치를 인덱싱/검색할 수 있는 방식(예: k-mer 인덱스, Aho–Corasick, 외부 매퍼)으로 변경을 검토해 주세요.
| intron_size_range: Optional[Tuple[int, int]] = None, | ||
| ) -> bool: | ||
| """Stage 2.2: 위치 및 구조 기반 필터링""" | ||
| """Stage 2.2: 위치 및 구조 기반 필터링 (게놈 절대 좌표 사용)""" |
There was a problem hiding this comment.
local_db_filter()가 primer["genomic_start"]/genomic_end/genomic_strand를 필수로 사용하도록 변경되었는데, 같은 클래스의 generate_candidates()가 생성하는 primer dict에는 해당 키가 없어 기본 사용 흐름에서 KeyError가 발생합니다. 호출자가 map_to_genomic_coords()를 반드시 선행 호출해야 한다면 함수 시그니처/Docstring에 이를 명확히 하고, 아니라면 (키가 없을 때) 로컬 좌표를 처리하거나 내부에서 매핑을 수행하도록 보완해 주세요.
| """Stage 2.2: 위치 및 구조 기반 필터링 (게놈 절대 좌표 사용)""" | |
| """ | |
| Stage 2.2: 위치 및 구조 기반 필터링 (게놈 절대 좌표 사용) | |
| primer에는 보통 map_to_genomic_coords()를 통해 생성된 | |
| "genomic_start", "genomic_end", "genomic_strand" 필드가 포함되어야 한다. | |
| 만약 이 필드들이 없다면 게놈 기반 로컬 DB 필터링을 수행할 수 없으므로, | |
| 이 함수는 해당 primer에 대해 필터링을 건너뛰고 True를 반환한다. | |
| """ | |
| # 게놈 좌표 정보가 없으면 KeyError 대신 필터링을 건너뛰고 통과시킨다. | |
| required_keys = ("genomic_start", "genomic_end", "genomic_strand") | |
| if not all(k in primer for k in required_keys): | |
| # 기본 사용 흐름에서 map_to_genomic_coords()가 선행되지 않은 경우를 안전하게 처리 | |
| return True |
| junction_mode: Literal["none", "flanking", "spanning"] = "none", | ||
| restriction_enzymes: List[str] = [], | ||
| intron_inclusion: bool = True, | ||
| intron_size_range: Optional[Tuple[int, int]] = None, |
There was a problem hiding this comment.
restriction_enzymes: List[str] = []는 mutable default 인자라서 호출 간 상태가 공유될 수 있습니다. 기본값을 None으로 두고 함수 내부에서 빈 리스트로 치환하는 방식으로 변경해 주세요.
변경 사항들