我有给定的染色体编号和位置(chr1 和位置 1599812)。我想使用 python 的 pysam 模块访问 bam 文件以获取仅特定区域 chr1 和位置 1599812 的读取数字信息。我尝试使用 pileup()
但它需要一系列位置而就我而言,我只想要一个特定的位置,而不是一系列这样的位置。
最佳答案
我认为 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)
Note that this works because, as noted in the pysam API glossary, pysam uses half-open intervals, so the range [1599812, 1599813) will include exactly one base-pair.
运行上面的代码会给你这样的结果:
> (array('L', [0]), array('L', [0]), array('L', [0]), array('L', [0]))
这是一个数组元组,分别包含覆盖该基因组位置的读数中 A、C、G 和 T 碱基的数量。如果您只是对映射到该特定基因组位置总数的读取数量感兴趣,则可以对该元组求和:
import numpy as np
print(np.sum(coverage))
关于python - 使用 Pysam 访问特定位置的 Bam 文件,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/30697271/