BAM file을 정렬 및 인덱싱한 뒤 bcftools(SAMtools와 같이 설치됨)을 실행한다. BAM file은 sort가 먼저 되어 있어야 인덱싱이 가능하다.
$ samtools faidx GCF_000017985.1_ASM1798v1_genomic.fna $ samtools index BL21_sorted.bam $ bcftools mpileup -f GCF_000017985.1_ASM1798v1_genomic.fna BL21.bam | bcftools call --ploidy 1 -mv -Ob -o calls.bcf
BCF 파일의 필터링 기법은 Variant calling and filtering 문서를 참조하라.
$ bcftools view -i '%QUAL>=20' calls.bcf | bcftools stats $ bcftools filter -sLowQual -g3 -G10 -e'%QUAL<10 || (RPB<0.1 && %QUAL<15) || (AC<2 && %QUAL<15)' calls.bcf > calls.bcf.filtered.vcf
염기서열 변이가 실제로 유전자 산물의 기능에 어떠한 영향을 미치는지를 예측하려면 SnpEff를 활용한다. 위 단계에서 만든 VCF 파일과 reference database가 필요하다. Ensemble genome database에서 제공하는 레퍼런스 자료를 이용하지 못하는 경우를 대비하여 read mapping에 사용하였던 유전체 파일과 annotation file(GenBank or GFF3)를 사용하여 직접 데이터베이스를 만들어야 한다. /usr/local/apps/snpEff/snpEff.config를 현재 디렉토리에 복사한 뒤 data.dir을 쓰기 권한이 있는 디렉토리로 변경해 놓는다(예: ~/snpEff/data/). 다음으로는 매핑에 사용한 FASTA file과 짝을 이루는 GenBank file을 입수하여 ~/snpEff/data/REL606/genes.gbk로 복사해 둔다(reference는 적당한 문자열로 설정한다). snfEff.config의 맨 끝에 다음을 추가한다.
REL606.genome : REL606 $ java -jar /usr/local/apps/snpEff/snpEff.jar build -genbank -v REL606
위 단계에서 생성한 VCF 파일의 첫번째 필드에 보이는 서열 ID에서 버전 번호(.1)를 제거한다. REL606의 경우는 다음과 같이 한다.
$ sed -e 's/^NC_012967.1/NC_012967/' calls.bcf.filtered.vcf > calls.bcf.filtered.vcf.mod
VCF 파일에 대한 annotation을 실시한다.
$ java -jar /usr/local/apps/snpEff/snpEff.jar REL606 calls.bcf.filtered.vcf.mod > calls.bcf.filtered.vcf.mod.ann
현 디렉토리에 생기는 snpEff_genes.txt 및 snpEff_summary.html 파일도 참조한다.
Picard는 대용량 시퀀싱 데이터와 SAM/BAM/CRAM 및 VCF 포맷의 데이터의 조작을 위한 명령행 도구이다. Broad Institute에서 개발한 GATK는 variant discovery와 genotyping에 초점을 맞춘 도구로서 사용자의 최종 목적에 맞는 best practice를 제공한다.
Snippy는 박테리아 유전체에 최적화된 variant calling 및 core genome alignment 도구이다. 주어진 reference genome sequence에 대하여 NGS sequence read를 매핑하고(bwa mem 사용) 이로부터 변이를 탐색하여주며, core SNP 영역까지 추출해 주는 매우 편리한 pipeline이다. 비교 자료가 contig로 주어지면 read simulation을 먼저 실시하여 이를 사용하여 분석에 들어가게 된다. 생성된 alignment 정보는 phylogenetic tree를 만드는데 활용할 수 있다.
한 쌍의 fastq file을 이용하여 snippy를 실행하는 사례는 다음과 같다.
$ snippy --cpus 16 --outdir mysnps --ref REL606.gbk --R1 SRR2014531_1.fastq --R2 SRR2014531_2.fastq
mysnps/snps.log 파일을 열어보면 실제 어떤 명령이 수행되었는지를 볼 수 있다.
samtools faidx reference/ref.fa bwa index reference/ref.fa mkdir -p reference/genomes && cp -f reference/ref.fa reference/genomes/ref.fa ln -sf reference/ref.fa . ln -sf reference/ref.fa.fai . mkdir -p reference/ref && gzip -c reference/ref.gff > reference/ref/genes.gff.gz snpEff build -c reference/snpeff.config -dataDir . -gff3 ref bwa mem -Y -M -R '@RG\tID:mysnps\tSM:mysnps' -t 16 reference/ref.fa /nas/project/MT03_manuals/snippy/SRR2014531_1.fastq /nas/project/MT03_manuals/snippy/SRR2014531_2.fastq | samclip --max 10 --ref reference/ref.fa.fai | samtools sort -n -l 0 -T /tmp/snippy.21328. --threads 16 -m 500M | samtools fixmate -m - - | samtools sort -l 0 -T /tmp/snippy.21328. --threads 16 -m 500M | samtools markdup -T /tmp/snippy.21328. -r -s - - > snps.bam samtools index snps.bam fasta_generate_regions.py reference/ref.fa.fai 151839 > reference/ref.txt freebayes-parallel reference/ref.txt 16 -p 2 -P 0 -C 10 --min-repeat-entropy 1.5 --strict-vcf -q 13 -m 60 --min-coverage 10 -F 0.05 -f reference/ref.fa snps.bam > snps.raw.vcf bcftools view --include 'FMT/GT="1/1" && QUAL>=100 && FMT/DP>=10 && (FMT/AO)/(FMT/DP)>=0' snps.raw.vcf | vt normalize -r reference/ref.fa - | bcftools annotate --remove '^INFO/TYPE,^INFO/DP,^INFO/RO,^INFO/AO,^INFO/AB,^FORMAT/GT,^FORMAT/DP,^FORMAT/RO,^FORMAT/AO,^FORMAT/QR,^FORMAT/QA,^FORMAT/GL' > snps.filt.vcf snpEff ann -noLog -noStats -no-downstream -no-upstream -no-utr -t -c reference/snpeff.config -dataDir . ref snps.filt.vcf > snps.vcf snippy-vcf_to_tab --gff reference/ref.gff --ref reference/ref.fa --vcf snps.vcf > snps.tab snippy-vcf_extract_subs snps.filt.vcf > snps.subs.vcf bcftools convert -Oz -o snps.vcf.gz snps.vcf bcftools index -f snps.vcf.gz bcftools consensus -f reference/ref.fa -o snps.consensus.fa snps.vcf.gz bcftools convert -Oz -o snps.subs.vcf.gz snps.subs.vcf bcftools index -f snps.subs.vcf.gz bcftools consensus -f reference/ref.fa -o snps.consensus.subs.fa snps.subs.vcf.gz rm -f snps.subs.vcf.gz snps.subs.vcf.gz.csi snps.subs.vcf.gz.tbi
동일한 reference에 대하여 여러 쌍의 sequence read(single end read나 assembled contig도 무관)를 반복하여 분석한 뒤 core genome SNP alignment를 추출해야 하는 경우에는 snippy-multi를 이용하는 것이 편리하다. Tab-delimited file(‘input.tab’)에 샘플과 read file(혹은 contig sequence file) 정보를 수록한 뒤 snippy-multi를 실행하면 실행 가능한 스크립트가 생성된다. Input.tab 파일의 사례는 다음과 같다. 파일은 압축 여부에 상관없이 사용할 수 있다.
# input.tab = ID R1 [R2] Isolate1 /path/to/R1.fq.gz /path/to/R2.fq.gz Isolate1b /path/to/R1.fastq.gz /path/to/R2.fastq.gz Isolate1c /path/to/R1.fa /path/to/R2.fa # single end reads supported too Isolate2 /path/to/SE.fq.gz Isolate2b /path/to/iontorrent.fastq # or already assembled contigs if you don't have reads Isolate3 /path/to/contigs.fa Isolate3b /path/to/reference.fna.gz
Input.tab 파일을 인수로 사용하여 snippy-multi를 실행하여 다음과 같이 실행용 스크립트 runme.sh를 만들고 이를 별도로 실행하면 된다. 아래의 실행 사례는 Genetic characterization and comparison of Clostridium botulinum isolates from botulism cases in Japan between 2006 and 2011 (Applied and Environmental Microbiology 2014) 논문에서 언급한 10개 균주의 일루미나 시퀀싱 자료(SRA DRA002253)를 이용한 것이다. Reference는 ATCC 3502(NC_009495)를 사용하였다.
$ snippy-multi input.tab --ref reference.gbk --cpus 8 > runme.sh $ sh ./runme.sh …(중간 생략)… Generating core.full.aln Creating TSV file: core.txt Running: snp-sites -c -o core.aln core.full.aln Have a suggestion? Tell me at http://github.com/tseemann/snippy/issues Done. $ cat core.txt # alignment/core-size statistics ID LENGTH ALIGNED UNALIGNED VARIANT HET MASKED LOWCOV Af650 4155278 4047142 106918 223 68 0 1150 Aichi2011 4155278 3438770 602285 62463 1247 0 112976 Ehime2011 4155278 581072 1113528 12161 119 0 2460559 Fukuoka2010 4155278 3556374 505508 5832 2243 0 91153 Ibaraki2007 4155278 3443909 515634 5822 767 0 194968 Iwate2007 4155278 3555971 497487 5825 3201 0 98619 Iwate2008 4155278 3631759 500401 5821 4033 0 19085 Miyagi2006 4155278 2572053 637580 6420 991 0 944654 Miyagi2006-01 4155278 1408615 757847 6215 335 0 1988481 NCTC2916 4155278 2880942 519180 5299 747 0 754409 Okayama2011 4155278 3009491 1062954 21982 3860 0 78973 Tochigi2008 4155278 3529946 592005 5285 1369 0 31958 Reference 4155278 4155278 0 0 0 0 0
Runme.sh 스크립트에서는 주어진 reference에 대하여 각 샘플을 snippy로 분석한 다음 모든 결과물에 대하여 snippy-core 명령을 실행하게 만든다. 바로 이 명령이 core genome SNP alignment를 산출하는 것이다. Snippy-core의 결과 파일은 전부 'core'라는 접두사를 갖는다. 정렬된 core genome의 FASTA 파일은 core.full.aln에 저장되고, 이를 이용하여 계통수를 그릴 수 있다. 이 과정을 실제로 수행하려면 alignment에 대한 후처리가 필요하다. Snippy의 부속 스크립트인 snippy-clean_full_align은 비정상적인 문자를 'N'으로 치환하고, gubbins는 recombination이 일어난 영역을 제거한다. Core.full.aln 파일에는 모든 taxon이 같은 염기를 갖는 invariant site를 포함하므로, snp-sites를 이용하여 whole-genome alignment에서 SNP가 존재하는 부분만을 추출하여 데이터의 분량을 줄일 수 있다1). Gubbins는 염기 치환 밀도가 높은 곳을 재조합이 일어난 영역으로 간주하여 이를 FASTA alignment에서 제거함은 물론, 재조합이 일어나지 않은 영역을 바탕으로 RAxML에 의한 트리 파일('.final_tree.tre')도 생성해 준다. 재조합 영역의 예측 결과는 EMBL과 GFF 포맷으로 별도 제공한다.
Core genome alignment를 후처리하여 계통수를 작성하는 방법은 다음과 같다.
$ snippy-clean_full_aln core.full.aln > clean.full.aln $ conda activate nanopore # gubbins 실행을 위해 conda nanopore 환경으로 진입 (nanopore) $ run_gubbins.py -p gubbins clean.full.aln # -p <prefix> (nanopore) $ conda deactivate # conda base 환경으로 돌아옴 $ snp-sites -c gubbins.filtered_polymorphic_sites.fasta > clean.core.aln $ fasttree -gtr -nt clean.core.aln > clean.core.tre
Run_grubbins.py를 거치고 나면 일부 샘플이 제거되는 일이 벌어지기도 한다. 이는 gap이 너무 많은 taxon을 제거하는 동작 때문일 것이다(–filter_percentage <%> 또는 -f <%>으로 설정; 기본 수치는 25). 따라서 같은 종에 속하는 균주라 하더라도 유전적 거리가 먼 관계에 있는 것은 gubbins로 처리하지 않는 것이 좋을 것이다. 다시 말해서 reference sequence에 대하여 샘플의 read를 매핑하였을 때 커버되지 않는 영역이 25%를 넘을 정도로 멀다면2), 이는 gubbins를 이용하여 recombination 영역을 검출하는 분석 작업의 의미가 없어진다.
Gubbins의 분석 결과는 Phandango 웹서버에 직접 드래그하여 넣으면 tree 구조와 recombination 위치를 같이 시각화할 수 있다. Phandango는 세균의 집단유전체 분석 결과를 대화식으로 표현하는 유용한 도구이다.
Reference와 이리 조립된 염기서열을 비교하여 변이를 탐색하는 방법에 대해서는 별도로 작성한 문서 Whole-genome alignment를 읽어보기 바란다.
박테리아의 유전체가 갖는 고유 특성(mosaic) 때문에 core genome을 기반으로 하는 SNP 분석으로는 만족할 만한 결과를 얻기 어렵다. 2021년 Genome Biology에 발표된 Pandora: nucleotide-resolution bacterial pan-genomics with reference graphs라는 논문에서는 시퀀싱된 genome을 reference의 재조합체로 근사하여 새로운 변이를 탐지하며, 여러 샘플을 pan-genotype으로 분류한다. 이 논문에서 활용한 varifier는 VCF로 주어지는 변이 발굴 결과를 평가하거나(subcommand 'vcf_eval'), 또는 VCF를 생성하는(subcommand 'make_truth_vcf') 도구이다.
make_truth_vcf는 두 개의 유전체 염기서열 파일(G1 & G2)을 인수로 갖는다. dnadiff와 minimap2/paftools를 각각 사용하여 pairwise SNP set를 생성한 뒤 합쳐서 거른다. 각 allele에 대하여 G1으로부터 좌우50 염기씩을 포함하는 probe를 만든 뒤 G2에 대하여 매핑한 뒤 mapping quality가 0이 되는 것은 paralog/duplicate/repeat에서 유래한 것으로 간주하여 제거한다. 마지막으로 매핑 후 allele 내부의 미스매치를 점검하여 확정한다.
usage: varifier make_truth_vcf [options] <truth_fasta> <ref_fasta> <outdir>