참조 서열에 대한 매핑(reference mapping) 및 시각화

BWA나 bowtie2와 같은 대중적인 short read aligner(mapper)를 써도 좋지만 prokaryotic genome 유래 데이터라면 SSAHA2 또는 smalt도 좋다. Read alignment 결과는 SAM 또는 BAM(specification documentation)으로 출력되고, 이는 SAMtools를 사용하여 서로 전환 가능하다. SAMtools는 alignment file의 편집과 매핑 관련 수치 추출 등 다양한 기능을 수행한다. 이 섹션에서는 NCBI SRA(Short Read Archive)에서 E. coli BL21의 sequence read file 일부를 직접 다운로드하여 이용하는 방법을 설명한다. Reference sequence는 E. coli B str. REL606의 유전체(GCF_000017985.1, ftp)를 사용하였다. 아래에 보인 사례에서는 FASTA file을 다운로드하였지만, GenBank 파일을 받아서 이로부터 FASTA 및 GFF3를 추출하여 사용하는 것도 좋다. FASTA 추출은 seqret(EMBOSS package) 혹은 bp_seqconvert.pl을, GFF 추출은 bp_genbank2gff3.pl을 사용한다. ‘bp_’로 시작하는 Perl 스크립트는 BioPerl(conda base environment로 설치)에서 제공하는 것이다. SRA 데이터를 다운받으려면 SRA-Toolkit 패키지의 fastq-dump가 필요하다.

샘플로 사용할 Illumina sequencing read는 E. coli BL21(TaKaRa)에서 유래한 SRX1021661의 단일 run 결과(SRR2014531)로부터 2백만 read(2 x 101 nt cycle)만을 추출한 것이며, read mapper로는 bowtie2를 사용한다. Bowtie2에는 paired read의 간격이나 방향 설정 등 입출력물의 형식 및 정렬 조건에 대한 다양한 옵션을 줄 수 있는데, 여기에서는 일루미나 paired end read에 해당하는 기본 조건(forward-reverse, no distance constraint)으로 하였다. Alignment의 기본 조건은 –end-to-end인데, 트리밍을 하지 않은 read라면 –local이 더 나을 수도 있다. 파라미터의 조합에 따라서 speed/sensitivity/accuracy이 달라지므로 자세한 사항은 bowtie2 매뉴얼 혹은 'Getting started' 섹션을 참조하라.

# 실습용 raw data의 설명은 https://www.ncbi.nlm.nih.gov/sra/SRR2014531 참조
$ fastq-dump -X 10000000 --split-files SRR2014531
$ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/017/985/GCF_000017985.1_ASM1798v1/GCF_000017985.1_ASM1798v1_genomic.fna.gz
# 실습용 reference sequence의 accession은 NC_012967.1 또는 CP000819.1이다.
$ gzip -d GCF_000017985.1_ASM1798v1_genomic.fna.gz
$ bowtie2-build GCF_000017985.1_ASM1798v1_genomic.fna REL606
# 다음의 명령은 --fr -I 0 -X 500을 기본 조건으로 이용한다.
# -S(표준 출력) 또는 –S <SAM file name>을 지정하지 않으면 BAM file로 출력된다.
$ bowtie2 -p 8 -x REL606 -1 SRR2014531_1.fastq -2 SRR2014531_2.fastq -S BL21.sam
# unpaired read가 있다면 -U fastq_file로 공급한다.
# SAM <-> BAM 파일 전환을 위해 samtools view 명령을 사용한다.
# samtools에서 -S는 입력 파일의 포맷이 SAM임을 지정하는 것이다.
# samtools 1.9에서는 입력 파일 포맷을 자동으로 검출하므로 –S 옵션은 무시된다.
# 아래에서는 하위 호환성(버전 0.1.19까지)을 위하여 –S를 삽입하였다.
# -b는 출력 파일 포맷이 BAM임을 지정한다.
# -o OUTFILE은 samtools 버전 1.9 이상부터 사용 가능하다.
$ samtools view -b -S -o BL21.bam BL21.sam

