파이썬
라이브러리
HTSeq
#
Find similar titles
- 최초 작성자
-
최근 업데이트
khcho@insilicogen.com
Structured data
- Category
- Programming
Table of Contents
HTSeq #
대용량 NGS데이터를 이용하여 분석할 때 많은 Samtools 둥 다양한 방법을 이용하고 있지만 개발자의 입장에서 좀 더 파이썬 친화적인 분석을 하고자 한다면 HTSeq 라이버러리를 이용하는 것이 한 방법일 수 있다. HTSeq은
등 다양한 분석 작업을 할 수 있다. 파이썬을 이용하여 분석 스크립트를 작성하거나 HTS 데이터를 다룬다면 그 궁합은 두말하면 잔소리다.
HTSeq 설치 #
먼저 사용을 위해서는 먼저 관련 팩키지를 설치 한 후 사용할 수 있다. 먼저 그 전에 이 팩키지를 사용할 수 있도록 관련 프로그램을 먼저 인스톨해야 한다. 관련 프로그램을 설치하는 것은 간단하다. 만약 데비안 계열을 사용한다면 다음과 같이 간단한 명령어로 빌드 관련 툴, 파이썬 및 관련 팩키지를 설치하면 된다.
sudo apt-get install build-essential python2.7-dev python-numpy python- matplotlib
만약 페도라, CentOS등의 레드햇 계열의 리눅스를 이용하고 있다면 마찬가지로 아래의 명령어를 이용하여 간단히 설치할 수 있다.
sudo yum groupinstall "Development Tools"
sudo yum install python-devel numpy python-matplotlib
HTSeq 사용 #
NGS 데이터 QC #
이제 HTSeq을 사용할 준비가 되었다. 사용을 위해서 먼저 관련 소스를 다운받아야 한다. 다음의 주소에서 https://pypi.python.org/packages/source/H/HTSeq/HTSeq-0.6.1.tar.gz#md5=b7f4f38a9f4278b9b7f948d1efbc1f05 를 다운받을 수 있다.
다운후에 압축을 풀고 디렉토리로 이동하여 설치 준비를 한다. 빌드 및 설치 과정은 다음과 같다.
python setup.py build
sudo python setup.py install
시스템에 따라서 혹은 계정에 따라서 root 권한이 필요할 수도 있다. 설치 후 관련 패키지를 호출하면 이제 대용량 NGS데이터를 다룰 준비가 되었다.
>>> import HTSeq
>>> fastq_file = HTSeq.FastqReader( "yeast_RNASeq_excerpt_sequence.txt", "solexa" )
먼저 파이썬을 실행 후 파일을 열어 관련 정보를 얻을 수 있다. 두개의 인자가 필요하다. 하나는 파일이름며 두번째는 퀄리티를 위한 인코딩영역을 지정한다. 여기에서는 솔렉사의 파이프라인에 의해서 생성된 파일을 지정한다. 이제 파일을 읽은 정보를 확인해보면 파일의 내용을 간략히 확인할 수 있다.
>>> fastq_file
<FastqReader object, connected to file name 'yeast_RNASeq_excerpt_sequence.txt'>
대용량의 데이터이기 때문에 파일의 10번째 항목까지 불러와서 정보를 확인한다. itertools 팩키지는 관련 파일을 읽어서 순환할 수 있는 구조를 만들어준다. for문을 이용하여 데이터를 읽어 출력해보면 다음과 같은 내용을 확인할 수 있다.
>>> import itertools
>>> for read in itertools.islice( fastq_file, 10 ):
... print read
CTTACGTTTTCTGTATCAATACTCGATTTATCATCT
AATTGGTTTCCCCGCCGAGACCGTACACTACCAGCC
TTTGGACTTGATTGTTGACGCTATCAAGGCTGCTGG
ATCTCATATACAATGTCTATCCCAGAAACTCAAAAA
AAAGTTCGAATTAGGCCGTCAACCAGCCAACACCAA
GGAGCAAATTGCCAACAAGGAAAGGCAATATAACGA
AGACAAGCTGCTGCTTCTGTTGTTCCATCTGCTTCC
AAGAGGTTTGAGATCTTTGACCACCGTCTGGGCTGA
GTCATCACTATCAGAGAAGGTAGAACATTGGAAGAT
ACTTTTAAAGATTGGCCAAGAATTGGGGATTGAAGA
위의 예제에서 마지막에 저장된 read의 정보를 확인해보면 SequenceWithQualities라는 객체 정보를 확인한다. 이는 하나의 리드 정보를 객체로 관리하기 때문에 리드명, 리드서열, 리드 퀄리티 정보를 각각 .xxx로 확인할 수 있다.
>>> read
<SequenceWithQualities object 'HWI-EAS225:1:10:1284:142#0/1'>
>>> read.name
'HWI-EAS225:1:10:1284:142#0/1'
>>> read.seq
'ACTTTTAAAGATTGGCCAAGAATTGGGGATTGAAGA'
>>> read.qual
array([33, 33, 33, 33, 33, 33, 29, 27, 29, 32, 29, 30, 30, 21, 22, 25, 25,25, 23, 28, 24, 24, 29, 29, 29, 25, 28, 24, 24, 26, 25, 25, 24, 24,
24, 24])
하나의 예로 읽어들인 리드들의 평균 퀄리티 값을 계산하기 위해 파이썬의 수학 패키지인 numpy를 이용하여 다음과 같이 간단히 리드의 길이를 계산한다. 다음 리드의 갯수로 나누어 주면 된다.
>>> import numpy
>>> len( read )
36
>>> qualsum = numpy.zeros( len(read), numpy.int )
이러한 정보를 이용하여 모든 리드들의 정보를 읽은 후 퀄리티 값과 리드의 갯수를 카운팅하여 평균을 내어보면 그 평균값을 구할 수가 있다. 이때 이러한 정보를 그래프 형태의 플롯으로 보고자 한다면 pyplot를 이용하며 그림으로 확인할 수도 있다.
>>> from matplotlib import pyplot
>>> pyplot.plot( qualsum / nreads )
[<matplotlib.lines.Line2D object at 0x...>]
>>> pyplot.show()
지금까지 우리는 아주 간단한 예제를 통하여 퀄리티 값을 추출하고 이를 이용하여 퀄리티 평균을 구하여 보았다. 만약 FastQ파일이 아니라 이미 분석한 SAM, BAM형식의 파일을 가지고 있다면 SAM_Reader, BAM_Reader를 이용하여 마찬가지로 동일한 결과를 얻을 수 있다.
RNASeq 분석 #
다음에는 HTSeq을 이용한 RNASeq분석을 다뤄보고자 한다.
