1

我有一个给定的染色体编号和位置（chr1 和位置 1599812）。我想使用 python 的 pysam 模块来访问 bam 文件，以仅获取特定区域 chr1 和位置 1599812 的读取数字信息。我尝试过使用pileup()，但它需要一系列位置，而在我的情况下，我只想要一个特定位置而不是这样的范围。

4

2 回答 2

2

我认为这不是pileup()您想要的-根据pysam API，此函数返回“基因组位置的迭代器”，特别是“返回与该区域重叠的'所有'读取。返回的第一个碱基将是第一个碱基第一个读“不一定”是查询中使用的区域的第一个碱基。”

您是说您想获取“读取数信息” - 即该特定位置的读取数，对吗？为此，count_coverage()应该做的工作。在你的情况下，我认为这段代码应该给你你正在寻找的答案：

import pysam

my_bam_file = '/path/to/your/bam_file.bam'
imported = pysam.AlignmentFile(my_bam_file, mode = 'rb')  # 'rb' ~ read bam
coverage = imported.count_coverage(
                  contig = '1',     # Chromosome ID; also might be "chr1" or similar 
                  start = 1599812,
                  stop = 1599813,
                  )
print(coverage)

请注意，这是有效的，因为如pysam API 词汇表中所述，pysam 使用半开区间，因此范围 [1599812, 1599813) 将仅包含一个碱基对。

运行上面的代码会给你这样的东西：

> (array('L', [0]), array('L', [0]), array('L', [0]), array('L', [0]))

它是一个数组元组，分别包含覆盖该基因组位置的读取中的 A、C、G 和 T 碱基数。如果您只是对映射到此特定基因组位置总数的读取数感兴趣，则可以在此元组中求和：

import numpy as np

print(np.sum(coverage))
于 2019-11-18T21:55:43.677 回答
1

如果您设置相同的开始和结束，则堆积将仅引用该特定位置。例如（纯 samtools）：

$ samtools mpileup -r chr1:808957-808957 YourFile.bam
chr1    808957  N   102 READSTRING READQUALITYSTRING

显示 102 个读数，覆盖 1 号染色体的 808957 位置。

于 2015-06-12T11:58:58.123 回答