我正在尝试制作一个程序,该程序将标识包含CDSSNP。它从两个文件填充两个字典,一个包含SNP,另一个包含GFF3文件。从GFF3文件填充的一则字典包含CDS名称以及它们作为元组的位置。

一个例子:

cds_pos = {'PRELSG_0019500_6': ('2091320', '2092988'), 'PRELSG_1338600_3': ('1542760','1542853'), 'PRELSG_0013000_1': ('1275531', '1275568')}


另一个字典包含染色体和SNP的位置,例如。

chrom_pos = {'PRELSG_13_v1': {'272093', '964287', '844454', '65770', '336211', '36660'}, 'PRELSG_12_v1': {'1270177', '1368630'}}


我的想法是遍历两个字典并进行成对比较,看看是否在CDS的间隔中找到了SNP位置。我尝试了下面的代码,但似乎没有用。

for chrom, snp_pos in chrom_pos.items():
    for cds, pos in cds_pos.items():
            if pos[0] <= str(snp_pos) <= pos[1]:
                print(cds)
                print(snp_pos)


对于我发现不起作用的事情。首先,没有什么满足区间if语句。其次,由于可以在多个染色体上找到SNP的位置,因此需要考虑这一点,我尝试使用chrom ==基因声明。但这似乎不起作用。

对于必须进行的任何想法和评论,我们将非常高兴。谢谢

编辑:

到目前为止,我的脚本看起来像:

cds_snp = defaultdict(set)
for chrom, snp_pos in chrom_pos.items():
        for cds, pos in cds_pos.items():
            ed_chrom = chrom[:9]
            ed_cds = cds[:9]
            if ed_cds == ed_chrom:
                for i in snp_pos: # Iterate through the set of snp positions
                    if int(pos[0]) <= int(i) <= int(pos[1]):
                        cds_snp[cds].add(i)

for i,j in sorted(cds_snp.items()):
    print(i)
    print('\n'.join(j))


我必须找到一种方法来评估输出是否正确,但似乎是可行的。

最佳答案

为了使此代码正常工作,需要解决一些问题:


  由于可以在多个染色体上找到SNP的位置,因此需要考虑这一点,我尝试使用chrom == gene语句。


您正在尝试将SNP与CDS中的位置链接。但是,您似乎没有CDS位置的染色体(只是名称)。在您发布的示例中,您的CDS名称为'PRELSG_0019500_6', 'PRELSG_1338600_3', 'PRELSG_0013000_1',染色体名称为'PRELSG_10_v1', 'PRELSG_12_v1'。这些都不是匹配的情况,老实说,它们看起来像是具有不同的格式,并且您永远不会遇到chrom == gene

CDS名称中是否存在一些识别信息,这些信息可以告诉您它位于哪个染色体上?如果存在,则可以从染色体名称中提取染色体编号(即从“ PRELSG_12_v1”中提取12个染色体编号),并将其与从CDS名称中提取的染色体编号进行比较。


  没有什么可以满足区间if语句


通过这种方式,我假设您的意思是if pos[0] <= str(snp_pos) <= pos[1]:

有一个简单的解释说明为什么它不起作用。首次从字典中提取特征时:

for chrom, snp_pos in chrom_pos.items():
    for gene, pos in gene_pos.items():


您提取snp_pos。这不是个人职位,而是一组职位。要遍历各个位置,您可以添加另一个循环:

for snpPos in snp_pos:


最后,为了正确起见,最好将您的interval语句保持为整数。即写为

if int(pos[0]) <= int(snpPos) <= int(pos[1]):

关于python - 将SNP映射到CDS的Python脚本,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/43325602/

10-12 23:12