网页视频抓取脚本(如何用python3编写一个文件序列文件中的基因功能)

　　网页视频抓取脚本(如何用python3编写一个文件序列文件中的基因功能)

　　操作说明

　　在基因组分析中，我们经常有这样的需求，需要从fasta文件中提取一些序列。有时这些序列是完整的序列，有时它们只是原创 fasta 文件中序列的一部分。尤其是数据量大的时候，用肉眼来选择序列会很困难，那么我们可以通过简单的编程来实现。

　　例如，一个物种的全基因组序列（0-refer/Bacillus_subtilis.str168.fasta）及其基因组gff注释文件（0-refer/Bacillus_subtilis.str168.gff）。这里假设我们对这个物种进行研究，通过gff注释文件中的基因功能描述字段定位一些特定的基因，加上相关数据的访问等。接下来，我们期望在全基因组序列中找到并提取这些基因fasta 文件根据 gff 文件中这些基因的位置描述，得到一个新的 fasta 文件，该文件只收录目标基因序列。

　　请使用python3编写一个可以做到这一点的脚本。

　　例子

　　示例脚本如下（见网盘附件“seq_select1.py”）。

　　为了达到上述目的，我们首先需要准备一个txt文件（以下简称list文件，示例list.txt可以在网盘附件中找到），根据记录的基因位置信息gff 文件，填写类似如下内容（列和列之间用制表符分隔）。

　　1.png

　　#将以下内容保存到list.txt

　　基因46 NC_000964.3 42917 43660 +

　　NP_387934.1 NC_000964.3 59504 60070 +

　　yfmC NC_000964.3 825787 826734 -

　　cds821 NC_000964.3 885844 886173 -

　　在第 1 列中，为要获得的新序列命名；

　　第二列，要获取的序列的原创序列ID；

　　第3列，原创序列中要获取的序列的起始位置；

　　第4列，原创序列中要得到的序列的终止位置；

　　在第 5 列中，要获得的序列位于原创序列的正 (+) 或负 (-) 链上。

　　然后根据输入文件编辑py脚本，即输入fasta文件和list文件中记录要获取的序列位置的内容。

　　打开fasta文件“Bacillus_subtilis.scaffolds.fasta”，使用循环逐行读取序列id和碱基序列，将每个序列的所有碱基组合成一个字符串；基本序列存储为字典（字典样式 {'id':'base'}）。

　　打开列表文件“list.txt”，读取内容，存入字典。字典的键是列表文件中第一列的内容；字典的值是列表文件第2-5列的内容，除以tab得到一个列表，收录4个字符分别代表列表文件第2-5列的信息）。

　　最后根据read list文件中的序列位置信息，从read基因组中截取目标基因序列。由于某些基因序列可能位于基因组的负链中，因此需要取反向互补序列。因此，首先定义一个函数rev()，以便在后续调用中获取反向互补序列。在输出序列名称时，还可以选择是否将序列的位置信息一起输出（name_detail = True/False）。

　　#!/usr/bin/env python3

　　# -*- 编码：utf-8 -*-

　　#初始通过命令

　　input_file = '枯草芽孢杆菌.str168.fasta'

　　list_file = 'list.txt'

　　output_file = 'gene.fasta'

　　name_detail = 真

　　##读取文件

　　#读取基因组序列

　　seq_file = {}

　　使用 open(input_file, 'r') 作为 input_fasta：

　　对于 input_fasta 中的行：

　　line = line.strip()

　　如果行 [0] == '>'：

　　seq_id = line.split()[0]

　　seq_file[seq_id] = ''

　　别的：

　　seq_file[seq_id] += 行

　　input_fasta.close()

　　#读取列表文件

　　list_dict = {}

　　使用 open(list_file, 'r') 作为 list_table：

　　对于 list_table 中的行：

　　如果 line.strip():

　　line = line.strip().split('\t')

　　list_dict[line[0]] = [line[1], int(line[2]) - 1, int(line[3]), line[4]]

　　list_table.close()

　　##截取序列并输出

　　#定义一个截取逆补的函数

　　def rev(seq):

　　base_trans = {'A':'T', 'C':'G', 'T':'A', 'G':'C', 'N':'N', 'a':'t' , 'c':'g', 't':'a', 'g':'c', 'n':'n'}

　　rev_seq = 列表（反转（seq））

　　rev_seq_list = [base_trans[k] for k in rev_seq]

　　rev_seq = ''.join(rev_seq_list)

　　返回（rev_seq）

　　#截断序列并输出

