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) には正確に 1 つの塩基対が含まれます。

上記のコードを実行すると、次のようになります。

> (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 に答える