python - 使用 Python 从 BAM 文件和 vcf 文件中提取具有不同位置的读取和配对（PYSAM 和 PYSAM，也可以使用 bamnostics）

Question

我正在使用 Python 并使用 PYSAM 和 PyVCF 库进行数据处理。嗨，假设我有一个 bam 文件和一个包含变体调用结果的 vcf 文件。我想提取其中有变体的读取和配对。假设 mate 中的 read 和 variant 有一个变体，我想提取这些 reads 并丢弃所有没有变体位置的 reads 或者它们的 mate 没有变体（即使 reads 有变体但 mate 没有变体）。之后，我想计算从结果中生成的其他 read_mate 模式的计数。

完整问题：我们必须使用 VCF 文件搜索 bam 文件，如果 vcf 文件上的位置与 bam 文件中的读取匹配，则检查是否在同一位置上是否存在与 vcf 相同的 alt，如果两者相同则考虑一下并保存并继续前进。在处理完基因的所有数据之后，（逐个基因），我们将保存所有那些在其中具有 vcf 位置的读取和配对。具有vcf位置的读取和配对将在一个变量中逐个基因保存，信息应保存如下：最初是一个带有基因名称的输出文件，contig name : <position,alt base>, <Position,alt base> 等等每次阅读和配对。下一步在保存那些reads和只有vcf位置的mates的信息后，注意：（丢弃那些mates没有vcf位置的reads），下一步是计算read support，

下一步：合并基本模式和噪音：例如，我们有以下模式 pos alt pos alt pos alt post alt pos alt pos alt pos alt（读取，配对模式）：55 G 58 C 63 G 75 T 87 A 87 T 95 A（读数，配对模式）：C 63 G 75 T 87 A 87 T 95 A（读数，配对模式）：55 G 58 C 63 G 75
（读数，配对模式）：55 G 58 C 63 G 75 T 87 A 87 T 95 A

现在所有这些模式应该合并成一个更大的模式：pos alt pos alt pos alt post alt pos alt pos alt（过滤模式）：55 G 58 C 63 G 75 T 87 A 87 T 95 A

底线：我只想获得这些读取和具有 vcf 位置的配对，并为每个读取及其具有 vcf 位置的配对制作一个模式。并丢弃所有其他没有 vcf 位置的人。

我在映射中遇到了一些问题，不知道如何使用雪茄串来映射我的数据并在与 vcf 位置匹配后将 DNA 序列与其进行比较以检查它是否具有碱基。任何帮助将不胜感激。

我尝试了以下代码：

import pysam
import vcf
#Main logics:
j=0

for ind,key in gencode_genes.iterrows():
#     print(key[0])
    i=0
#below is signature of bamnostic.fetch() function     
# def fetch(self, contig=None, start=None, stop=None, region=None,tid=None, until_eof=False, multiple_iterators=False,reference=None, end=None):
#     for read in samfile.fetch(key[0],None,None,None,None,True):
    for read in samfile.fetch(key['seqname'],key['start'],key['end'],None,None,True):
        
#         if read's position and mate's position is same than it skips
#         if(read.pos == read.next_reference_start):
#             continue
        b=vcf_df.query('POS>={} and POS<={}'.format(read.reference_start,read.reference_end))
        a=b.sort_values(by = 'POS')
#         if any position in vcf file found which matches with read names start and end position. then it proceeeds
        if(not a.empty):
            if read.read_name not in result.keys():
                
                result[read.read_name]=[]

依此类推，长代码但不是所需的输出。谢谢

标签： pythonbioinformaticspysambam