SAM 파일의 대부분을 구성하는 read alignment 필드에서 두 번째 필드(flag, 0-255 범위의 숫자)는 read의 mapping 상태를 비트 단위로 나타낸다. 특정 flag를 만족하는 alignment를 추출하려면 samtools view -f INT를, 만족하지 않는 것을 추출하려면 samtools -view -F INT 명령을 사용하면 된다. Reference에 반복 서열이 존재하면 하나의 read가 여러 곳에 붙어서 각각 별도의 alignment를 생성하게 되므로(primary vs. secondary alignment), mapped read의 수치와 관련된 작업을 하려면 주의가 필요하다. 뿐만 아니라 chimeric alignment인 경우에도 하나의 read가 representative 및 supplementary alignment로 여러 차례 나타날 수 있다. SAM flags를 해독하려면 ‘samtools flags 4’와 같이 입력하거나 Decoding SAM flags 웹사이트를 이용하라. 예를 들어 0x04는 segment unmmaped를 의미한다. ‘0x’는 16진수임을 명시하기 위한 것이다(예: 10 = 0xa). getinsersize.py는 매핑 데이터를 분석하여 라이브러리 크기에 대한 정보를 제공한다. 이 스크립트에 대한 상세한 설명은 Estimating paired-end read insert length from SAM/BAM files를 참조하라.

$ samtools flagstat BL21.bam
$ samtools view –f 0x04 BL21.bam | wc -l     # unmapped read의 수
$ samtools view –F 0x04 BL21.bam | wc –l     # mapped read의 수
# -c 옵션을 주면 해당되는 레코드, 즉 read의 수를 반환한다.
$ samtools view –f 0x04 –c BL21.bam           # unmapped read의 수
$ samtools view –F 0x40 BL21.bam | cut –f1 | sort | uniq | wc –l   # uniquely mapped read의 수
$ samtools view BL21.bam | getinsertsize.py –   #library insert size의 측정

getinsertsize.py 스크립트가 출력하는 read span 정보가 정확히 무엇을 의미하는가? SAM definition에 의하면 read span은 매핑된 왼쪽 read의 시작에서 오른쪽 read의 끝까지를 의미한다. 다시 말해서 {paired read의 간격 + (2 x read length)}에 해당한다(출처).

Smalt는 mapping을 하기 전에 reference sequence의 hash table을 만드는 인덱싱을 먼저 해야 한다. Hashed word의 크기를 -k 옵션으로 미리 지정하지 않으면 13을 택하게 된다. Bacterial genome의 인덱싱에는 -k 11 -s 1이 적당하다. -s는 sampling step size를 의미한다. 매핑 단계에서 쓰레드를 지정하려면 –n NUM 옵션을 사용하라.

$ smalt index –k 11 –s 1 referencek11s1 reference.fa
# 인덱싱 결과 referencek11s1.sma 및 referencek11s1.smi 파일이 생긴다
$ smalt map –f bam –o mapping.bam referencek11s1 read_1.fq read_2.fq

다음으로 평균적인 read depth(또는 coverage)를 산출하는 방법을 소개한다. Depth란 de novo assembly 또는 read mapping으로 재구성된 염기 서열의 주어진 위치에 몇 개의 read가 정렬하고 있는지를 수치로 나타낸 것이다1. De novo assembly에서 sequencing에 의한 유전체의 평균 coverage를 구하는 가장 간단한 방법은 read 수(N)와 평균 read length(L)를 곱한 것을 유전체 크기(G)로 나눈 것, 즉 N x L/G를 산출하는 것이다. Read mapping을 통하여 SAM/BAM 파일이 만들어진 상태라면 samtools를 사용하여 보다 정확한 계산을 할 수 있다. 이를 위해서는 read alignment를 reference 서열 내에서의 mapping position으로 정렬하는 일(samtools sort)이 선행되어야 한다. samtools sort –n은 read name을 기준으로 정렬한다.

