将SNP映射到CDS的Python脚本

Question

我正在尝试制作一个程序来识别包含光盘的SNP。它从两个文件中填充两个字典，一个包含SNP，另一个包含GFF3文件。从GFF3文件中填充的一个dicts包含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'}}

我的想法是迭代这两个切分，并进行成对的比较，看看SNP位置是否在CDS的间隔中找到。我尝试了下面的代码，但它似乎不起作用。

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))

我必须找到一种方法来评估输出是否正确，但这似乎是合理的。

Answer 1

为了使代码正常工作，需要解决以下几个问题：

由于SNP的位置可以在多条染色体上找到，因此需要考虑到这一点，我用chrom ==基因语句进行了尝试。

您正在尝试将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]): 

Answer 2

也许你的问题过于简单化了。

我们可以使用mysql或UCSC表浏览器下载所有CDS区域，请参阅这个例子文章。

然后使用床身工具与SNP床文件进行重叠合并。