NCR+可变电荷块3——NCB/cell绘图1
文献method参考：
蛋白质序列数据从uniprot中获取
https://www.uniprot.org/uniprotkb/P46013/entry
https://www.uniprot.org/uniprotkb/P06748/entry、
1，电荷分布计算：
Charge distribution was calculated as the sum of the charges (Arg and Lys, +1; Glu and Asp, −1; phospho-Ser and phospho-Thr, −2) in the indicated window range. A charge block was designated when the area of the charge plot (window size: 35 amino acids) was larger than 20.
（1）电荷分布是计算指定窗口window内净电荷之和，所以dict对应关系是1个window对应1个电荷总和，所以电荷分布图实际上就是绘制（x，y）=（窗口编号，该窗口净电荷总和）的折线图
①首先计算给定窗口内每个氨基酸的净电荷是简单的：
碱性带正电aa：KR，但是没有组氨酸H？
酸性带负电aa：DE
至于磷酸化的氨基酸没有数据可以验证（我们只有序列数据，没法获取，需要实际in vivo数据，也许可以从蛋白质晶体数据库上获取？）
可以另外开一个磷酸化的字典dict，存放不同的aa字符与带电荷的对应关系，其实这样的PTM的字典可以开很多
**②问题是计算电荷分布的时候窗口开多大？
**是method里指定的win=35吗？
这里很容易被误导，实际上不是35。
35指的是charge plot已经绘制出来之后，然后计算特定区域内charge plot的area面积，实际上应该需要进行积分。
有个问题，如何理解这里的area>20，首先要理解block，肯定是单一电荷性质的，所以肯定是看的积分之后的净area，比如说是计算对应的正电荷的正面积之后，以及负电荷的负面积，然后计算这个窗口的电荷总和，所以实际上就是电荷面积计算的一个净net值的绝对值。
**所以在计算charge plot中的window窗口应该设置为多大？
**
可以看到实际上每一幅图指定的window窗口大小尺寸都是25
但是在提供的附近材料（审稿意见答复）中又说是35
当然仅仅只是获取charge plot的话窗口定义多大多小都不影响，因为只是展示：
看复现效果：Ki-67
16个RD重复单元，每个大概在110个aa左右
我们只看R12
实际上可以看到对应紫色那块，也就是K167R 12
实际上大小是112
示例代码如下：
def calculate_charge_distribution(sequence, phospho_sites, window_size=35):"""计算蛋白质电荷分布Args:sequence: 蛋白质氨基酸序列 (字符串)phospho_sites: 磷酸化位点列表 (例如：[10, 25, 50])，位点编号从1开始，蛋白质序列index是1-indexwindow_size: 滑动窗口大小 (整数)Returns:一个列表，包含每个窗口的净电荷。"""charges = {'R': 1, 'K': 1, 'E': -1, 'D': -1, 'S': 0, 'T': 0} #因为原文考虑到了丝氨酸以及苏氨酸的磷酸化与否对蛋白质电荷块分布的影响，所以这里需要对丝氨酸S以及苏氨酸T设置2个dict#但是正电荷没有考虑到组氨酸H，所以这里没有设置组氨酸的正电荷phospho_charges = { 'S': -2, 'T':-2}#这个PTM翻译后修饰的字典可以根据实际情况进行修改，进行扩展sequence = sequence.upper()# 将序列转换为大写，以确保一致性charge_distribution = []# 初始化一个空列表来存储每个窗口的净电荷for i in range(len(sequence) - window_size + 1):# 遍历的窗口编号，从0开始window = sequence[i:i + window_size]# 获取当前窗口的序列total_charge = 0# 初始化当前窗口的总电荷for j, aa in enumerate(window):# enumerate(iterable, start=0)，用于在遍历可迭代对象（如列表、字符串等）时，同时获取元素的索引和值，索引默认从0开始# j是索引从0开始，aa是氨基酸字符pos = i+j+1 #对应氨基酸的真实物理位置1-based，注意i是range从0开始（0:第1个窗口，j是enumrate返回的元组的索引，从0开始），所以pos=i+j+1#编号i=0的窗口、索引j=0的序列，实际上是第1个氨基酸pos=1charge = charges.get(aa, 0)# 需要计算电荷的氨基酸都在charges字典里，如果不在字典里，返回0if pos in phospho_sites:charge = phospho_charges.get(aa,0)#如果当前氨基酸在磷酸化位点列表里，则用磷酸化后的电荷替换掉原来的电荷total_charge += charge# 将氨基酸的电荷加到当前窗口的总电荷中charge_distribution.append(total_charge)# 将当前窗口的总电荷添加到charge_distribution列表中return charge_distribution# 返回电荷分布列表# 示例用法:
sequence = "KTTKIACKSPQPDPVDTPASTKQRPKRNLRKADVEEEFLALRKRTPSAGKAMDTPKPAVSDEKNINTFVETPVQKLDLLGNLPGSKRQPQTPKEKAEALEDLVGFKELFQTP" #替换为你的氨基酸序列
phospho_sites = [] #替换为你的磷酸化位点
charge_distribution = calculate_charge_distribution(sequence, phospho_sites)
print(charge_distribution)#绘制图形，需要用到matplotlib库
import matplotlib.pyplot as plt
plt.plot(charge_distribution)
plt.xlabel("Residue")
plt.ylabel("Charge")
plt.title("Charge plot")
plt.show()如果窗口是35内计算的净电荷，那R12：
如果是25窗口大小
很显然不对，因为RD12序列就长112，如果是使用窗口的话，实际上是无法延伸到100多位置的氨基酸序列的
所以实际上都和实物图对不上
其实看下面这张图就很清晰了：
这里是105到250都有，所以肯定不是这么画的
——》如果只是单纯看实物图的话，其实实物图肯定是每一个氨基酸位置都考虑进去的，
并且每个位置上的氨基酸序列肯定是计算了不止一个电荷的（所以x轴肯定又不是表征的单个氨基酸序列位置的信息，应该是某个集合或者是窗口下的位置信息之和）
问题就在于如果考虑了所有窗口内的氨基酸序列之后，
那几张又如何计算整个窗口下的
——》索性查看引用这篇NCB的所有文献里的charge plot怎么画的：
https://www.bioinformatics.nl/cgi-bin/emboss/charge
如果我只是输入同样的R12序列，结果是：
确实是能够计算到每一个位置，且电荷显然不是只计算一个，既然肯定有方法计算
我们重点关注的是下面NCPR的图，
we calculated the net charge per residue (NCPR) using a sliding window of 10 residues
没有关于电荷的图，但是倒是有一个有趣的图
总之我们能够看到，计算电荷块，看的多是整体区域的带电性质，也就是NPCR，一般不会用来看DE或者KRH数目的；
所以应该先搞清楚NCPR概念，然后再绘制charge plot（有一个参考来源即可）
参考2010的PNAS：
f+定义为一段给定aa序列中正电荷残基比例（理论上比例应该是这样算的：
f+=（带正电荷残基数x该残基所带正电荷电荷数）/该aa序列长度，
同理f-，
因为KRH和DE都是单位带电荷数，所以实际上就是看aa数目比例
以下图第一段序列为例：
全长29，正电荷氨基酸数9，负电荷为3，所以
f+=9/29=0.31
f-=3/29=0.10
其实就是将正负电荷均分到窗口/全长序列的每个残基中，比如说第一段序列相当于每个残基有31%的正电，同时带10%的负电，净电荷就是21%的正电荷
对应很多文献滑动窗口下的charge曲线实际上有正有负，定义NCPR=f±f-
（2）如何由电荷分布charge plot推导出电荷块charge block：
实际上就是如何表征电荷块
A charge block was designated when the area of the charge plot (window size: 35 amino acids) was larger than 20.
然后单个净电荷如何应用于计算电荷块：
设置一个窗口大小为35aa，但是是什么大于20呢：
带净电荷的aa数目？不合理，万一正的10个，负的10个 55开呢？
sum之后的净电荷值（应该是绝对值）？合理——》进一步评估
（3）必要时可可以进行多序列比对：
用的还是Clustal Omega（ https://www.ebi.ac.uk/Tools/msa/clustalo/ ）
（4）自定义的电荷模式参数：charge patterning parameters
文献里最终使用的是Blc参数
Dseg：
Blc：
另外还可以看一下reporting summary，实际上看到使用的数据分析方面的工具都是一些普通的python库
一些更细节的地方可以深入分析：一些问题都被peer提出来了
2，电荷模式参数计算：
实际使用参数评估过程中有很多参数可以使用：
kappa参数
SCD参数
Dseg参数
Blc参数
==================================================================
鉴于NCB文献中对于block的模糊定义，还是借鉴cell文献的method：
NCPR依据上面的定义我们能算了，然后窗口就是开10为大小，
window=10，如果NCPR>=5，则为酸性/负电block；
NCPR<=5，则为碱性/正电block
魔改一下我之前写的calc_NCR.py程序
import sys #提供参数
import csv #处理分隔符文件，不一定是csv
import argparse #用于参数解析，帮助手册
import subprocess #用于执行shell命令
# import os #用于文件操作# 定义一个函数来解析FASTA文件,主要是获得1个蛋白质的ID和氨基酸序列对应的字典sequences
def parse_fasta(fasta_file): with open(fasta_file, 'r') as file:#使用with语句打开文件，'r'表示以只读模式打开sequences = {} current_id = Nonecurrent_sequence = []#初始化了一个空字典sequences来存储蛋白质序列:一个变量current_id来存储当前处理的蛋白质的ID，以及一个列表current_sequence来存储当前蛋白质的氨基酸序列for line in file:line = line.strip()#遍历文件中的每一行，并使用strip()方法去除每行的换行符，这样后面append的时候所有序列会连在一起if line.startswith('>'):#检查当前行是否以'>'字符开头，如果是，那么这一行就是一个新的蛋白质记录的开始if current_id:sequences[current_id] = ''.join(current_sequence)current_sequence = []#如果current_id不为空（即之前已经处理过一个蛋白质），那么将当前蛋白质的序列（current_sequence）添加到sequences字典中，并重置current_sequence列表parts = line[1:].split()#对于每一行，从">"之后开始切割，然后获取['sp|Q6ZN18|AEBP2_HUMAN', 'Zinc', 'finger', 'protein', 'AEBP2', 'OS=Homo', 'sapiens', 'OX=9606', 'GN=AEBP2', 'PE=1', 'SV=2']id_info = parts[0].split('|')current_id = id_info[2] #实际上获取的sp|Q6ZN18|AEBP2_HUMAN中第3个分隔符字段才是对应的蛋白质名字，当然也可以对照GN字段后面的基因名字else:current_sequence.append(line)#对应前面的if，如果当前行不是以'>'字符开头，那么这一行就是一个蛋白质的氨基酸序列，将其添加到current_sequence列表中，直接添加并合并if current_id:sequences[current_id] = ''.join(current_sequence)#确保当前序列已经有一个ID（即 current_id 不为空，然后直接存储对应的序列到字典中return sequencesdef sliding_window_for_charged_block(sequence, phospho_sites=[],window_size=10, threshold=5):"""滑动窗口算法，用于查找蛋白质序列中NCPR绝对值大于阈值的窗口Args:sequence (str): 蛋白质序列phospho_sites (list): 磷酸化位点列表，默认为空,因为main主函数中是批量化处理,但显然每一个蛋白质序列都有差异的磷酸化位点，如果要分开计算，这个函数模块可以保留下来，但需要修改main函数window_size (int): 窗口大小,默认为10threshold (float): NCPR阈值,默认为5Returns:list: 合并后的窗口起点和终点列表,即找出的block"""charges = {'R': 1, 'K': 1,'H':1, 'E': -1, 'D': -1, 'S': 0, 'T': 0} #因为原文考虑到了丝氨酸以及苏氨酸的磷酸化与否对蛋白质电荷块分布的影响，所以这里需要对丝氨酸S以及苏氨酸T设置2个dict#NCB文献里正电荷没有考虑到组氨酸H，但是cell里有，还是补上phospho_charges = { 'S': -2, 'T':-2}#这个PTM翻译后修饰的字典可以根据实际情况进行修改，进行扩展results_acidic = [] #用于保存找到的酸性block的起点和终点，结果起点终点对保存在列表中results_basic = [] #用于保存找到的碱性block的起点和终点for i in range(len(sequence) - window_size + 1):# 遍历的窗口编号，从0开始 window = sequence[i:i + window_size]# 获取当前窗口的序列，注意取的是str切片，终点取不到，但是终点-起点+1-1total_charge = 0# 初始化当前窗口的总电荷for j, aa in enumerate(window):# enumerate(iterable, start=0)，用于在遍历可迭代对象（如列表、字符串等）时，同时获取元素的索引和值，索引默认从0开始# j是索引从0开始，aa是氨基酸字符pos = i+j+1 #对应氨基酸的真实物理位置1-based，注意i是range从0开始（0:第1个窗口，j是enumrate返回的元组的索引，从0开始），所以pos=i+j+1#编号i=0的窗口、索引j=0的序列，实际上是第1个氨基酸pos=1charge = charges.get(aa, 0)# 需要计算电荷的氨基酸都在charges字典里，如果不在字典里，返回0if pos in phospho_sites:charge = phospho_charges.get(aa,0)#如果当前氨基酸在磷酸化位点列表里，则用磷酸化后的电荷替换掉原来的电荷total_charge += charge# 计算窗口内的总电荷if abs(total_charge) > threshold:# 窗口内的绝对值大于阈值，则判定为charged blockif total_charge > 0:results_basic.append((i+1, i + window_size))# 该起点终点对保存在正电碱性block列表中，(i+1, i + window_size)保存的数据类型是元组tuple，且是实际1-based的physical位置else:results_acidic.append((i+1, i + window_size))# 该起点终点对保存在酸性block列表中，(i+1, i + window_size)保存的数据类型是元组tuple# 合并重叠的窗口merged_results_acidic = []merged_results_basic =[] #同前面，对酸性负电以及碱性正电氨基酸的窗口进行合并，设置初始化的空列表# 用于保存已经合并的窗口起点终点元组对，注意每次比较都是比较当前元组的起点和上一个已经合并的元组（也就是merged_results[-1])的终点#下面是合并酸性带负电的窗口if results_acidic:# 如果酸性负电窗口results列表不为空merged_results_acidic = [results_acidic[0]]# 将第一个窗口加入merged_results列表中，作为后续比对的起点for current in results_acidic[1:]:# 遍历results_acidic列表中剩余的窗口last = merged_results_acidic[-1]# 取出已经合并的窗口的最后一个窗口的元组对if current[0] <= last[1] + 1: # 如果当前窗口与最后一个窗口有重叠，也就是当前窗口的起点小于等于已经合并的最后一个窗口的终点+1（overlap以及相邻）merged_results_acidic[-1] = (last[0], max(last[1], current[1]))# 更新已经合并的最后一个窗口的终点为当前窗口和已经合并的最后一个窗口的终点的最大值，总之就是更新最后一个窗口，起点不变，终点更新else:merged_results_acidic.append(current)# 如果当前窗口与最后一个窗口没有重叠，直接将当前窗口加入merged_results列表中，也就是当前窗口并没有重叠#下面是合并碱性带正电的窗口if results_basic:merged_results_basic = [results_basic[0]]for current in results_basic[1:]:last = merged_results_basic[-1]if current[0] <= last[1] + 1:merged_results_basic[-1] = (last[0], max(last[1], current[1]))else:merged_results_basic.append(current)return merged_results_acidic, merged_results_basic# 返回合并之后的酸性负电窗口和碱性正电窗口，各自是合并后的窗口起点和终点列表，且每一个元素是元组，后续需要解包赋值def main(fasta_file, output_file,window_size, threshold, delimiter='\t'): #默认输出为tsv格式sequences = parse_fasta(fasta_file)#调用parse_fasta函数来解析FASTA文件，并返回一个字典sequences，其中键是蛋白质的ID，值是蛋白质的氨基酸序列（字符串）#处理正电荷block、负电荷block区域输出文件with open(output_file, 'w') as csvfile:#使用with语句打开一个输出文件（需要自己提供参数命名，'w'表示以写入模式打开writer = csv.writer(csvfile, delimiter=delimiter)#创建一个csv.writer对象，用于将数据写入csv文件，delimiter='\t'表示使用制表符作为字段分隔符writer.writerow(['Protein_Name', 'charged_block_start', 'charged_block_end', 'Sequence', 'charged_block_type'])#写入表头，即列名for protein_name, sequence in sequences.items():#sequences是一个字典，键是蛋白质的ID，值是蛋白质的氨基酸序列，现在从值中获取charged block区域起点以及终点的元组对列表merged_results_acidic, merged_results_basic= sliding_window_for_charged_block(sequence.replace('\n', ''),phospho_sites=[],window_size=window_size, threshold=threshold) #.replace('\n', '')双重保证，如果前面fasta序列中没有去除换行符for results in merged_results_acidic:writer.writerow([protein_name, results[0], results[1], sequence[results[0]-1:results[1]], 'Acidic'])for results in merged_results_basic:writer.writerow([protein_name, results[0], results[1], sequence[results[0]-1:results[1]], 'Basic'])#将蛋白质的ID、charged block区域的起点和终点，以及对应的序列切片,以及是酸性还是碱性，写入csv文件中，注意切片是0-based的，而真实位置是1-based的#后续可以导入到R中使用tidyverse#如果需要统计信息并写入到{output_file}.stats文件中，可以参考我原本calc_NCR.py的代码，这里不再赘述 #使用 argparse 模块来实现类似于 Linux 工具的命令行参数解析（包括生成帮助手册等）if __name__ == "__main__":parser = argparse.ArgumentParser(description="FASTA文件处理工具")parser.add_argument("--fasta_file", help="输入的FASTA文件路径")parser.add_argument("--delimiter", help="输出文件的分隔符（如$'\t'等）")parser.add_argument("--window_size", type=int, default=10, help="滑动窗口大小")parser.add_argument("--threshold", type=float, default=5, help="NCPR阈值,注意是>=")parser.add_argument("--output_file", help="输出可变电荷块文件路径")# parser.add_argument("--help", action="help", help="显示帮助信息")，不需要手动添加 --help 选项，因为 argparse 模块已经自动添加了一个 --help 选项来显示帮助信息qargs = parser.parse_args()main(args.fasta_file, args.output_file, args.delimiter, args.window_size, args.threshold)
去掉符号修改之后是
import sys #提供参数
import csv #处理分隔符文件，不一定是csv
import argparse #用于参数解析，帮助手册
import subprocess #用于执行shell命令
# import os #用于文件操作# 定义一个函数来解析FASTA文件,主要是获得1个蛋白质的ID和氨基酸序列对应的字典sequences
def parse_fasta(fasta_file): with open(fasta_file, 'r') as file:#使用with语句打开文件，'r'表示以只读模式打开sequences = {} current_id = Nonecurrent_sequence = []#初始化了一个空字典sequences来存储蛋白质序列:一个变量current_id来存储当前处理的蛋白质的ID，以及一个列表current_sequence来存储当前蛋白质的氨基酸序列for line in file:line = line.strip()#遍历文件中的每一行，并使用strip()方法去除每行的换行符，这样后面append的时候所有序列会连在一起if line.startswith('>'):#检查当前行是否以'>'字符开头，如果是，那么这一行就是一个新的蛋白质记录的开始if current_id:sequences[current_id] = ''.join(current_sequence)current_sequence = []#如果current_id不为空（即之前已经处理过一个蛋白质），那么将当前蛋白质的序列（current_sequence）添加到sequences字典中，并重置current_sequence列表parts = line[1:].split()#对于每一行，从">"之后开始切割，然后获取['sp|Q6ZN18|AEBP2_HUMAN', 'Zinc', 'finger', 'protein', 'AEBP2', 'OS=Homo', 'sapiens', 'OX=9606', 'GN=AEBP2', 'PE=1', 'SV=2']id_info = parts[0].split('|')current_id = id_info[2] #实际上获取的sp|Q6ZN18|AEBP2_HUMAN中第3个分隔符字段才是对应的蛋白质名字，当然也可以对照GN字段后面的基因名字else:current_sequence.append(line)#对应前面的if，如果当前行不是以'>'字符开头，那么这一行就是一个蛋白质的氨基酸序列，将其添加到current_sequence列表中，直接添加并合并if current_id:sequences[current_id] = ''.join(current_sequence)#确保当前序列已经有一个ID（即 current_id 不为空，然后直接存储对应的序列到字典中return sequencesdef sliding_window_for_charged_block(sequence, phospho_sites=[],window_size=10, threshold=5):"""滑动窗口算法，用于查找蛋白质序列中NCPR绝对值大于阈值的窗口Args:sequence (str): 蛋白质序列phospho_sites (list): 磷酸化位点列表，默认为空,因为main主函数中是批量化处理,但显然每一个蛋白质序列都有差异的磷酸化位点，如果要分开计算，这个函数模块可以保留下来，但需要修改main函数window_size (int): 窗口大小,默认为10threshold (float): NCPR阈值,默认为5Returns:list: 合并后的窗口起点和终点列表,即找出的block"""charges = {'R': 1, 'K': 1,'H':1, 'E': -1, 'D': -1, 'S': 0, 'T': 0} #因为原文考虑到了丝氨酸以及苏氨酸的磷酸化与否对蛋白质电荷块分布的影响，所以这里需要对丝氨酸S以及苏氨酸T设置2个dict#NCB文献里正电荷没有考虑到组氨酸H，但是cell里有，还是补上phospho_charges = { 'S': -2, 'T':-2}#这个PTM翻译后修饰的字典可以根据实际情况进行修改，进行扩展results_acidic = [] #用于保存找到的酸性block的起点和终点，结果起点终点对保存在列表中results_basic = [] #用于保存找到的碱性block的起点和终点for i in range(len(sequence) - window_size + 1):# 遍历的窗口编号，从0开始 window = sequence[i:i + window_size]# 获取当前窗口的序列，注意取的是str切片，终点取不到，但是终点-起点+1-1total_charge = 0# 初始化当前窗口的总电荷for j, aa in enumerate(window):# enumerate(iterable, start=0)，用于在遍历可迭代对象（如列表、字符串等）时，同时获取元素的索引和值，索引默认从0开始# j是索引从0开始，aa是氨基酸字符pos = i+j+1 #对应氨基酸的真实物理位置1-based，注意i是range从0开始（0:第1个窗口，j是enumrate返回的元组的索引，从0开始），所以pos=i+j+1#编号i=0的窗口、索引j=0的序列，实际上是第1个氨基酸pos=1charge = charges.get(aa, 0)# 需要计算电荷的氨基酸都在charges字典里，如果不在字典里，返回0if pos in phospho_sites:charge = phospho_charges.get(aa,0)#如果当前氨基酸在磷酸化位点列表里，则用磷酸化后的电荷替换掉原来的电荷total_charge += charge# 计算窗口内的总电荷if abs(total_charge) >= threshold:# 窗口内的绝对值大于等于阈值，则判定为charged blockif total_charge > 0:results_basic.append((i+1, i + window_size))# 该起点终点对保存在正电碱性block列表中，(i+1, i + window_size)保存的数据类型是元组tuple，且是实际1-based的physical位置else:results_acidic.append((i+1, i + window_size))# 该起点终点对保存在酸性block列表中，(i+1, i + window_size)保存的数据类型是元组tuple# 合并重叠的窗口merged_results_acidic = []merged_results_basic =[] #同前面，对酸性负电以及碱性正电氨基酸的窗口进行合并，设置初始化的空列表# 用于保存已经合并的窗口起点终点元组对，注意每次比较都是比较当前元组的起点和上一个已经合并的元组（也就是merged_results[-1])的终点#下面是合并酸性带负电的窗口if results_acidic:# 如果酸性负电窗口results列表不为空merged_results_acidic = [results_acidic[0]]# 将第一个窗口加入merged_results列表中，作为后续比对的起点for current in results_acidic[1:]:# 遍历results_acidic列表中剩余的窗口last = merged_results_acidic[-1]# 取出已经合并的窗口的最后一个窗口的元组对if current[0] <= last[1] + 1: # 如果当前窗口与最后一个窗口有重叠，也就是当前窗口的起点小于等于已经合并的最后一个窗口的终点+1（overlap以及相邻）merged_results_acidic[-1] = (last[0], max(last[1], current[1]))# 更新已经合并的最后一个窗口的终点为当前窗口和已经合并的最后一个窗口的终点的最大值，总之就是更新最后一个窗口，起点不变，终点更新else:merged_results_acidic.append(current)# 如果当前窗口与最后一个窗口没有重叠，直接将当前窗口加入merged_results列表中，也就是当前窗口并没有重叠#下面是合并碱性带正电的窗口if results_basic:merged_results_basic = [results_basic[0]]for current in results_basic[1:]:last = merged_results_basic[-1]if current[0] <= last[1] + 1:merged_results_basic[-1] = (last[0], max(last[1], current[1]))else:merged_results_basic.append(current)return merged_results_acidic, merged_results_basic# 返回合并之后的酸性负电窗口和碱性正电窗口，各自是合并后的窗口起点和终点列表，且每一个元素是元组，后续需要解包赋值def main(fasta_file, output_file,window_size, threshold): #默认输出为tsv格式sequences = parse_fasta(fasta_file)#调用parse_fasta函数来解析FASTA文件，并返回一个字典sequences，其中键是蛋白质的ID，值是蛋白质的氨基酸序列（字符串）#处理正电荷block、负电荷block区域输出文件with open(output_file, 'w') as csvfile:#使用with语句打开一个输出文件（需要自己提供参数命名，'w'表示以写入模式打开writer = csv.writer(csvfile, delimiter='\t')#创建一个csv.writer对象，用于将数据写入csv文件，delimiter='\t'表示使用制表符作为字段分隔符writer.writerow(['Protein_Name', 'charged_block_start', 'charged_block_end', 'Sequence', 'charged_block_type'])#写入表头，即列名for protein_name, sequence in sequences.items():#sequences是一个字典，键是蛋白质的ID，值是蛋白质的氨基酸序列，现在从值中获取charged block区域起点以及终点的元组对列表merged_results_acidic, merged_results_basic= sliding_window_for_charged_block(sequence.replace('\n', ''),phospho_sites=[],window_size=window_size, threshold=threshold) #.replace('\n', '')双重保证，如果前面fasta序列中没有去除换行符for results in merged_results_acidic:writer.writerow([protein_name, results[0], results[1], sequence[results[0]-1:results[1]], 'Acidic'])for results in merged_results_basic:writer.writerow([protein_name, results[0], results[1], sequence[results[0]-1:results[1]], 'Basic'])#将蛋白质的ID、charged block区域的起点和终点，以及对应的序列切片,以及是酸性还是碱性，写入csv文件中，注意切片是0-based的，而真实位置是1-based的#后续可以导入到R中使用tidyverse#如果需要统计信息并写入到{output_file}.stats文件中，可以参考我原本calc_NCR.py的代码，这里不再赘述 #使用 argparse 模块来实现类似于 Linux 工具的命令行参数解析（包括生成帮助手册等）if __name__ == "__main__":parser = argparse.ArgumentParser(description="接受uniprot输入的蛋白质原始fasta序列文件(可以是一群蛋白质)，计算每个窗口下的NCPR值，输出NCPR值大于等于阈值的窗口(也就是酸性带正电氨基酸block，或碱性带负电氨基酸block)的起点和终点")parser.add_argument("--fasta_file", help="输入的蛋白质FASTA文件路径")# parser.add_argument("--delimiter", help="输出文件的分隔符（如$'\\t'等）,默认是制表符", default='\t')parser.add_argument("--window_size", type=int, default=10, help="滑动窗口大小，默认是10")parser.add_argument("--threshold", type=float, default=5, help="NCPR阈值，默认是5，注意是>=")parser.add_argument("--output_file", help="输出可变电荷块文件路径")# parser.add_argument("--help", action="help", help="显示帮助信息")，不需要手动添加 --help 选项，因为 argparse 模块已经自动添加了一个 --help 选项来显示帮助信息qargs = parser.parse_args()main(args.fasta_file, args.output_file, args.window_size, args.threshold)
python3 calc_charged_block.py --fasta_file /mnt/sdb/zht/project/uniprot/uniprot_870.fasta --window_size 10 --threshold 5 --output_file uniprot_870_NCPR_10_5.tsv
大概效果如下：
实际上是可以和IDR等各种bed区域信息都统一起来
下一步就是结构预测，或者是提取感兴趣的区域做个结构表位分析（motif，奇美拉/pymol中搜索motif表征之类的）