$ samtools sort -o BL21_sorted.bam BL21.bam
# mapping부터 sort까지 한번에 실시하기
$ bowtie2 -p 8 -x REL606 -1 SRR2014531_1.fastq -2 SRR2014531_2.fastq | samtools sort > BL21_sorted.bam
$ samtools depth -a BL21_sorted.bam | awk '{c++;s+=$3}END{print s/c}' # awk문 안에서 c 변수를 계산하지 말고 '{s+=$3}END{print s/NR}'이라고 해도 된다.

samtools depth의 -a 옵션은 zero coverage를 포함한 reference의 모든 영역을 추출하기 위한 것이다. 위에서 소개한 명령어는 reference sequence가 몇 개의 contig로 나뉘어 있는지 관계 없이 평균 read depth를 출력한다. BBmap 패키지의 pileup.sh 명령어는 reference를 이루는 각 서열에 대한 길이, 평균 fold coverage, covered percent, covered base, mapping된 read의 방향 등 상세한 정보를 tab-delimited format으로 출력한다.

$ pileup.sh in=mapping_sorted.bam out=stats.txt

각 contig에 대하여 몇 개의 read가 mapping되었는지를 간략하게 살펴보고 싶다면 다음 명령어를 사용한다.

$ samtools idxstats mapping_sorted.bam

뒤에서 소개할 bedtools의 genomecov를 사용하면 coverage와 관련한 보다 복잡한 계산을 할 수 있다.

Read alignment의 시각화에는 IGV(Integrated Genomics Viewer) 또는 tablet을 사용하라.

프로그램에서 열게 될 BAM file은 sort 및 index가 되어 있어야 한다. 또한 reference genome sequence(FASTA file)와 유전자 정보(annotation) 파일도 필요하다. 유전자 정보는 보통 GFF3 파일로 준비해야 한다. IGV에서는 Genomes→Create .genome File… 메뉴에서 FASTA file과 gene file(GFF3)을 입력하여 사용자의 홈 디렉토리에 각 reference에 대한 .genome 파일을 만들게 되며, 이렇게 만들어진 파일은 다른 사용자에게도 공유할 수 있다.

다음 명령어는 다소 복잡해 보이지만 파일의 다운로드, 압축해제 및 전환 과정을 파이프로 연결함으로써 중간 단계에서 파일을 저장하는 과정을 생략한 것이다. Prokaryotic genome이라면 GFF3 파일에서 gene과 CDS feature만 남기는 것이 편리하다. NCBI에서 다운로드한 GenBank 파일을 bp_genbank2gff3.pl로 처리하면 gene/mRNA/CDS/exon feature가 전부 생기는 반면 같은 사이트에서 다운로드한 GFF3 파일에는 CDS feature만 존재한다.

$ curl --output - ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/017/985/GCF_000017985.1_ASM1798v1/GCF_000017985.1_ASM1798v1_genomic.gbff.gz | zcat - | bp_genbank2gff3.pl -in stdin -out stdout > temp.gff3
$ awk '$3!="mRNA" && $3!="exon"' temp.gff3 > REL606.gff3; rm temp.gff3
(또는)
$ curl --output - ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/017/985/GCF_000017985.1_ASM1798v1/GCF_000017985.1_ASM1798v1_genomic.gff.gz | zcat - > REL606.gff3
$ samtools index BL21_sorted.bam

Mapping에 사용한 FASTA file(NCBI에서 직접 다운로드)의 서열 ID는 NC_12967.1의 형태이나 GFF3 파일의 첫번째 컬럼에 보이는 서열 ID에서는 버전 번호(.1)가 없으므로 GFF3 파일을 수정하여 이를 서로 일치시키도록 한다. 그러나 GenBank 파일로부터 seqret(EMBOSS)로 만들어낸 FASTA 파일의 서열 ID에는 버전 번호가 없다.

$ sed -r 's/^NC_012967/NC_012967.1/' REL606.gff3 > REL606.mod.gff3

만약 FASTA file에 여러 개의 서열이 포함되어 있다면 위 방법으로는 버전 번호를 한번에 제거하기 어렵다. 따라서 sed에서 좀 더 정교한 치환을 실시해야 한다. 아래의 사례에서는 서열 ID에서 가장 마지막에 오는 점 이후의 모든 문자를 제거한다.

