实验篇——基因家族成员鉴定代码实现可能性(二)_怎么利用blast搜索家族成员

实验篇——基因家族成员鉴定代码实现可能性(二)

---

前言

联系上一章的使用TBtools软件以及NCBI官网完成的基因家族成员鉴定，虽然能够完成鉴定，但是繁琐重复的步骤依然很多，最主要的是要人自己主动去操作，使用一个参考序列完成一次鉴定已经很费时了，还要使用许多参考序列去鉴定（是一个大工程）。
而我想的是能否用代码去实现来代替人的重复操作。我主要有以下三个思路。嗯…,不幸的是，前面两种我进行不下去了，而最后一种虽然能进行下去，但是这无疑是一种笨方法，而且不太稳定(具体原因会在下文解释).

一、思路一

这是我的第一个思路
我当时想的是，进行基因家族鉴定主要的是进行blast对比，那我就要去得到能执行blast对比的代码，下面是我想实现的代码(不完善，放弃了)

import glob
import os
from io import StringIO
import Bio.ExPASy
import Bio.SwissProt
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
import pandas as pd

def split_gene_file(source_file, output_folder, genes_per_file):
    df = pd.read_csv(source_file, sep='\t')
    # 计算需要创建的文件数
    num_files = len(df) // genes_per_file + 1
    os.makedirs(output_folder, exist_ok=True)
    for i in range(num_files):
        # 计算每个文件对应的基因索引范围
        start = i * genes_per_file
        end = start + genes_per_file
        # 创建一个新的DataFrame用于存储当前文件的基因数据
        df_subset = df.iloc[start:end]
        output_file = f"{output_folder}/gene_file_{i + 1}.csv"
        # 将数据保存到新文件中
        df_subset.to_csv(output_file, index=False)
split_gene_file("weiyi_shili.pep", "gene1", 1000)

# 读取参考序列
with open("cankaoshili.fasta", "r", encoding="utf-8-sig") as f:
    cankao_seq = list(SeqIO.parse(f, "fasta"))

def jiyinjiazu_folder(folder_path):
    # 获取文件夹下的所有文件路径
    file_paths = glob.glob(os.path.join(folder_path, "*"))

    for file_path in file_paths:
        # 获取待对比的基因组数据
        known_seqfile = file_path
        known_sequences = []
        with open(known_seqfile, "r", encoding="utf-8") as f:
            known_seqfile = list(SeqIO.parse(f, "fasta"))
        for seq_record in known_seqfile:
            known_sequence = str(seq_record.seq).replace("\n", "")
            known_sequences.append(known_sequence)
       
        # 执行BLAST比对
        result_handle = NCBIWWW.qblast("blastp", "nr", known_sequences, entrez_query=cankao_seq, expect=10, word_size=4,
                                       matrix_name="BLOSUM80")
        result_str = result_handle.read()
        result_handle.close()

        # 将BLAST结果解析为可迭代对象
        result_io = StringIO(result_str)
        blast_records = NCBIXML.parse(result_io)
        
        # 根据保守结构域筛选基因家族成员
        family_members = []
        for blast_record in blast_records:
            print("Alignment:")
            for alignment in blast_record.alignments:
                print(alignment)
                for hsp in alignment.hsps:
                    print(hsp)
                    # 提取SwissProt ID
                    swissprot_id = alignment.title.split("|")[1]

                    # 获取SwissProt记录
                    handle = Bio.ExPASy.get_sprot_raw(swissprot_id)
                    record = Bio.SwissProt.read(handle)

                    # 遍历保守结构域
                    for feature in record.features:
                        if feature.type == "DOMAIN":
                            # 比对参考序列和基因的结构域
                            if hsp.sbjct.find(str(feature.location)) > -1:
                                # 将基因鉴定为参考序列所属基因家族成员
                                family_member = {
                                    'Title': alignment.title,
                                    'Length': alignment.length,
                                    'Alignment': hsp.sbjct
                                }
                                family_members.append(family_member)

        # 输出家族成员列表
        print("Family members:")
        for member in family_members:
            print(member)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

