Python 将 SNP 映射到 CDS 的脚本

Question

我正在尝试制作一个程序来识别 CDS containing SNP. It populates two dictionaries from two files, one containing the SNP and the other a 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'}}

我的想法是遍历两个dict，进行两两比较，看是否在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)

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

如果有任何想法和评论必须继续，我们将很高兴。谢谢

编辑：

到目前为止，我的脚本如下所示：

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

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

since the position of the SNP can be found on several chromosomes, this needs to be taken into account, which I tried with the chrom == gene statement.

您正在尝试 link 在 CDS 中具有位置的 SNP。但是，您似乎没有 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 名称中提取的染色体编号进行比较。

nothing satisfy the interval if statement

据此，我假设您指的是行 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:

最后，为了正确起见，最好将间隔语句保留为整数。即写成

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

Answer 2

可能把你的问题简单化了。

我们可以使用 mysql 或 UCSC table 浏览器下载所有 CDS 区域，参见 this example post。

然后使用 bedtools 对 SNP 床文件进行 merge on overlap。

Python 将 SNP 映射到 CDS 的脚本

Python script to map SNP to CDS

python

bioinformatics