Pysam #

Pysam은 SAM과 BAM 포맷파일을 쉽게 읽고 쓸 수 있는 파이썬 라이브러리 이다. SAM/BAM 파일은 short read sequence 데이터의 alignment 정보를 담고 있는 파일 형식이다. 현재(2017년 11월 기준) 가장 최신 버전은 0.13 이다. 아래와 같이 pip 를 이용하여 설치 할 수 있다. 더 자세한 도움말은 이 곳 에서 제공한다.

$ pip install pysam

bam/sam 파일을 AlignmentFile 객체로 읽어오기 #

pysam의 AlignmentFile() 메소드를 통해 BAM/SAM 파일을 읽을 수 있다. 이 메소드를 통해 읽어온 BAM/SAM 파일은 AlignmentFile 객체로 반환된다.

#!python
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")

AlignmentFile 객체에서 제공하는 fetch(), pileup()등의 메소드를 통해 읽어온 BAM/SAM 파일을 원하는 형태로 읽어올 수 있다. 단 이러한 메소드를 사용하기 위해서는 samtool를 통해 생성한 bam파일의 인덱스파일(.bai)가 같은 경로에 존재해야한다. 다음 명령어를 통헤 BAM 파일의 인덱스를 생성할 수 있다.

$ samtools index ex1.bam ex1.bai

파일쓰기는 "rb"대신 "wb"를 사용한다.

#!python

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")
pairedreads = pysam.AlignmentFile("allpaired.bam", "wb", template=samfile)
for read in samfile.fetch():
     if read.is_paired:
             pairedreads.write(read)

pairedreads.close()
samfile.close()

pysam.AlignmentFile.fetch() 메소드 #

fetch()메소드는 레퍼런스 특정 영역의 맵핑된 리드(read)들을 불러온다.

for read in samfile.fetch('chr1', 100, 120):
    print read
samfile.close()

실행결과는 다음과 같다.

EAS56_57:6:190:289:82       0       99      <<<7<<<;<<<<<<<<8;;<7;4<;<;;;;;94<;     69      CTCAAGGTTGTTGCAAGGGGGTCTATGTGAACAAA     0       192     1
EAS56_57:6:190:289:82       0       99      <<<<<<;<<<<<<<<<<;<<;<<<<;8<6;9;;2;     137     AGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCAC     73      64      1
EAS51_64:3:190:727:308      0       102     <<<<<<<<<<<<<<<<<<<<<<<<<<<::<<<844     99      GGTGCAGAGCCGAGTCACGGGGTTGCCAGCACAGG     99      18      1
...

pysam.AlignmentFile.pileup() 메소드 #

pileup()메소드는 fetch()와 비슷하나 지정한 레퍼런스 범위 내에서 alignment를 불러올때 read단위가 아닌 base단위로 불러온다는 점이 다르다. pileup()은 PileupColumn이라고 하는 특정 위치에 맵핑되는 reads들의 1base의 한 뭉치를 iteration 시켜준다.

#!python
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

예제와 같이 sampile.pileup()을 forloop를 통해 iteration시키면 한번 loop를 돌때마다 PileupColumn객체를 반환한다. 이 객체는 레퍼런스의 single base에 맵핑되는 모든 reads에 대한 정보를 가지고 있다. PileupColumn.pileups 속성을 통해 각각의 read들에 접근할 수 있다. 위 예제의 실행결과는 다음과 같다.

coverage at base 99 = 1
    base in read EAS56_57:6:190:289:82 = A

coverage at base 100 = 1
    base in read EAS56_57:6:190:289:82 = G

coverage at base 101 = 1
    base in read EAS56_57:6:190:289:82 = G

coverage at base 102 = 2
    base in read EAS56_57:6:190:289:82 = G
    base in read EAS51_64:3:190:727:308 = G
...

Suggested Pages #

• 0.025 Java
• 0.025 SAMtools
• 0.025 컴파일
• 0.025 PyPy
• 0.025 Jython
• 0.025 CPython
• 0.025 객체 지향
• 0.025 C#
• 0.025 귀도 반 로섬
• 0.025 BWA
• More suggestions...