Table of Contents
Whole-genome alignment
유전체 서열의 쌍을 1:1로 비교하여 dot plot을 만드는 GUI 프로그램으로는 gepard가 있다. 단일 염기서열을 자체 비교하면 반복 구조라든가 서열의 양 끝에 나타나는 겹침을 검출할 수 있다. 입력 파일이 multi FASTA인 경우 파일 내의 순서 그대로 서열을 이어 붙여서 dot plot을 생성한다. Windows 환경에서 쓰려면 여기의 것을 다운로드하여 사용하는 것이 낫다.
GUI 방식으로 실행하려면 다음과 같이 입력한다.
$ java -jar -Xms1024m -Xmx1024m /usr/local/apps/gepard/dist/Gepard-1.40.jar
명령행 방식으로 실행하려면 최소한 서열 2종과 매트릭스, 그리고 출력 파일(default: PNG 포맷)을 지정해야 한다.
$ java -cp /usr/local/apps/gepard/dist/Gepard-1.40.jar org.gepard.client.cmdline.CommandLine -seq1 genome_1.fasta -seq2 genome_2.fasta -matrix /usr/local/apps/gepard/resources/matrices/edna.mat -outfile out.png
MUMmer and dnadiff
MUMmer(v.3)는 유전체 서열 정렬 프로그램의 원조이자 대명사라고 해도 과언이 아니다. 버전 4부터는 다룰 수 있는 유전체 크기에 제한이 없고 다중 쓰레드를 지원한다. 두 개의 입력 서열 파일로부터 MUM(Maximal Unique Matches)를 추출하여 dot plot을 보여주거나(mummer), MUM의 클러스터로부터 alignment를 생성하며(nucmer), 정렬된 영역의 위치와 길이, % identity 등의 정보를 추출하는 show-coords 등의 다양한 부속 유틸리티로 구성되어 있다. MUMmer는 자체적인 GUI는 갖고 있지 않지만 출력물인 delta 파일을 mummerplot으로 전환하면 gnuplot의 입력이 되어 화면으로 표시를 하거나 이미지 파일로 출력을 하도록 배려해 놓았다. MUMmer는 실행 속도가 빠르기 때문에 SNP 및 ANI 계산 등에서 BLAST를 대신하여 쓰이기도 한다. 매우 방대한 기능의 프로그램이므로 매뉴얼과 활용사례를 보면서 차근차근히 익히도록 한다.
MUMmer의 부속 wrapper 스크립트인 dnadiff는 script로서 nucmer와 관련 유틸리티를 실행하여 두 서열의 차이에 관련한 종합적인 리포트를 작성해 준다. 예를 들어서 같은 종에 속하는 두 균주의 유전자를 비교하거나, 동일 생명체로부터 유래한 sequence read를 서로 다른 de novo assembler를 이용하여 만든 contig 서열에 어떤 차이가 있는지를 알아보고자 할 때 유용하다. -p <prefix>로 지정하지 않으면 'out'을 접두사로 하는 일련의 결과 파일이 만들어진다. Alignment, difference 및 SNP에 관한 종합 리포트는 out.report 파일이다. 두 유전체 서열의 identity(%) 수치는 [Alignments]의 1-to-1 섹션에 나오는 AvgIdentity 라인을 참조하라.
$ dnadiff reference.fa query.fa $ dnadiff -d (Delta file)
dnadiff가 만들어내는 결과 파일은 다음과 같다. '-p|--prefix'로 지정하지 않으면 'out'을 기본 prefix로 사용한다. show-diff는 MUMmer 3.0 공식 페이지에서 설명을 제공하지는 않는다('man show-diff'를 입력하면 매뉴얼은 나온다).
- out.report - Summary of alignments, differences and SNPs
- out.delta - Standard nucmer alignment output
- out.1delta - 1-to-1 alignment from delta-filter -1
- out.mdelta - M-to-M alignment from delta-filter -m
- out.1coords - 1-to-1 coordinates from show-coords -THrcl .1delta
- out.mcoords - M-to-M coordinates from show-coords -THrcl .mdelta
- out.snps - SNPs from show-snps -rlTHC .1delta
- out.rdiff - Classified ref breakpoints from show-diff -rH .mdelta
- out.qdiff - Classified qry breakpoints from show-diff -qH .mdelta
- out.unref - Unaligned reference sequence IDs and lengths
- out.unqry - Unaligned query sequence IDs and lengths
리포트 파일(out.report)의 내용은 다음과 같다. Refenece(GCF_000063585.1)와 query(GCF_000017045.1) 전부 complete genome sequence이며, reference는 16 kb의 플라스미드가 있다. 리포트 파일 내용의 설명은 dnadiff REAME 파일을 참조하라.
/home/hyjeong/temp_20200624/ATCC_3502.fna /home/hyjeong/temp_20200624/Hall.fna
NUCMER
[REF] [QRY]
[Sequences]
TotalSeqs 2 1
AlignedSeqs 1(50.00%) 1(100.00%)
UnalignedSeqs 1(50.00%) 0(0.00%)
[Bases]
TotalBases 3903260 3760560
AlignedBases 3724663(95.42%) 3719169(98.90%)
UnalignedBases 178597(4.58%) 41391(1.10%)
TotalSeqs 1 1
AlignedSeqs 1(100.00%) 1(100.00%)
UnalignedSeqs 0(0.00%) 0(0.00%)
<이하 생략>
SNP 파일을 VCF로 전환하려면 all2vcf 유틸리티의 mummer 명령을 이용하라.
두 염기서열 FASTA file의 identity만을 계산하고 싶다면 identity.sh 스트립트를 사용하라. 내가 만든 스크립트였던가, 혹은 어디서 퍼 온 것인가? 잘 기억이 나지 않는다.
$ cat identity.sh
#!/bin/sh
nucmer $1 $2
delta-filter -1 out.delta > filtered.delta
show-coords -THrcl filtered.delta > filtered.coords
awk '{product_sum += ($5 * $7); len += $5} END {print product_sum / len}' filtered.coords > out.identity
echo -n 'Percent identity: '
cat out.identity
$ identity.sh A.fa B.fa
Percent identity: 99.9693
NucDiff 및 기타
Nucmer를 사용한 두 서열 set의 비교 결과로부터 어떤 차이점이 있는지를 체계적으로 알아보고자 한다면 NucDiff를 추천한다. NucDIff는 SNP는물론 structural rearrangement를 유형별로 정리하여 GFF와 VCF 파일로 결과를 제공한다. 다음은 Clostridium botulinum ATCC 3502를 reference로 간주하여 Hall A 균주의 유전체 서열을 비교하는 NucDiff 실행 방법이다. 마지막 두 인수는 출력 디렉토리와 결과 파일의 앞에 붙는 접두사이다. 총 9개의 결과 파일은 out_dir/results 디렉토리에 저장된다. 전체 요약은 test_stat.out 파일에 수록된다.
$ nucdiff ATCC_3502.fna Hall.fna out_dir test
NucDiff가 정의하는 ‘차이’는 Global, Local 및 Structural의 세 가지 범주로 구분되며, 상세한 정의는 위키 페이지에서 상세히 설명하였다.
Mauve는 여러 유전체 서열을 정렬하여 시각적인 결과를 보여주는 GUI 프로그램이다.
dnadiff 결과로부터 균주 특이적 영역 찾기
두 미생물 유전체(A & B)를 서로 비교하여 특이적인 영역을 찾는 방법을 생각해 보자. A 균주의 유전체는 RefSeq 등 데이터베이스에서 다운로드 가능한 상태이고, B 균주는 아직 공개된 정보가 없어서 실험실에서 일루미나 기법으로 short read를 생성했다고 가정하자. A 유전체를 reference로 삼아서 B의 read를 매핑하면 zero coverage region으로부터 A 특이적인 유전체 염기서열을 찾을 수 있고, unmapped read를 회수하여 조립하면 B 특이적인 염기서열을 찾을 수 있다. 매핑시 percent identity threshold는 95% 정도로 높게 잡는 것이 무난할 것이다. 선정된 균주 특이적 영역은 상대방 유전체에 대하여 BLASTN 검색을 하여 매치하지 않음을 최종적으로 점검하는 것이 바람직하다.
비교 대상 균주가 전부 contig 상태의 유전체 정보로만 존재한다면, 이로부터 fake read를 생성하여 마찬가지의 방법으로 분석을 실시하면 될 것이다. 그러나 MUMmer(dnadiff)를 이용하여 유전체 단위의 직접 비교를 실시한 뒤 unaligned region에 대한 정보(위치 및 염기서열)를 추출하는 방법을 알아보도록 하자.
Dnadiff를 이용한 두 서열 세트간의 비교에서 정렬 위치에 대한 정보를 확인하려면 out.mcoords 파일을 참조하면 된다. out.1coords 파일은 왜 이러한 목적에 적합하지 않은지 각자 생각해 보라. out.mcoords 파일은 MUMmer-NUCmer의 정렬 결과(out.mdelta)를 패키지에 포함된 유틸리티인 show-coords가 정리하여 reference coordinate 위치에 맞추어 출력한 것이다. Show-coords 실행에 사용된 옵션은 ‘-THrcl’에 해당한다. Out.mcoords 파일은 총 13개의 컬럼으로 이루어지며, 각각이 의미하는 바는 다음과 같다.
- [S1] start of the alignment region in the reference sequence
- [E1] end of the alignment region in the reference sequence
- [S2] start of the alignment region in the query sequence
- [E2] end of the alignment region in the query sequence
- [LEN 1] length of the alignment region in the reference sequence
- [LEN 2] length of the alignment region in the query sequence
- [% IDY] percent identity of the alignment
- [LEN R] length of the reference sequence
- [LEN Q] length of the query sequence
- [COV R] percent alignment coverage in the reference sequence
- [COV Q] percent alignment coverage in the query sequence [FRM] reading frame for the reference and query sequence alignments respectively
- [TAGS] the reference and query FastA IDs respectively. All output coordinates and lengths are relative to the forward strand of the reference DNA sequence.
이 파일에서 query의 위치 정보에 해당하는 3번째와 4번째 컬럼의 값을 각 라인(alignment)에 대하여 뽑아낸 뒤 query data set을 구성하는 각 서열에 대하여 병합한 다음, 그것에 대한 complementary region을 추출하여 염기서열로 추출하면 된다. 3번째와 4번째 컬럼은 서열 상에서 구간의 시작과 끝을 의미하는 데이터이므로 BED 포맷으로 전환한 뒤 bedtools를 이용하여 병합 등의 연산을 간편하게 수행할 수 있다. 먼저 out.mcoords로부터 query sequence set의 alignment 위치 정보를 추출하여 BED 파일을 만들어 보자. GFF 포맷과 마찬가지로 위치 정보는 항상 작은 수가 앞으로 나오도록 해야 하므로, out.mcoords의 3번과 4번 컬럼을 비교하여 작은 수가 BED 파일의 두번째 컬럼이 되도록 조정해야 한다. Query sequence를 reverse complementary로 만들어야 정렬하는 경우가 여기에 해당한다. 그러고 나서 bedtools merge 기본 명령으로 병합을 하자. 병합을 하려면 반드시 각 서열(염색체 혹은 contig)과 시작 위치에 따라서 미리 정렬해 두어야 한다.
$ awk -vOFS="\t" '{if($3<$4) print $13,$3,$4;else print $13,$4,$3}' out.mcoords > query.bed
$ sort -k1,1 -k2,2n query.bed > query.bed.sorted
$ bedtools merge -i query.bed.sorted > query.bed.sorted.merged
$ cat query.bed.sorted.merged
NC_009698.1 1 68108
NC_009698.1 70053 248022
NC_009698.1 248024 349397
...
NC_009698.1 2523054 2679762
NC_009698.1 2682911 3088750
NC_009698.1 3088857 3886916
병합한 구간을 염색체 전체에 대해서 빼면(bedtools complement 이용) query 특이적인 염기서열 의 구간 정보를 얻는다. 이 작업에는 query 서열의 길이를 정의한 파일(query.genome)이 필요하다. EMBOSS에서 제공하는 infoseq 명령어의 출력을 awk로 다듬으면 query.genome 파일을 쉽게 만들 수 있다. 아래의 사례에서는 bedtools complement를 실행한 다음 특이적 영역의 길이($3 - $2 + 1)가 너무 짧은 것은 필터링하였다. 최소 길이의 기준은 100 bp로 하였다.
$ infoseq Hall_A.fna | awk -vOFS="\t" '!/^USA/ {print $3, $6}' > query.genome
$ cat query.genome
NC_009698.1 3760560
$ bedtools complement –i query.bed.sorted.merged -g query.genome > query-specific.bed
$ awk –vOFS="\t" '$3-$2 > 99' query-specific.bed > query-specific.bed.filt
$ cat query-specific.bed.filt
NC_009698.1 68108 70053
..
NC_009698.1 3088750 3088857
특이적 염기서열의 위치 정보를 이용하여 EMBOSS seqret 명령으로 서열을 추출하면 된다.
$ cat query-specific.bed.filt NC_009698.1 68108 70053 .. NC_009698.1 3088750 3088857 $ cat query-specific.bed.filt | while read a b c > do > seqret -sbegin $b -send $c -auto -stdout Hall_A.fna:$a >> query-specific.fa > done
FASTA file과 BED file을 이용하여 서열을 추출하려면 seqtk subseq 명령어를 써도 된다. Seqtk는 추출된 서열의 ID에 위치 정보를 포함시키는 점이 편리하다.
$ seqtk subseq Hall_A.fna query-specific.bed.filt > query-specific.fa
이미 BED 파일에 대하여 짧은 구간은 제거하였지만, 염기서열을 추출하여 파일을 만든 뒤에 다른길이 조건에 대하여 필터링을 할 수도 있다. 대단히 빠른 FASTA/Q file 조작용 툴킷을 자처하는 SeqKit를 사용하면 길이를 기준으로 필터를 적용할 수 있다. 'seqkit seq -
--min-len) 150'은 최소 길이 150 bp인 서열만을 추출하는 것이다. -M 옵션을 이용하면 특정 길이 이하의 서열을 추출하는 반대의 동작을 한다. 유사한 명령으로 seqtk seq -L INT infile.fa(or infile.fq) 명령을 실행하면 INT보다 짧은 서열을 입력 FASTA/Q 파일에서 제거한다.
$ seqkit seq -m 150 query-specific.fa > filtered.fa $ seqtk seq -L 150 query-specific.fa > filtered.fa
유전체 염기서열의 특정 구간 정보를 BED 파일로 확보하였다면 여기에 어떤 기능을 갖는 유전자가 위치하고 있는지를 알고 싶을 것이다. GenBank 파일에서 유전자 정보를 수록한 BED 파일을 만든 뒤, 이를 위에서 만든 query-specific.bed(.filt)와 함께 bedtools intersect 명령어에 인수로 제공하면 된다. 첫 번째 BED 파일에는 유전자의 locus tag과 기능 정보를 4번째 컬럼부터 삽입해 두면 최종 결과 파일에도 이 정보가 상속되므로 후속 작업에 매우 편리하다. gbkInfo.pl을 활용하면 첫 번째 BED 파일을 만들 수 있다. 두 BED 파일의 경계에 걸치는 유전자를 처리하기 위하여 bedtools intersect 명령 실행 시 '-f 0.3' 옵션을 주었다. 이는 유전자 정보를 담고 있는 첫 번째 BED 파일의 정보에 대하여 최소 30%가 겹쳐야 함을 의미한다. 이는 단지 하나의 사례이므로 적절히 바꾸어도 된다. 예를 들어 '-f 1.0'으로 설정하면 specific region 안에 전 영역이 완전히 포함되는 유전자들만 추출될 것이다.
$ gbkInfo.pl Hall_A.gbk # Hall_A.gbk.txt 파일 생성(query genome)
# seqID, start, end, locus tag, product 컬럼을 추출한다.
# 필요하다면 다른 컬럼을 추가해도 된다.
# product 컬럼이 빈 것은 pseudo를 의미한다.
$ awk -F"\t" -vOFS="\t" '!/^#/{print $2,$7,$8,$4,$11}' Hall_A.gbk.txt > Hall_A.gene.bed
# Hall_A.gene.bed의 seqID에는 소수점 이하의 버전 번호가 붙어있지 않다.
# 따라서 query-specific.bed(.filt)의 첫 번째 컬럼에서 이를 제거해야 한다.
$ cp query-specific.bed query-specific.bed.mod
$ sed –I ‘s/\.1//g’ query-specific.bed.mod
$ bedtools intersect –f 0.3 –a Hall_A.gene.bed –b query-specific.bed.mod
CP000727 441348 442319 CLC_0444 conserved hypothetical protein
CP000727 497058 498272 CLC_0493 putative ABC transporter, permease protein
CP000727 498262 498978 CLC_0494 ABC transporter, ATP-binding protein
…(이하 생략)
minimap2의 활용
이 항목은 참조 서열에 대한 매핑(reference mapping) 및 시각화에서 다루는 것이 원칙적으로 맞을 것이다. 그러나 reference에 대하여 다수의 짧은 read를 붙이는 것이 아니라 이미 조립이 완료된 긴 염기서열을 1:1로 비교하는 것을 주된 내용으로 하므로, whole-genome alignment 항목에서 설명하는 것이 더 낫다고 생각한다.
minimap(논문)은 원래 long read와 reference의 빠르고도 대략적인 alignment를 위해서 만들어진 프로그램이다(심지어 20%의 차이까지 수용하는!). 이는 BWA의 개발자인 Heng Li가 개발하였다. 현재 쓰이는 minimap2(2018년도 논문)는 genome과 spliced nucleotide의 정렬 목적으로 쓰일 수 있게 더욱 발전하였다. Heng Li의 블로그에서 minimap2로의 전환 및 BWA의 미래에 관하여 쓴 글(링크)을 참조하라. 또한 공식 문서의 Getting Started를 통해서 minimap2가 여러 상황에서 어떻게 쓰이는지 감상해 보기 바란다.
2018년도 논문에서는 minimap2의 특성을 다음과 같이 소개하였다.
Minimap2 is a versatile mapper and pairwise aligner for nucleotide sequences. It works with short reads, assembly contigs and long noisy genomic and RNA-seq reads, and can be used as a read mapper, long-read overlapper or a full-genome aligner. Minimap2 is also accurate and efficient, often outperforming other domain-specific alignment tools in terms of both speed and accuracy.
minimap2의 여러 사용례 중에서 미생물 유전체를 1:1로 비교하여 SNP를 추출하는 방법을 정리해 보았다. 이를 더욱 확장하면 정렬이 이루어지지 않는 영역 정보를 추출하는 것도 가능할 것이다.
minimap2 misc 문서의 Calling variants from haploid assemblies 항목을 주로 참조하였다. minimap2의 여러 옵션에 대해서는 따로 공부를 하는 것이 좋을 것 같다. PAF(a Pairwise mApping Format) 파일의 설명은 여기를 참조하라. paftools.js는 minimap2과 같이 설치되었으므로 어디에 있는지 잘 찾아보기 바란다. 아래의 실행 사례는 위에서 dnadiff 분석의 샘플로 사용했던 보툴리늄 균주 2종의 유전체를 사용한 것이다.
# keeping this file is recommended; --cs required! $ minimap2 -cx asm5 -t8 --cs ref.fa asm.fa > asm.paf # sort by reference start coordinate $ sort -k6,6 -k8,8n asm.paf > asm.srt.paf $ k8 paftools.js call asm.srt.paf > asm.var.txt 3686527 reference bases covered by exactly one contig 794 substitutions; ts/tv = 3.671 40 1bp deletions 25 1bp insertions 0 2bp deletions 3 2bp insertions 2 [3,50) deletions 5 [3,50) insertions 3 [50,1000) deletions 9 [50,1000) insertions 6 >=1000 deletions 4 >=1000 insertions $ head -n 5 asm.var.txt R NC_009495.1 0 56650 V NC_009495.1 8776 8776 1 60 - tttaaatgaagaaaa NC_009698.1 8776 8791 + V NC_009495.1 52826 52827 1 60 t c NC_009698.1 52841 52842 + V NC_009495.1 55514 55515 1 60 a g NC_009698.1 55529 55530 + V NC_009495.1 55533 55534 1 60 t c NC_009698.1 55548 55549 +
asm.var.txt에서 R로 시작하는 줄은 하나의 query contig가 덮는 영역을, V로 시작하는 줄은 variant를 뜻한다. paftools.js가 variant calling을 할 때 취하는 alignment의 최소 길이가 정해져 있으므로 short read를 활용할 때에는 조정이 필요하다.