$ sed –r 's/^(/+)\..+$/\1/' file.gff3 > file.mod.gff3

이렇게 만들어진 FASTA file과 GFF3 및 BAM file을 IGV에서 시각화해 본다. tablet에서도 비슷한 방법을 적용할 수 있을 것이다.

Samtools depth는 reference의 각 염기 위치 혹은 구간에 대한 read depth 정보를 라인 단위로 출력한다. 따라서 read depth가 일정 기준을 만족하지 않는 위치를 구간으로 묶어서 출력하려면 bedtools를 사용하는 것이 편리하다. Bedtools는 유전체 분석에서 흔히 다루게 되는 파일인 BAM. BED, GFF/GTF, VCF를 대상으로 하여 구간(genomic interval)에 대한 다양한 조작을 하는 도구이다. 여기에는 intersect, merge, count, complement, shuffle 등이 포함된다. Coverage 관련 계산에 사용되는 bedtools의 명령어는 genomecov이다.

다음의 예제는 BAM file에서 read coverage가 9 미만인 곳을 출력한 다음 병합하여 구간으로 출력하는 방법을 보여주고 있다.

$ bedtools genomecov –bga –ibam infile.bam |
> awk ‘$4 < 9’ | bedtools merge –i - > RESULT.bed

Bedtools를 사용하면 참조 서열이 read mapping에 의하여 어느 정도나 cover되었는지를 계산할 수 있다. 다음 명령행은 BAM 파일에서 read depth가 0인 영역(zero)과 그렇지 않은 영역(nonzero)의 길이를 합산한 다음 이로부터 % of reference genome covered를 산출하는 사례이다(Biostars tutorial). BAM 파일이 입력되는 경우 -g reference.fasta 옵션은 필요하지 않다.

$ zero=$(bedtools genomecov -ibam infile.bam -g reference.fasta -bga | awk '$4==0 {bpCountZero+=($3-$2)} {print bpCountZero}' | tail -1)
$ nonzero=$(bedtools genomecov -ibam infile.bam -g reference.fasta -bga | awk '$4>0 {bpCountNonZero+=($3-$2)} {print bpCountNonZero}' | tail -1)
# Round up to 6 decimal places and print
$ percent=$(bc <<< "scale=6; ($nonzero / ($zero + $nonzero))*100")
$ echo $percent

Quality가 좋지 않거나 sample에만 특이적인 부분에서 유래한 read 및 오염에서 유래한 read는 unmapped read로 남게 된다. 이것만을 별도로 추려내면 sample의 상태를 파악할 수 있는 자료로 쓸 수 있다. Samtools를 이용하여 BAM file로부터 unmapped read를 fastq 파일로 추출하는 방법을 알아본다.

$ samtools view –bf 0x04 BL21.bam > unmapped.bam
$ samtools bam2fq unmapped.bam > unmapped.fq
# unmapped.bam을 경유하지 않고 원본 BAM file에서 직접 추출하기
$ samtools bam2fq –f 0x04 BL21.bam > unmapped.fq

unmapped.fq 파일에는 쌍을 이룬 read와 그렇지 않은 것들이 한데 섞여있는 상태이다. 만약 pair를 이룬 read와 그렇지 않은 것(singleton)을 별도의 파일로 분리하고 싶다면 khmer package의 split-paired-reads.py를 이용하면 된다.

$ split-paired-reads.py -0 orphan.fq -1 unmmaped_1.fq -2 unmapped_2.fq unmapped.fq

Fastq 파일 내의 read 순서가 뒤섞였거나 pair 관계가 깨진 것을 고치려면 BBmap package의 repair.sh를 사용해도 된다. 다음 두 명령어의 결과는 같다.

$ repair.sh in=unmapped.fq out=fixed.fq outs=singletons.fq repair
$ bbsplitpairs.sh in=unmapped.fq out=fixed.fq outs=singletons.fq fint

