2019년 5월 26일 일요일

파이썬으로 BAM 파일 파싱하기, pysam

최근에는 파이썬으로 작업을 최대한 하려고 노력하는 편인데,


  1. 작성이 쉽고,
  2. 기록을 남기기 좋으며,
  3. 유연하기 때문이다.

그런데 여러 작업들을 하다보면 OS에서 실행해야 하는 tool을 파이썬으로 래핑 할 때, 조금 난감할 때가 있다. 

아무래도 가장 좋은 방법은 아무래도 중간 파일 없이 하나의 스크립트로 메모리 상에서 작업을 끝내는 건데, 

외부 tool의 실행 결과를 python의 subprocess 같은 모듈로 받아오면 실행 시간이라던가, msg cache 메모리 부족 같은 생각지도 못한 에러를 마주 할 때가 있다.


다행히도, Samtools view의 경우는 pysam으로 완전 대체 가능!
(simplesam이라는 모듈도 있는데, 이건 걍 samtools 대신 실행시켜주는 모듈)

정말 편하다고 느낀 것은 
보통 parsing을 하다보면 각 라인마다 일정한 구분자(탭, 쉼표 등)로 split하고 각각의 field를 다시 정의하는 과정이 필요한데,
pysam의 경우에는 보다 그 작업이 간편하다. 

다음 코드를 보자, 원하는 값을 정확한 이름(이거 중요하다)으로 받아올 수 있는 이런 편한 모듈을 만들다니.(감격)

1
2
3
4
5
6
7
8
9
10
import pysam
iterable_bam = pysam.AlignmentFile(loc_bam, "rb")
for i, line in enumerate(iterable_bam):
    chrom = line.reference_name
    pos = line.reference_start  #0-based
    strand = line.is_reverse #boolean
    CIGAR = line.cigarstring
    CIGAR_tuples = line.cigartuples
cs

pysam의 AlignmentFile함수는 입력한 경로의 BAM 파일을 iterable 객체로 반환하고, 

순환문을 이용해서 받아온 객체의 각 item은 pysam의 AlignmentSegment 클래스로써,
여러 종류의 함수 및 속성을 지원한다.

특히 마음에 드는것은 위의 예시에도 있는 cigartuples로,

7I8M12N3D와 같은 복잡한 CIGAR string을 (1, 7) (0, 8).. 식으로 
(operation_number, nucleotide_number)로 반환해줘서 낑낑대면서 직접 하나 하나 경우의 수 따져가며 맞출 필요가 없다는 것.

다음 링크에 pysam에 관한 자세한 설명이 있으니 참고