python提取目标片段

将每份材料组装的 contig 长度大于 500 bp 的提取出来，

1. 首先，打开组装结果文件，该文件一般为fasta格式，其中包含了所有的contig序列。
2. 读取fasta文件，将每个contig的序列和长度存储到一个字典中。
3. 遍历字典中的每个contig，筛选出长度大于500 bp的contig，并将其序列写入一个新的fasta文件中。
   以下是一个Python代码示例，用于实现上述步骤：

from Bio import SeqIO

# 输入文件路径
input_file = "assembly_result.fasta"
# 输出文件路径
output_file = "contigs_gt_500bp.fasta"

# 存储contig序列和长度的字典
contigs = {}

# 读取fasta文件，将contig序列和长度存储到字典中
for record in SeqIO.parse(input_file, "fasta"):
    contig_seq = str(record.seq)
    contig_length = len(contig_seq)
    contigs[record.id] = (contig_seq, contig_length)

# 筛选出长度大于500 bp的contig，并将其写入新的fasta文件中
selected_contigs = {contig_id: contig_info for contig_id, contig_info in contigs.items() if contig_info[1] > 500}

with open(output_file, "w") as f:
    for contig_id, contig_info in selected_contigs.items():
        f.write(f">{contig_id}\n{contig_info[0]}\n")

对于单条 contig，如果连续有 300 bp 比对上的片段，并且比对的一致性大于 90%的区域（参数为 -I 90 -l 300 -q），我们定义为比对上的片段。如果有连续 500 bp 为未比对上的序列，将这些序列提取出来，和未比对的片段合并。

要提取连续500 bp为未比对上的序列，并与未比对的片段合并，可以使用samtools软件进行操作。以下是一个示例操作流程：

1. 首先，运行samtools的view命令，将比对结果（BAM格式）转换为SAM格式：

samtools view -h alignment.bam > alignment.sam

2. 使用samtools的calmd命令对SAM文件进行碱基比对序列的修改：

samtools calmd -b alignment.sam reference.fasta > alignment_calmd.bam

其中，reference.fasta是参考基因组的FASTA文件。

3. 运行samtools的view命令，将修改后的比对结果转换回SAM格式：

samtools view -h alignment_calmd.bam > alignment_calmd.sam

4. 使用samtools的bedcov命令计算每个contig上的覆盖度：

samtools bedcov reference.bed alignment_calmd.sam > coverage.txt

其中，reference.bed是参考基因组的BED文件。

5. 读取coverage.txt文件，提取连续300 bp比对上的片段，并且比对的一致性大于90%：

import re

# 定义参数
min_length = 300
min_identity = 0.9

# 存储提取的比对上的片段
aligned_fragments = []

# 读取coverage.txt文件
with open("coverage.txt", "r") as f:
    for line in f:
        contig, start, end, _, coverage = line.strip().split("\t")
        coverage = int(coverage)
        length = int(end) - int(start)

        # 判断是否满足条件
        if length >= min_length and coverage / length >= min_identity:
            aligned_fragments.append((contig, int(start), int(end)))

# 对提取的比对上的片段按照位置进行排序
aligned_fragments.sort(key=lambda x: (x[0], x[1]))

# 提取未比对上的序列
unaligned_sequences = []
prev_contig = ""
prev_end = 0

for fragment in aligned_fragments:
    contig, start, end = fragment

    if contig != prev_contig or start > prev_end + 1:
        # 提取未比对上的片段
        unaligned_start = prev_end + 1
        unaligned_end = start - 1
        unaligned_sequence = (prev_contig, unaligned_start, unaligned_end)
        unaligned_sequences.append(unaligned_sequence)

    prev_contig = contig
    prev_end = end

# 处理最后一个片段
contig_length = end # 假设最后一个片段是contig的最后一段
if contig != prev_contig or end < contig_length:
    # 提取未比对上的片段
    unaligned_start = prev_end + 1
    unaligned_end = contig_length
    unaligned_sequence = (prev_contig, unaligned_start, unaligned_end)
    unaligned_sequences.append(unaligned_sequence)

# 合并未比对的片段
merged_unaligned_sequence = []
prev_contig = ""
prev_end = 0

for sequence in unaligned_sequences:
    contig, start, end = sequence

    if contig != prev_contig or start > prev_end + 1:
        # 开始新的未比对片段
        merged_unaligned_sequence.append((contig, start, end))
    else:
        # 合并连续的未比对片段
        merged_unaligned_sequence[-1] = (contig, merged_unaligned_sequence[-1][1], end)

    prev_contig = contig
    prev_end = end

# 打印合并后的未比对片段
for sequence in merged_unaligned_sequence:
    contig, start, end = sequence
    print(f"{contig}\t{start}\t{end}")

在上述示例代码中，需要将alignment.bam、reference.fasta和reference.bed的路径修改为实际的文件路径。同时，也可以根据实际需要调整min_length和min_identity参数的值。

示例代码首先读取coverage.txt文件，提取连续300 bp比对上的片段，并且比对的一致性大于90%的片段。然后，按照位置对这些片段进行排序，并提取出未比对上的序列。最后，将连续的未比对片段合并，并打印合并后的结果。