Help
파이썬
라이브러리
pySam
#
Find similar titles
- (rev. 6)
- ㅎnㅂrㄹrㄱi
Structured data
- Category
- Programming
Table of Contents
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
...
