Skip to content

Posts 파이썬으로 메타지놈 기본 프로파일링 따라하기 #
Find similar titles

Structured data

About
pandas
Metagenomics

메타지놈(Metagenomics) NGS 데이터가 있다. 이를 파이썬 pandas로 간단하게 프로파일링(어노테이션)하는 스크립트를 만들어 보았다. 어셈블리(Sequence assembly)한 후 컨티그(contig)로 분석할 수 도 있지만, 생물종별 유전자별 프로파일링한다는 의미로 각 리드(read)들을 직접 BLAST해보고자 한다. 생물정보분석 초심자들에게 도움을 주고자 알기쉽게 정리해 보았다.

Table of Contents

전략 #

FASTQ 데이터를 FASTA로 변환한 뒤, 곧바로 BLAST 구동하고, 생물종 식별번호와 함께 표 형태로 출력한다. pandas로 이 표를 읽어 생물종별, 유전자별로 groupby하고 데이터 현황을 살핀다.

파이썬과 관련모듈 설치 #

파이썬은 최신버전 3.4.2를 쓴다. Python 3도 이젠 충분히 쓸만한 시대다. 더이상 관리되지 않는 파이썬2는 이제 그만. HTSeq이 Python 3을 지원하지 않는 점이 맘에 걸리지만, 내가 만드는 코드가 2.7에도 구동될 수 있도록 프로그래밍하면 그만이다.

OS X 사용자라면 파이썬도 Homebrew로 설치해 쓰는 것이 좋다.

$ brew install python3
$ python3
Python 3.4.3 (default, Feb 25 2015, 21:28:45)
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.56)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>

파이썬 모듈 설치는 pip를 쓰는 것이 제일 좋다. 필요하면 pyvenv(virtualenv)와 조합하여 별도의 독립환경안에서 사용할 수 도 있다. Homebrew에선 Python 3용 pip를 pip3이란 이름으로 관리한다. pip3으로 pandas를 설치해보자.

$ pip3 install pandas
Collecting pandas
  Downloading pandas-0.15.2-cp34-cp34m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (6.8MB)
    100% |################################| 6.8MB 43kB/s
Collecting pytz>=2011k (from pandas)
  Downloading pytz-2014.10-py2.py3-none-any.whl (477kB)
    100% |################################| 479kB 479kB/s
Collecting numpy>=1.7.0 (from pandas)
  Downloading numpy-1.9.2-cp34-cp34m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (3.7MB)
    100% |################################| 3.7MB 88kB/s
Collecting python-dateutil>=2 (from pandas)
  Downloading python_dateutil-2.4.1-py2.py3-none-any.whl (188kB)
    100% |################################| 188kB 420kB/s
Collecting six>=1.5 (from python-dateutil>=2->pandas)
  Downloading six-1.9.0-py2.py3-none-any.whl
Installing collected packages: six, python-dateutil, numpy, pytz, pandas

Successfully installed numpy-1.9.2 pandas-0.15.2 python-dateutil-2.4.1 pytz-2014.10 six-1.9.0

위 명령어 한방에 의존성있는 모듈 NumPy, pytz, python-dateutil, six가 모두 설치된다.

NGS 데이터 받기 #

직접 실험한 NGS 데이터도 있겠지만, 이 예제에선 공개된 데이터를 써보자. 공개된 시퀀싱 원본데이터는 국제적으로 공인된 기관의 SRA 서비스를 통해 받아볼 수 있다. 나의 NGS 데이터도 특별한 절차를 거쳐 등록할 수 있다. NCBI, EBI, DDBJ 각 기관에서 GenBank처럼 공유하여 서비스한다. NCBI에선 특히 SRA Toolkit이란 프로그램을 제공하며, 이를 통해 명령행에서 쉽게 NGS 데이터를 관리할 수 있다.

OS X를 사용한다면 이 역시 Homebrew로 쉽게 설치할 수 있다.

$ brew install sratoolkit

SRA Toolkit을 설치하면 바로 fastq-dump 프로그램을 쓸 수 있다. SRA 식별번호를 입력하여 해당 FASTQ 파일을 내려받는다. 이 예제에선 2010년 공개된 SRX019702: Metagenomics of human microbiome를 이용해보자. SRX019702는 SRR04165~7 4개의 Run 정보로 구성되어 있다. 

$ fastq-dump  SRR041657 --split-files --fasta 60
$ ls
SRR041657_1.fasta    SRR041657_2.fasta

위 명령은 Run 데이터 SRR041657를 FASTA 형식으로 받는다. "--split-files" 옵션은 pair 별로 파일을 분리하라는 의미이다. 이 명령이 수행되면 pair 별로 파일이 생성된다.

BLAST 구동하기 #