Unmapped read pair를 추출하려면 어떻게 하는 것이 좋을까? 바로 위에 보인 사례에서는 일단 모든 unmapped read를 뽑아낸 다음 이를 수선하여 pair를 이룬 것만을 최종적으로 얻었다. 다른 방법으로는 BAM 파일에서부터 unmapped read pair를 직접 꺼내는 것도 가능하다. 0x04(segment unmmaed)와 0x08(next segment in the template unmapped)를 합쳐서 십진수 12(0x0c)를 flag로 주면 read와 그 mate가 전부 unmapped 상태인 것을 골라서 추출하게 된다.

$ samtools view –bf 12 BL21.bam > unmapped_pair.bam
$ samtools bam2fq -1 paired1.fq -2 paired2.fq -0 /dev/null –s /dev/null unmapped_pair.bam
# 또는 중간 BAM 파일을 경유하지 않고 한번에 실행
$ samtools bam2fq -1 paired1.fq -2 paired2.fq -0 /dev/null –s /dev/null –f 12 BL21.bam

자료에 따라서는 -f 12 -F 256으로 설정하기도 한다. 256(0x100)은 not primary alignment라는 뜻이므로 -F 256은 이를 부정하여 primary alignment가 된다. 그러나 unmapped read에게는 이것이 의미가 없고, 실제로 테스트한 결과 –f 12 단독으로 실행하거나 –f 12 –F 256 복합으로 실행한 경우에 차이는 없었다.

Read unmapped에 대한 SAM flag이 4(samtools view -f 0x04)이므로 이 조건에 반대되는 모든 read(-F 0x04)가 mapped read에 해당한다. 그러나 properly mapped reads, 즉 mapping된 mate의 거리와 간격이 라이브러리의 크기의 평균적 분포를 만족하는 것을 고르려면 -f 0x03(0x01 for read paired; 0x02 for read mapped in proper pair; paired read의 경우 두 flag는 동시에 쓰여야 함)을 적용해야 한다. 후속 분석 작업의 목적에 따라서는 mate 중 어느 하나만 mapping이 되었어도 read pair를 전부 추출하고 싶은 경우도 있을 것이다. SAM 및 BAM 필터링의 사례에 대해서는 David Fredman의 웹페이지를 참조하라.

여러 조건을 이용하여 복수의 BAM 파일을 추출하였다면 samtools merge 명령을 이용하여 하나로 합친 뒤 samtools bam2fq 명령을 이용하여 fastq 파일로 추출하면 된다. 위의 표에서 보인 사례에서 [5]와 [6]는 paired read까지 만들어 놓으면 같은 결과가 되므로 어떤 flag을 쓰든 관계가 없다. 다음의 사례는 BAM 파일을 탐색하여 mate 중 어느 한쪽만 mapping이 되었더라도 read pair를 추출하여 한 쌍의 fastq 파일로 출력하는 방법을 보이고 있다.

$ samtools view -b -F 12 BL21.bam > mapped_paired.bam
$ samtools view -b -f 4 -F 8 BL21.bam > unmap_map.bam
$ samtools view -b -f 8 -F 4 BL21.bam > map_unmap.bam
$ samtools merge at_least_one_mapped.bam mapped_paired.bam unmap_map.bam map_unmap.bam mapped_paired.bam
# Read name을 기준으로 정렬을 해야(-n) 나중에 read pair를 별도의 파일로 잘 분리해 낼 수 있다.
$ samtools sort -n at_least_one_mapped.bam -o at_least_one_mapped_sorted.bam
# -N 옵션은 read name에 /1 또는 /2를 추가한다.
$ samtools bam2fq –N -1 paired_1.fq -2 paired_2.fq at_least_one_mapped_sorted.bam

추출할 read의 목록을 먼저 만든 다음에 원본 fastq 파일로부터 추출해 내는 방법도 가능하다. 이때에는 seqtk가 필요하다.

$ samtools view -F 12 BL21.bam | cut -f1 | sort | uniq > list_both
# 다음 명령어에서는 -f 4 -F 8로 실행해도 무방하다.
$ samtools view -f 8 -F 4 BL21.bam | cut -f1 | sort | uniq > list_single
$ cat list_both list_single > list_all
$ cp list_all list_1
$ cp list_all list_2
$ sed -i 's/$/\/1/' list_1
$ seqtk subseq reads_1.fq list_1 > pairedNew_1.fq
$ sed -i 's/$/\/2/' list_2
$ seqtk subseq reads_2.fq list_2 > pairedNew_2.fq

