203 lines
7.1 KiB
Python
203 lines
7.1 KiB
Python
|
import numpy as np
|
|||
|
|
import pandas as pd
|
|||
|
|
from typing import Tuple, Optional
|
|||
|
|
from Bio import SeqIO
|
|||
|
|
from Bio.Seq import Seq
|
|||
|
|
|
|||
|
|
def fscanr(blastx_output: pd.DataFrame,
|
|||
|
|
mismatch_cutoff: float = 10,
|
|||
|
|
evalue_cutoff: float = 1e-5,
|
|||
|
|
frameDist_cutoff: float = 10) -> pd.DataFrame:
|
|||
|
|
"""
|
|||
|
|
identify PRF sites from BLASTX output
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
blastx_output: BLASTX output DataFrame
|
|||
|
|
mismatch_cutoff: mismatch threshold
|
|||
|
|
evalue_cutoff: E-value threshold
|
|||
|
|
frameDist_cutoff: frame distance threshold
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
pd.DataFrame: DataFrame containing PRF site information
|
|||
|
|
"""
|
|||
|
|
blastx = blastx_output.copy()
|
|||
|
|
|
|||
|
|
blastx.columns = ["qseqid", "sseqid", "pident", "length", "mismatch",
|
|||
|
|
"gapopen", "qstart", "qend", "sstart", "send",
|
|||
|
|
"evalue", "bitscore", "qframe", "sframe"]
|
|||
|
|
|
|||
|
|
blastx = blastx[
|
|||
|
|
(blastx['evalue'] <= evalue_cutoff) &
|
|||
|
|
(blastx['mismatch'] <= mismatch_cutoff)
|
|||
|
|
].dropna()
|
|||
|
|
|
|||
|
|
freq = blastx['qseqid'].value_counts()
|
|||
|
|
multi_hits = freq[freq > 1].index
|
|||
|
|
blastx = blastx[blastx['qseqid'].isin(multi_hits)]
|
|||
|
|
|
|||
|
|
blastx = blastx.sort_values(['qseqid', 'sseqid', 'qstart'])
|
|||
|
|
|
|||
|
|
prf_list = []
|
|||
|
|
for i in range(1, len(blastx)):
|
|||
|
|
curr = blastx.iloc[i]
|
|||
|
|
prev = blastx.iloc[i-1]
|
|||
|
|
|
|||
|
|
if (curr['qseqid'] == prev['qseqid'] and
|
|||
|
|
curr['sseqid'] == prev['sseqid'] and
|
|||
|
|
curr['qframe'] != prev['qframe'] and
|
|||
|
|
curr['qframe'] * prev['qframe'] > 0):
|
|||
|
|
|
|||
|
|
if curr['qframe'] > 0 and prev['qframe'] > 0:
|
|||
|
|
frame_start = prev['qend']
|
|||
|
|
frame_end = curr['qstart']
|
|||
|
|
pep_start = prev['send']
|
|||
|
|
pep_end = curr['sstart']
|
|||
|
|
strand = "+"
|
|||
|
|
elif curr['qframe'] < 0 and prev['qframe'] < 0:
|
|||
|
|
frame_start = prev['qstart']
|
|||
|
|
frame_end = curr['qend']
|
|||
|
|
pep_start = curr['send']
|
|||
|
|
pep_end = prev['sstart']
|
|||
|
|
strand = "-"
|
|||
|
|
else:
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
q_dist = frame_end - frame_start - 1
|
|||
|
|
s_dist = pep_end - pep_start
|
|||
|
|
fs_type = q_dist + (1 - s_dist) * 3
|
|||
|
|
|
|||
|
|
if (abs(q_dist) <= frameDist_cutoff and
|
|||
|
|
abs(s_dist) <= frameDist_cutoff // 3 and
|
|||
|
|
-3 < fs_type < 3):
|
|||
|
|
|
|||
|
|
prf_list.append({
|
|||
|
|
'DNA_seqid': curr['qseqid'],
|
|||
|
|
'FS_start': frame_start,
|
|||
|
|
'FS_end': frame_end,
|
|||
|
|
'Pep_seqid': curr['sseqid'],
|
|||
|
|
'Pep_FS_start': prev['send'] + 1,
|
|||
|
|
'Pep_FS_end': curr['sstart'],
|
|||
|
|
'FS_type': fs_type,
|
|||
|
|
'Strand': strand
|
|||
|
|
})
|
|||
|
|
|
|||
|
|
if not prf_list:
|
|||
|
|
print("No PRF events detected!")
|
|||
|
|
return pd.DataFrame()
|
|||
|
|
|
|||
|
|
prf = pd.DataFrame(prf_list)
|
|||
|
|
|
|||
|
|
for col in ['DNA_seqid', 'Pep_seqid']:
|
|||
|
|
for pos in ['FS_start', 'FS_end']:
|
|||
|
|
loci = prf[col] + '_' + prf[pos].astype(str)
|
|||
|
|
prf = prf[~loci.duplicated()]
|
|||
|
|
|
|||
|
|
return prf
|
|||
|
|
|
|||
|
|
def extract_prf_regions(mrna_file: str, prf_data: pd.DataFrame) -> pd.DataFrame:
|
|||
|
|
"""
|
|||
|
|
从mRNA序列中提取PRF位点周围的序列
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
mrna_file: mRNA序列文件路径 (FASTA格式)
|
|||
|
|
prf_data: FScanR输出的PRF位点数据
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
pd.DataFrame: 包含399bp序列的DataFrame
|
|||
|
|
"""
|
|||
|
|
mrna_dict = {rec.id: str(rec.seq)
|
|||
|
|
for rec in SeqIO.parse(mrna_file, "fasta")}
|
|||
|
|
|
|||
|
|
results = []
|
|||
|
|
for _, row in prf_data.iterrows():
|
|||
|
|
seq_id = row['DNA_seqid']
|
|||
|
|
if seq_id not in mrna_dict:
|
|||
|
|
print(f"警告: {seq_id} 未在mRNA文件中找到")
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
sequence = mrna_dict[seq_id]
|
|||
|
|
strand = row['Strand']
|
|||
|
|
fs_start = int(row['FS_start'])
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
if strand == '-':
|
|||
|
|
sequence = str(Seq(sequence).reverse_complement())
|
|||
|
|
|
|||
|
|
# 只提取399bp序列，33bp由predictor内部截取
|
|||
|
|
full_seq = extract_window_sequences(sequence, fs_start)[1]
|
|||
|
|
|
|||
|
|
results.append({
|
|||
|
|
'DNA_seqid': seq_id,
|
|||
|
|
'FS_start': fs_start,
|
|||
|
|
'FS_end': int(row['FS_end']),
|
|||
|
|
'Strand': strand,
|
|||
|
|
'399bp': full_seq,
|
|||
|
|
'FS_type': row['FS_type']
|
|||
|
|
})
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"处理 {seq_id} 时出错: {str(e)}")
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
return pd.DataFrame(results)
|
|||
|
|
|
|||
|
|
def extract_window_sequences(seq: str, position: int) -> Tuple[Optional[str], Optional[str]]:
|
|||
|
|
"""
|
|||
|
|
从指定位置提取分析窗口序列
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
seq: 输入DNA序列
|
|||
|
|
position: 当前分析位置 (FS_start)
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
Tuple[str, str]: (33bp序列, 399bp序列) - 已调整为与训练模型匹配的长度
|
|||
|
|
"""
|
|||
|
|
# 确保位置在密码子边界上（整数倍的3）
|
|||
|
|
frame_position = position - (position % 3)
|
|||
|
|
|
|||
|
|
# 计算33bp窗口的起止位置 (GB模型)
|
|||
|
|
half_size_small = 33 // 2
|
|||
|
|
start_small = frame_position - half_size_small
|
|||
|
|
end_small = frame_position + half_size_small + (33 % 2) # 添加余数以处理奇数长度
|
|||
|
|
|
|||
|
|
# 计算399bp窗口的起止位置 (CNN模型)
|
|||
|
|
half_size_large = 399 // 2
|
|||
|
|
start_large = frame_position - half_size_large
|
|||
|
|
end_large = frame_position + half_size_large + (399 % 2) # 添加余数以处理奇数长度
|
|||
|
|
|
|||
|
|
# 提取序列并填充
|
|||
|
|
seq_small = _extract_and_pad(seq, start_small, end_small, 33)
|
|||
|
|
seq_large = _extract_and_pad(seq, start_large, end_large, 399)
|
|||
|
|
|
|||
|
|
return seq_small, seq_large
|
|||
|
|
|
|||
|
|
def _extract_and_pad(seq: str, start: int, end: int, target_length: int) -> str:
|
|||
|
|
"""提取序列并用N填充"""
|
|||
|
|
if start < 0:
|
|||
|
|
prefix = 'N' * abs(start)
|
|||
|
|
extracted = prefix + seq[:end]
|
|||
|
|
elif end > len(seq):
|
|||
|
|
suffix = 'N' * (end - len(seq))
|
|||
|
|
extracted = seq[start:] + suffix
|
|||
|
|
else:
|
|||
|
|
extracted = seq[start:end]
|
|||
|
|
|
|||
|
|
# 确保序列长度正确
|
|||
|
|
if len(extracted) < target_length:
|
|||
|
|
# 从中心填充
|
|||
|
|
pad_left = (target_length - len(extracted)) // 2
|
|||
|
|
pad_right = target_length - len(extracted) - pad_left
|
|||
|
|
extracted = 'N' * pad_left + extracted + 'N' * pad_right
|
|||
|
|
elif len(extracted) > target_length:
|
|||
|
|
# 从序列两端等量截取
|
|||
|
|
excess = len(extracted) - target_length
|
|||
|
|
trim_each_side = excess // 2
|
|||
|
|
extracted = extracted[trim_each_side:len(extracted)-trim_each_side]
|
|||
|
|
|
|||
|
|
return extracted
|
|||
|
|
|
|||
|
|
def prepare_cnn_input(sequence: str) -> np.ndarray:
|
|||
|
|
"""prepare CNN model input"""
|
|||
|
|
base_to_num = {'A': 1, 'T': 2, 'G': 3, 'C': 4, 'N': 0}
|
|||
|
|
seq_numeric = [base_to_num.get(base, 0) for base in sequence.upper()]
|
|||
|
|
return np.array(seq_numeric).reshape(1, len(sequence), 1)
|