这个代码是我通过不断对chatgpt提问，然后整理得到的，首先第一个定义的函数是执行将我要导入的带鉴定的基因序列文件分成小文件以适应我使用的IDE 的要求(因为文件太大，无法执行后续操作，会报错,如下所示)
，我开始并不知道，故我先执行的是后面的blast对比等一系列流程，后面报错，才又修改了代码。

  raise ValueError(f"Error message from NCBI: {msg}")
ValueError: Error message from NCBI: Your total query length is greater than allowed on the BLAST webserver. You can either reduce the size to 100,000 or less and try again

1
2
3

所以我将大文件分成了几百个小文件，做到这一步我就想着应该行不通了，然后我执行了一下，果不其然，运行的极慢，这肯定不行，这是我直接放弃的主要原因。另外还有的原因是，这样得到的似乎与我前篇文章的基因家族成员鉴定流程不太一致，就算能运行，后面的代码我也不知道该如何实现。
思路一，pass

二、思路二

我想的是既然自己编写blast对比那一系列操作无法实现，那为何不直接用现成的，我想在python中能否实现与生物软件TBtools 的接通，直接通过代码去操纵该软件。
然后，我又去chatgpt了一下，它给出了如下代码

import TBtools as tb
# 创建TBtools对象
tbtools = tb.TBtools()

# 导入基因序列数据
fasta_file = 'gene_sequences.fa'
tbtools.ImportSeq(fasta_file)

# 序列比对
aligned_file = 'gene_alignment.aln'
tbtools.Align('Muscle', 'protein', fasta_file, aligned_file)

# 进行家族成员鉴定
family_file = 'gene_family.txt'
tbtools.AnnotateOrthologs(aligned_file, family_file)

# 输出鉴定结果
results = tbtools.GetAnnotationResults(family_file)

# 打印鉴定结果
for result in results:
    print(result)

# 可以使用其他功能和方法进行可视化或进一步的分析
...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

它是先安装TBtools这个模块，原以为很简单，但我就是被卡住了，我点击安装，它显示安装失败，然后我在终端使用pip安装，它显示已经安装，我就重启了几下pycharm,结果还是不行，折腾了几下就放弃了 （难道是我的IDE有问题，不清楚)。
握不住的沙，不如扬了它。

三、思路三