입력용 FASTA 파일이 준비됐으니, 이젠 BLAST를 수행해보자. NCBI BLAST+도 Homebrew로 설치된다. 생물정보 분석하기 매우 편해지는 시대. 예전에 매번 NCBI FTP 사이트들어가서 다운받았던 기억이...

$ brew install blast

아, BLAST 검색용 DB는 물론 NCBI FTP 사이트 들어가서 찾아야 한다. 이번 예제에서는 정리된 단백질 데이터베이스인 SWISS-PROT을 써보자. 이 DB는 한참 예전에 UniProt 이름으로 통합됐는데, NCBI FTP에선 여전히 swissprot이란 이름으로 사용된다.

$ wget ftp://ftp.ncbi.nlm.nih.gov/blast/db//swissprot.tar.gz
$ tar zxvf swissprot.tar.gz

BLAST 실행은 다음 명령을 이용한다. 

$ blastx -db swissprot -max_taget_seqs 1 -evalue 0.001 -outfmt '6 std sgi staxids' < SRR041657_1.fasta > SRR041657_1.blout

중요한 옵션들이 있다. "-max_target_seqs 1" 옵션은 베스트 매치(1등)만 쓰겠다는 의미이다. 프로파일링 목적이므로 베스트 매치만 쓰고자 한다. "-outfmt '6 std sgi staxids" 옵션은 6번 출력용 기본 std와 GenBank id인 gi, NCBI Taxon id를 함께 출력하겠다는 의미이다. (자세한 사용법은 "blastn -help" 명령으로 알 수 있다.)

생물종별 유전자별 프로파일링 #

간략한 프로파일 스크립트는 다음과 같다.

argparse 모듈로 명령행 도움말도 바로 출력해 볼 수 있다.

$ python3 metagenome_profiler.py -h
usage: metagenome_profiler.py [-h] fasta

positional arguments:
  fasta       input FASTA file

optional arguments:
  -h, --help  show this help message and exit

FASTA 파일 리드수는 아무래도 꽤 큰 파일이다보니 파이썬에서 직접 읽기 꺼려진다. Stack Overflow 답변들을 보니 "wc -l" subprocess가 가장 나아 보인다. 위 코드의 len_fasta 함수는 "cat fastafile | grep '>' | wc -l" Shell명령을 파이썬으로 구동한 것이다.

(뭔가 검색해서 알아내면 트윗하는 습관...)

생물종별 프로파일링은 다음 한줄이다.

#!python
taxon_counts = subdata.groupby('staxids2').count().sort('qseqid', ascending=False)

pandas가 없었다면 이 부분을 어떻게 만들어야 했을까 싶다. 생물종 id로 groupby하고 갯수를 세어 높은 것 부터 정렬한다가 직관적으로 한줄로 표시된다. pandas 만세!

위 스크립트 파일을 실행한다. 실행인수로 FASTA 파일을 적어주면 FASTA 파일에서 전체 리드수를 확인하고, ".blout" 파일을 읽어 프로파일 데이터를 출력한다. 

$ python3 metagenome_profiler.py SRR041657_1.fasta
Number of reads: 163642
Number of annotations: 1574
Number of taxons: 361
Top 10 taxons
 - id 226186 has 304 reads
 - id 435590 has 174 reads
 - id 272559 has 118 reads
 - id 224308 has 71 reads
 - id 295405 has 57 reads
 - id 83333 has 55 reads
 - id 435591 has 28 reads
 - id 242619 has 26 reads
 - id 3702 has 23 reads
 - id 71421 has 21 reads  
Number of genes: 1190
Top 10 genes
 - id 81740756 has 32 reads
 - id 597502304 has 8 reads
 - id 33301170 has 7 reads
 - id 81315088 has 7 reads
 - id 37538299 has 6 reads
 - id 1709284 has 5 reads
 - id 54035840 has 5 reads
 - id 55976660 has 5 reads
 - id 81444500 has 5 reads
 - id 61227671 has 5 reads

위 데이터는 맛배기로 ".blout"의 일부만 출력한 것이므로 실제 결과와 다르다. 이 예제를 그대로 돌리면 BLAST만도 꽤 오래 걸린다. 큰 서버에서 별도의 분산작업계획하에 돌려 보길 추천한다. 여러 샘플데이터가 있으면 각각 프로파일링 수행 후, 서로간에 얼마나 다른지 비교해 볼 수 있겠다.

그나저나 저사람 장속에 현저하게 많이 발견되는 생물종은 Bacteroides thetaiotaomicron VPI-5482(226186)과 Bacteroides vulgatus ATCC 8482(435590)이며, 많이 발견되는 유전자는 저 1번 생물종의 SUSC_BACTN (Starch-utilization system protein C)이다. 사람마다 저 패턴이 얼마나 많이 다를지 심히 궁금해지는 대목이다.

Suggested Pages #

Other Posts #

comments powered by Disqus
web biohackers.net
0.0.1_20140628_0