　　output_fasta = 打开（输出文件，'w'）

　　对于 list_dict.items() 中的键、值：

　　如果名称_详细信息：

　　print('>' + key, '[' + value[0], value[1] + 1, value[2], value[3] + ']', file = output_fasta)

　　别的：

　　打印（'>' + 键，文件 = output_fasta）

　　seq = seq_file['>' + value[0]][value[1]:value[2]]

　　如果值 [3] == '+'：

　　打印（序列，文件= output_fasta）

　　elif 值 [3] == '-'：

　　seq = rev(seq)

　　打印（序列，文件= output_fasta）

　　output_fasta.close()

　　编辑好脚本后，运行并输出一个新的fasta文件“gene.fasta”，其序列就是我们要得到的目标基因序列。

　　2.png

　　延期：

　　网盘附件“seq_select.py”是一个python3脚本，增加了命令行传输线，可以直接在shell中对目标文件进行I/O处理。该脚本可以指定一个输入fasta序列文件和一个记录要提取的序列位置的列表文件，新的输出fasta文件就是提取的序列。

　　#!/usr/bin/env python3

　　# -*- 编码：utf-8 -*-

　　#导入模块，初始传递命令、变量等

　　导入参数解析

　　parser = argparse.ArgumentParser(description = '\n该脚本用于截取基因组中特定位置的序列，并额外输入截取序列信息的列表文件', add_help = False, usage = '\npython3 seq_select.py -i [ input.fasta] -o [output.fasta] -l [list]\npython3 seq_select.py --input [input.fasta] --output [output.fasta] --list [list]')

　　required = parser.add_argument_group('必需选项')

　　可选 = parser.add_argument_group('可选')

　　required.add_argument('-i', '--input', metavar = '[input.fasta]', help = '输入文件, fasta 格式', required = True)

　　required.add_argument('-o', '--output', metavar = '[output.fasta]', help = '输出文件，fasta 格式', required = True)

　　required.add_argument('-l', '--list', metavar = '[list]', help = 'record "新序列名称/原序列的序列ID/序列起始位置/序列结束位置/正链( +) 或负链 (-)"，以制表符分隔，必填 = True)

　　optional.add_argument('--detail', action = 'store_true', help = '如果该参数存在，则在输出fasta的每个sequence id中显示序列在原创fasta中的位置信息', required = False)

　　optional.add_argument('-h', '--help', action = 'help', help = '帮助信息')

　　args = parser.parse_args()

　　##读取文件

　　#读取基因组序列

　　seq_file = {}

　　使用 open(args.input, 'r') 作为 input_fasta：

　　对于 input_fasta 中的行：

　　line = line.strip()

　　如果行 [0] == '>'：

　　seq_id = line.split()[0]

　　seq_file[seq_id] = ''

　　别的：

　　seq_file[seq_id] += 行

　　input_fasta.close()

　　#读取列表文件

　　list_dict = {}

　　使用 open(args.list, 'r') 作为 list_file：

　　对于 list_file 中的行：

　　如果 line.strip():

　　line = line.strip().split('\t')

　　list_dict[line[0]] = [line[1], int(line[2]) - 1, int(line[3]), line[4]]

　　list_file.close()

　　##截取序列并输出

　　#定义一个截取逆补的函数

　　def rev(seq):

　　base_trans = {'A':'T', 'C':'G', 'T':'A', 'G':'C', 'a':'t', 'c':'g' , 't':'a', 'g':'c'}

　　rev_seq = 列表（反转（seq））

　　rev_seq_list = [base_trans[k] for k in rev_seq]

　　rev_seq = ''.join(rev_seq_list)

　　返回（rev_seq）

　　#截断序列并输出

　　output_fasta = open(args.output, 'w')

　　对于 list_dict.items() 中的键、值：

　　如果 args.detail:

　　print('>' + key, '[' + value[0], value[1] + 1, value[2], value[3] + ']', file = output_fasta)

　　别的：

　　打印（'>' + 键，文件 = output_fasta）

　　seq = seq_file['>' + value[0]][value[1]:value[2]]

　　如果值 [3] == '+'：

　　打印（序列，文件= output_fasta）

　　elif 值 [3] == '-'：

　　seq = rev(seq)

　　打印（序列，文件= output_fasta）

　　output_fasta.close()

　　使用上面示例中的测试文件，脚本运行如下。

　　#python3 seq_select.py -h

　　python3 seq_select.py -i Bacillus_subtilis.str168.fasta -l list.txt -o gene.fasta --detail

　　3.png