其实这算是思路二的衍生，既然在python中我不会直接用代码实现基因家族成员鉴定，那就按照步骤一步一步来，毕竟我想实现的只是由机器代替我完成许多重复工作，摸摸鱼而已。我想的是让计算机自己使用鼠标去点击以完成流程，那如何实现呢。
其实这一思路有点无脑，最主要的是获取鼠标所在的界面坐标，然后点击（其实就跟我自己点击一样)。
这样一想确实是简单，但没成想，确是噩梦的开端 。(这一自动化程序代码的准备太费时了，而且还很不稳定，）因为要进入NCBI官网页面，而有时不知是因为宿舍网不好，还是我自己的电脑问题，嗯…,有时会好久都没打开，一直在转圈。出现这种情况，人操作的话，很容易解决，只要关掉再打开就行；而机器的话，反正在我的代码中是不能处理这种情况的，如果遇到这一情况，后面的代码就全运行错误了，(崩溃)

下面是我初始简陋的代码(妄想一套流程全走完)

import os
import time

import pyautogui
import pyperclip
from selenium import webdriver


def get_mouse_click_position():
    pyautogui.click()  # 等待用户点击
    return pyautogui.position()  # 返回鼠标当前的坐标位置


# 设置文件夹路径
folder_path = r"D:\python\PycharmProjects\pythonProject\sys\jd\output"
# 获取文件夹中的所有文件
files = os.listdir(folder_path)
# 获取每个文件的绝对路径
file_paths = [os.path.abspath(os.path.join(folder_path, file)) for file in files]
sorted_list = sorted(file_paths, key=lambda x: int(x.split('_')[2].split('.')[0]))
print(sorted_list)

daibidui_xulie = r"E:\shengxingfenxi\jiyinjiazujianding\Lindera_glauca.pep"
for i in range(len(sorted_list)):
    cankaoxulie = sorted_list[i]  # 参考序列
    # 设置TBtools程序的安装路径
    tbtools_path = "D:\TBtools\TBtools.exe"
    # 设置NCBI官网的URL
    ncbi_url = "https://www.ncbi.nlm.nih.gov/"
    # 创建Chrome浏览器实例
    driver = webdriver.Edge()

    print("===============第一次BLAST比对================")
    pyautogui.press('win')
    pyautogui.typewrite(tbtools_path)
    pyautogui.press('enter')
    time.sleep(20)
    tbtools_window = pyautogui.getWindowsWithTitle("TBtools")[0]
    tbtools_window.activate()

    pyautogui.moveTo(661, 45)
    click_position = get_mouse_click_position()
    print(click_position)
    pyautogui.typewrite("Several Sequences to a Big File [Commonly Used]")
    
    """
    
    
	步骤基本一致，此处就省去了，可自行添加


    """
    time.sleep(150)
    pyautogui.moveTo(962, 541)
    click_position = get_mouse_click_position()
    print(click_position)
    print("****************恭喜，第一次BLAST对比结束*******************")

    print("================删除重复项=================")
    filename = f"E:\shengxingfenxi\jiyinjiazujianding\Lindera_glauca {i} blast"
    gene_names = []
    with open(filename, 'r') as file:
        for line in file:
            line = line.strip()
            if not line.startswith("#"):
                fields = line.split("\t")
                gene_name = fields[1]  # 获取第二列的数据，表示基因名称

                gene_names.append(gene_name)
    gene_names = list(set(gene_names))
    print(len(gene_names))
    os.makedirs(r"E:\1 blast id fa", exist_ok=True)  # 创建文件夹

    new_filename = os.path.join(r"E:\1 blast id fa", f"{i}  blast id fa")

    with open(new_filename, 'w') as f:
        for item in gene_names:
            f.write(item + '\n')
    print("=============得到删除重复项的文件================")
    

    print("==============得到第一次比对后的AA序列==============")
    time.sleep(20)
    pyautogui.moveTo(79, 43)
    click_position = get_mouse_click_position()
    print(click_position)
	"""
	
	此处亦省去


	"""
    print("=============恭喜，已得到AA序列===============")
    driver.switch_to.window(driver.window_handles[-1])

    # 在NCBI官网
    driver.get(ncbi_url)
    print("=========================第二轮BLAST比对=====================")
    time.sleep(40)
    pyautogui.moveTo(859, 658)
    click_position = get_mouse_click_position()
    print(click_position)
	"""
	
	省去

	"""
	
    time.sleep(70)
    pyautogui.moveTo(549, 191)
    click_position = get_mouse_click_position()
    print(click_position)

    print("============恭喜，已得到xml文件=====================")
    driver.switch_to.window(driver.window_handles[0])
    print("==================跳转回TBtools软件，得到AA序列===============")

    time.sleep(10)
    pyautogui.moveTo(14, 8)
    click_position = get_mouse_click_position()
    print(click_position)

  	
  	"""


	省去


  	"""
    time.sleep(25)
    pyautogui.moveTo(962, 547)
    click_position = get_mouse_click_position()
    print(click_position)

    filename = r"E:\BEAPF4VM013-Alignment.xls"
    gene_names = []
    with open(filename, 'r') as file:
        for line in file:
            if line == 0:
                continue
            line = line.strip()
            fields = line.split("\t")
            gene_name = fields[0]
            if gene_name == "Query_def":
                continue
            gene_names.append(gene_name)

    gene_names = list(set(gene_names))
    print(len(gene_names))

    os.makedirs(r"E:\2 blast id fa", exist_ok=True)
    new_filename = os.path.join(r"E:\2 blast id fa", "K00626 2 blast id fa")  # 文件路径和名称

    with open(new_filename, 'w') as f:
        for item in gene_names:
            f.write(item + '\n')
	print("===========得到第二轮blast比对后的基因AA序列==============")
    time.sleep(20)
    click_position = get_mouse_click_position()
    print(click_position)

    """


	省去

    """

    time.sleep(10)
    click_position = get_mouse_click_position()
    print(click_position)

    print("================================================")
    driver.switch_to.window(driver.window_handles[-1])


	print("=============进行保守结构域筛选================")
    # 访问NCBI官网
    driver.get(ncbi_url)

    time.sleep(30)
    pyautogui.moveTo(66, 654)
    click_position = get_mouse_click_position()
    print(click_position)

    """

	省去

    """
    time.sleep(5)
    pyautogui.keyDown("shift")

    time.sleep(10)
    pyautogui.moveTo(277, 443)
    click_position = get_mouse_click_position()
    print(click_position)

    pyautogui.keyUp("shift")

    time.sleep(10)
    pyautogui.moveTo(244, 478)
    click_position = get_mouse_click_position()
    print(click_position)

    time.sleep(30)
    pyautogui.moveTo(356, 488)
    click_position = get_mouse_click_position()
    print(click_position)

    time.sleep(10)
    pyautogui.moveTo(422, 75)
    click_position = get_mouse_click_position()
    print(click_position)
	
    pyautogui.hotkey("ctrl", "c")
    time.sleep(2)
    # 获取剪贴板内容
    clipboard_content = pyperclip.paste()
    # 打印剪贴板内容
    print(clipboard_content)
    # 定义要保存的文件路径
    file_path = 'clipboard.txt'

    # 将剪贴板内容写入文件
    with open(file_path, 'a') as file:
        file.write(clipboard_content + "\n")

    # 打印保存的文件路径
    print(f'已保存剪贴板内容到文件：{file_path}')

    time.sleep(10)  # 等待页面加载完成
    click_position = get_mouse_click_position()
    print(click_position)

    driver.quit()
    pyautogui.hotkey('alt', 'f4')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239

我测试该代码时，是分开来一步一步执行的，按照的是 实验篇——基因家族成员鉴定（一) 一文中的2.5 至2.9步骤。如果有需要可以自己添加或减少所需步骤。
如果顺利的话，最后会得到一个存储最后保守结构域筛选的那个页面的URL的文件。在循环的遍历下会得到许多那个界面的URL。(现在最大的问题是不太顺利，我哭死)
只能说这个代码的可实现性不高，访问NCBI官网太不稳定了，而且不能踏错一步。

