python :: pysam 모듈을 이용한 bam, vcf 파일 다루기

python에서 bam이나 vcf 파일을 직접 다루기에는 정보도 많고 한줄한줄 코딩하기가 어렵다.

pysam이란 모듈을 이용하면 좀 더 쉽고 빠르게 bam 파일이나 vcf 파일을 다룰 수 있기 때문에 pysam 모듈 사용법에 대해 알아보려고한다.

 

1. pysam이란?

pysam은 SAM/BAM, VCF 형식의 파일을 읽고 조작하는 파이썬 모듈이다.

이 모듈은 samtools의 기능을 간단하게 조작도 가능하며 tabix에 대한 인터페이스도 포함하고 있다.

 

2. pysam 설치하기

pysam을 설치하는 방법은 두가지가 있는데, conda 를 이용해서 설치하거나 직접 pip install 명령어로 설치하면된다.

 

<conda로 설치하기>

conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
conda install pysam

 

 

<pip 으로 설치하기>

pip install pysam

 

 

3. pysam을 이용한 bam파일 다루기

pysam document에 잘 나와있긴 하나 기록용으로 vcf 파일 다루는법에 대해 설명하겠다.

document는 아래 링크를 참고하면 된다.

https://pysam.readthedocs.io/en/latest/api.html

 

먼저 pysam을 import 해온다.

사용하려는 bam파일의 경로를 ex1.bam 위치에 넣어주고 bam파일이므로 rb 즉 read binary 옵션을 써준다.

 

• r: 파일을 읽기 전용으로 파일은 연다.
• b: BAM 형식의 바이너리 파일을 처리.

 

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

 

bam파일에서 특정 범위의 read를 가져와 출력하기

 

chr1 염색체의 100번부터 120번 염기까지의 범위에 해당하는 read를 가져오고싶을때 아래와 같이 사용해준다.

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

samfile.close()

 

bam파일에서 paired read만 저장하기

bam파일에서 paired read만 가져와서 사용하고 싶을때 아래와 같은 코드를 사용한다.

allpaired.bam라는 bam파일이 최종 output이고 해당 bam파일에 ex1.bam에서의 paired read를 저장한다.

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()

 

4. pysam을 이용한 vcf 파일 다루기

pysam을 이용해서 vcf 파일을 가져오는 방법은 매우 간단하다.

VariantFile을 이용해서 불러오면 헤더나 컬럼명을 따로 정해주지 않아도 vcf를 읽어올 수 있다.

test.vcf라는 vcf파일을 vcf_in 이라는 변수로 불러와보겠다.

VCF 파일 읽어오기

from pysam import VariantFile

vcf_in = VariantFile("test.vcf")  # auto-detect input format
vcf_out = VariantFile('-', 'w', header=vcf_in.header)

for rec in vcf_in.fetch('chr1', 100000, 200000):
    vcf_out.write(rec) ##chr1:100000-200000 범위의 변이가 저장됨

 

VCF 파일에서 특정 컬럼의 정보 가져오기

 

vcf파일에서 chromosome, position, ref, alt 정보를 확인하고 싶은 경우rec.chrom 처럼 . 뒤에 컬럼명을 써주면 된다.

for rec in vcf.fetch():
  # print(rec.header)
  chrom=rec.chrom
  pos=str(rec.pos)
  ref=rec.ref
  alt=rec.alts[0]
  # print(rec)
  if len(ref) > 1 or len(alt) > 1 :
    print(ref, alt)

 

기존에 VCF파일을 읽어오기 위해 for문으로 한줄씩 코딩을 했다면 모듈을 이용하면 좀 더 쉽게 파일을 다룰 수 있다.