본 장의 시작 부분에서 fastq-dump를 이용한 SRA data 다운로드 방법을 간략하게 설명하였다. 만약 SRA와 연계된 메타데이터 파일이 필요하거나 웹브라우저 환경의 Run Selector를 이용해야 하는 경우, 또는 아마존 웹 서비스(Amazon Web Service, AWS)를 통한 다운로드가 필요하다면 NCBI의 공식 문서인 Download SRA sequences from Entrez search results를 참조하라.

Single run에 대한 데이터를 열람하여 다운로드하려면 우선 Run Browser에 하나의 Run accession(예: SRR2014531)을 입력하여 데이터에 접근한 다음, Reads 탭을 클릭한다. Filter 창에 조건을 입력하거나 혹은 빈 상태로 Filtered Download를 클릭한 뒤, 이어서 나타나는 창에서 Download Format을 설정한 뒤 Download를 클릭하면 된다. HTTP를 통한 다운로드이므로 한번에 하나의 데이터만 받을 수 있다.

복수의 SRA Experiment(예: https://www.ncbi.nlm.nih.gov/sra/SRP116709)로 이루어진 데이터의 경우 원하는 Experiment에만 체크를 하여 오른쪽 위의 ‘Send to’를 클릭하여 나타나는 Run Selector로 데이터를 보낸다. 모든 Experiment를 선택하려면 체크를 하지 않은 채로 진행하면 된다. SRA Run Selector에 나타나는 Select 패널에서 Accession List를 클릭하면 SRR_Acc_List.txt 파일을 다운로드할 수 있다. 이 파일에 포함된 Run Accession을 사용하여 .sra 파일을 일괄적으로 다운로드한 뒤(SRA-Toolkit의 prefetch 사용) fastq-dump를 이용하여 fastq로 최종 전환하는 방법은 Download_From_SRA 문서에 설명이 되어있다. fastq-dump는 (1) NCBI에서 직접 데이터 파일(.sra)을 정해진 위치에 다운로드한 뒤 fastq로 전환하는 것과 (2) 미리 받아놓은 .sra 파일에서 fastq 파일을 추출하는 것 전부 가능하다. 다음의 예제에서는 GNU parallel을 사용하여 SRR_Acc_List.txt 파일에 수록된 Run accession을 명령어로 공급하는 기법도 쓰였다. –j 2 옵션을 주면 동시에 두 개의 커맨드를 실행하게 된다. Prefech는 ~/ncbi/public/sra/에 모든 .sra 파일을 저장하며, fastq-dump는 이로부터 fastq 파일을 추출하여 현재 위치의 fastq 디렉토리에 이를 저장한다.

$ PATH=/usr/local/apps/sratoolkit.2.9.6-centos_linux64/bin/:$PATH
$ parallel -j 1 prefetch {} ::: $(cat SRR_Acc_List.txt)
$ parallel -j 1 fastq-dump --skip-technical -F --split-files -O fastq {} ::: $(cat SRR_Acc_List.txt) 

SRA의 Run Browser에서 Data access 탭을 선택한 다음 URL을 클릭하여 직접 다운로드를 할 수도 있다. 아래 그림처럼 NCBI Location 오른쪽의 URL을 클릭하면 SRR8981517라는 파일을 받게 된다.

이를 ~/ncbi/public/sra/SRR8981517.sra라는 이름으로 저장한 뒤 다음 명령어를 써서 fastq-dump를 실행한다.

$ fastq-dump --split-files SRR8981517
Read 1431064 spots for SRR8981517
Written 1431064 spots for SRR8981517

사실은 앞서 설명했듯이 ‘fastq-dump --split-files SRR8981517’라고만 입력하여 실행을 해도 .sra 파일 다운로드와 fastq 추출 작업이 자동으로 연이어서 진행된다.

Sequencing depth and coverage: key consideration in genomic analysis. Nature Reviews Genetics 15, 121–132 (2014)