四、运用的优化

在思路三中，最主要的问题是不稳定，那我就想是不是因为步骤太多了，那我就不奢求一套流程完成了。我分为几套流程完成，经过我测试这是可行的，保证了进程的稳定性。
以下是一个完成步骤2.5 至 2.7的代码，还算稳定(我先获取了17个参考序列，结果它在循环到第15个时，因为输入的文件路径自动给我由英文以中文打出，结果自然因为文件路径输入错误，TBtools软件报错，之后就进行不下去了，我只能将循环的range(len(sorted_list))改为range(15,17）,让它再跑一次。)
理论上来说，这个代码是能实现的，但是却出现了我上面出现的问题，我也不太清楚为何。

import os
import time

import pyautogui


def get_mouse_click_position():
    pyautogui.click()  # 等待用户点击
    return pyautogui.position()  # 返回鼠标当前的坐标位置


# 设置TBtools程序的安装路径
tbtools_path = "D:\TBtools\TBtools.exe"
# 设置NCBI官网的URL
ncbi_url = "https://www.ncbi.nlm.nih.gov/"

# 设置文件夹路径
folder_path = r"D:\python\PycharmProjects\pythonProject\sys\jd\output"
# 获取文件夹中的所有文件
files = os.listdir(folder_path)
# 获取每个文件的绝对路径
file_paths = [os.path.abspath(os.path.join(folder_path, file)) for file in files]
sorted_list = sorted(file_paths, key=lambda x: int(x.split('_')[2].split('.')[0]))
print(sorted_list)

daibidui_xulie = r"E:\shengxingfenxi\jiyinjiazujianding\Lindera_glauca.pep"
for i in range(len(sorted_list)):
    cankaoxulie = sorted_list[i]  # 参考序列
    if i == 0:
        pyautogui.press('win')
        pyautogui.typewrite(tbtools_path)
        pyautogui.press('enter')
        time.sleep(20)
        tbtools_window = pyautogui.getWindowsWithTitle("TBtools")[0]
        tbtools_window.activate()

        # 打开 Several Sequences to a Big File [Commonly Used]
        time.sleep(5)
        pyautogui.moveTo(158, 42)
        click_position = get_mouse_click_position()
        print(click_position)

       """

		省去
		
       """
    print("===============第一次BLAST比对================")
    time.sleep(5)
    pyautogui.moveTo(320, 70)
    click_position = get_mouse_click_position()
    print(click_position)

	"""
	省去

	"""    pyautogui.typewrite(f"E:\shengxingfenxi\jiyinjiazujianding\Lindera_glauca {i} blast")

    """
    省去
    """
    time.sleep(120)
    pyautogui.moveTo(954, 548)
    click_position = get_mouse_click_position()
    print(click_position)
    print("****************恭喜，第一次BLAST对比结束*******************")

    print("================删除重复项=================")
    filename = f"E:\shengxingfenxi\jiyinjiazujianding\Lindera_glauca {i} blast"
    gene_names = []
    with open(filename, 'r') as file:
        for line in file:
            line = line.strip()
            if not line.startswith("#"):
                fields = line.split("\t")
                gene_name = fields[1]  # 获取第二列的数据，表示基因名称

                gene_names.append(gene_name)
    gene_names = list(set(gene_names))
    print(len(gene_names))
    os.makedirs(r"E:\1 blast id fa", exist_ok=True)  # 创建文件夹

    new_filename = os.path.join(r"E:\1 blast id fa", f"{i}  blast id fa")

    with open(new_filename, 'w') as f:
        for item in gene_names:
            f.write(item + '\n')
    print("=============得到删除重复项的文件================")

    print("==============得到第一次比对后的AA序列==============")

    time.sleep(5)
    pyautogui.moveTo(578, 65)
    click_position = get_mouse_click_position()
    print(click_position)

   """
   省去
   """

    time.sleep(10)
    pyautogui.moveTo(963, 546)
    click_position = get_mouse_click_position()
    print(click_position)
    print("=============恭喜，已得到AA序列===============")

    time.sleep(3)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

通过这一代码，能循环遍历所有存储参考序列的文件，然后分别得到对应的经第一次BLAST比对后的文件。我测试了一下我用了17个参考序列，完成的时间差不多要一个半小时。（还是有点慢，虽然不需要人操作）
其实，这一代码还可以继续优化，因为我发现我要输入待鉴定的序列文件始终是"E:\shengxingfenxi\jiyinjiazujianding\Lindera_glauca.pep",而且只要点击一次初始化就行，我只要在第一次循环时完成前面两步，后续循环完全可以跳过，如此一来又省去了一大部分时间。
至于后续步骤，自然也可分步实现
如此一来，更加清晰明确了每一步骤，之后步骤的代码就不展示了，我么也可以根据自己的需要去设置。反正这一思路是十分简单的。

总结

这几天，一直在探索基因家族成员鉴定代码实现可能性，思路一，思路二都因为各种原因被pass,只有思路三这种无脑且费时的操作能进行下去。这也算是没有办法的办法吧。
其实，在我看来，我失败了，我并没有通过代码快速实现任务。顶多是由机器代替了我，而且不是全部代替，花费的时间甚至更长，而且最后的分步执行程序也不是太稳定，可能依然会出错。
好了，今天的总结就到这里吧。
黄鹤一去不复返，白云千载空悠悠。

-2023-7-19 实验篇