[유전자분석] samtools 유전자 데이터 컨트롤하기

지난 포스팅에 유전자 데이터와, samtools 을 설치하였다.

이제 samtools 를 사용해보자.

 

실습에 사용할 데이터는 아래 사이트에서 제공하고 있다.

Index of /goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq

 

샘플 파일 다운로드

위 사이트에서 아무거나 다운해보자.

$ wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam

 

BAM -> SAM 변환

다운받은 바이너리형태인 bam 파일을 아스키코드로 된 sam파일로 변환해보자.

samtools view {파일명}.bam > output.sam

$ samtools view wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam > output.sam

 

SAM 파일 보기

간단히 head 명령어로 파일을 읽어보았다.

(아스키코드기 때문에 사람이 읽을 수 있는 형태)

$ head output.sam
SOLEXA-1GA-2_2_FC20EMB:5:251:979:328    0       chr1    10145   25      36M     *       0       0       AACCCCTAACCCTAACCCTAACCCTAACCCTAAACT    hhhhHcWhhHTghcKA_ONhAAEEBZE?H?CBC?DA    NM:i:1  X1:i:1  MD:Z:33A2
SOLEXA-1GA-2_2_FC20EMB:5:102:214:278    0       chr1    10148   25      36M     *       0       0       CCCTAACCCTAACCCTAACCCTAACCCTAACCTAAC    hbfhhhXUYhT_ULZdLRTKNIMIKGLJCHFFJQJN    NM:i:0  X0:i:1  MD:Z:36
SOLEXA-1GA-2_2_FC20EMB:5:195:284:685    16      chr1    10149   25      36M     *       0       0       CCAAACACTAACCCTAACCCTAACCCTAACCTAACC    ><>B@>?>?D>>?B?D>DBC?E@BDHAKCEKERLOO    NM:i:1  X1:i:1  MD:Z:29A3A2
SOLEXA-1GA-2_2_FC20EMB:5:35:583:827     0       chr1    10150   25      36M     *       0       0       CTAACCCTAAACCTAACCCTAACCCTAACCTAACCA    hhW_X]MXNOHQQWMILHGIFMJGJLCFGGJAKIEH    NM:i:1  X1:i:1  MD:Z:10A24A0
SOLEXA-1GA-2_2_FC20EMB:5:248:130:724    0       chr1    10152   25      36M     *       0       0       AACCCTAACCCTAACCCCAACCCTAACCTAACCCTA    hchPhc___cWS[VRObRXDTOUSJLXOA@LGFMC@    NM:i:1  X1:i:1  MD:Z:17C18
SOLEXA-1GA-2_2_FC20EMB:5:236:644:107    16      chr1    10154   25      36M     *       0       0       CCCTAACCCTAACCCTAACCCTAACCTAACCCTAAC    ICF=A@BHGEGGDFLIKKYYGD^CahMaShfhNhhh    NM:i:0  X0:i:1  MD:Z:36
SOLEXA-1GA-2_2_FC20EMB:5:165:628:70     16      chr1    10155   25      36M     *       0       0       CCTAACCCTAACCCTAACCTTTACATAACCCTAACC    >D?BAEAA?E=JGGBK>FKDGFJPVHSFTT\QQMch    NM:i:1  X1:i:1  MD:Z:11A2T1T19
SOLEXA-1GA-2_2_FC20EMB:5:108:485:455    16      chr1    10156   25      36M     *       0       0       CTAACCCTAACCCTAACCCTAACATAACCCTAACCC    CFVIUIIIbGPROGhRhhhIhhhhhhhhhhhhfhhh    NM:i:1  X1:i:1  MD:Z:12A23
SOLEXA-1GA-2_2_FC20EMB:5:240:501:237    0       chr1    10158   25      36M     *       0       0       AACCCTAACCCTAACCCTAACCTAACCCTAACCATA    hhhchg_ORNbX]RMZLREQMNTNFLDPMDDDEDKL    NM:i:1  X1:i:1  MD:Z:33A2
SOLEXA-1GA-2_2_FC20EMB:5:258:882:389    16      chr1    10215   25      36M     *       0       0       CTAACCCTAAACCTAACCCCTAACCCTAACCCTAAA    IMZO>?EGDRRhhdXKVNRZKhhhhhhhNhhhhhhh    NM:i:1  X1:i:1  MD:Z:25A10

 

SAM -> BAM 변환

samtools view -Sb {파일명}.sam > output.bam

$ samtools view -Sb output.sam > output.bam

[E::sam_parse1] no SQ lines present in the header

이 에러가 떠서 찾아보니 sam 파일을 생성할 때 -h 옵션으로 헤더를 생성해야한다고 한다.

$ samtools view -h wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam > output-header.sam $ samtools view -Sb output-header.sam > output.bam

 

BAM 파일이 만들어졌으면 정렬하고 인덱스를 저장할 .bai를 만들어야한다.

$ samtools sort output.bam
$ samtools index output.bam

 

BAM 파일 읽기

검색해보면 많은 사용법들이 나오니 찾아보길 바란다.

간단히 Total read 하는 것을 가져왔다.

$ samtools idxstats wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam | cut -f 1,3 | more

[E::idx_find_and_load] Could not retrieve index file for 'wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam'
samtools idxstats: fail to load index for "wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam", reverting to slow method
chr1    179687
chr2    133152
chr3    111106
chr4    77087
chr5    94008
chr6    96140
chr7    96286
chr8    66273
chr9    77261
chr10   72889
chr11   88933
chr12   90121
chr13   43075
chr14   56976
chr15   56850
chr16   71169
chr17   90704
chr18   32388
chr19   90121
chr20   46817
chr21   18025
chr22   42468
chrM    3555
chrX    49776
*       0

 

마치며

이렇게 유전자 데이터를 확인해보았다.

이제 IGV로 유전자 데이터를 그래픽하게 표현하고 분석하면 된다.

 

처음보는 데이터들이고, 생명공학을 이해하긴 어렵겠지만 오픈소스를 활용해 접할 수 있어 좋은 기회가 되었다.