Genome Informatics Laboratory at KRIBB

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Microbial variant calling ======
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 파일의 필터링 기법은 [[https://samtools.github.io/bcftools/howtos/variant-calling.html|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
  
===== Variant annotation =====
염기서열 변이가 실제로 유전자 산물의 기능에 어떠한 영향을 미치는지를 예측하려면 [[https://pcingola.github.io/SnpEff/|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와 GATK(Genome Analysis Tool Kit) =====
[[https://broadinstitute.github.io/picard/|Picard]]는 대용량 시퀀싱 데이터와 SAM/BAM/CRAM 및 VCF 포맷의 데이터의 조작을 위한 명령행 도구이다. Broad Institute에서 개발한 [[https://software.broadinstitute.org/gatk/|GATK]]는 variant discovery와 genotyping에 초점을 맞춘 도구로서 사용자의 최종 목적에 맞는 best practice를 제공한다. 


===== Snippy를 이용한 간편한 변이 탐색 =====

[[https://github.com/tseemann/snippy|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 ([[https://journals.asm.org/doi/10.1128/AEM.02134-14|Applied and Environmental Microbiology 2014]]) 논문에서 언급한 10개 균주의 일루미나 시퀀싱 자료(SRA [[https://www.ncbi.nlm.nih.gov/sra/?term=DRA002253|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'으로 치환하고, [[https://github.com/nickjcroucher/gubbins|gubbins]]는 recombination이 일어난 영역을 제거한다. Core.full.aln 파일에는 모든 taxon이 같은 염기를 갖는 invariant site를 포함하므로, [[https://github.com/sanger-pathogens/snp-sites|snp-sites]]를 이용하여 whole-genome alignment에서 SNP가 존재하는 부분만을 추출하여 데이터의 분량을 줄일 수 있다((snippy-multi가 만든 runme.sh 스크립트를 실행하면 core genome alignment(‘core.full.aln’ 파일)로부터 자동적으로 SNP alignment(‘core.aln’ 파일)를 추출한다. 정확한 계통수를 그리는 것이 목적이라면 gubbins를 사용하여 recombination이 일어난 부분을 제거하는 것이 바람직할 것이다.)). 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%를 넘을 정도로 멀다면((Snippy가 각 샘플에 대하여 작성한 .aligned.fa 파일에서 depth=0이거나(‘-‘로 표시) 0 < depth < -mincov인 영역(‘N’으로 표시)는 core gene alignment 파일로 그대로 저장되고, gubbins에서는 전부 missing data로 취급된다. ‘-mincov’의 기본 수치는 10이다.)), 이는 gubbins를 이용하여 recombination 영역을 검출하는 분석 작업의 의미가 없어진다. 

Gubbins의 분석 결과는 [[http://jameshadfield.github.io/phandango/#/|Phandango 웹서버]]에 직접 드래그하여 넣으면 tree 구조와 recombination 위치를 같이 시각화할 수 있다. Phandango는 세균의 집단유전체 분석 결과를 대화식으로 표현하는 유용한 도구이다.

===== MUMmer(dnadiff) 또는 minimap2을 이용한 변이 탐색 =====
Reference와 이리 조립된 염기서열을 비교하여 변이를 탐색하는 방법에 대해서는 별도로 작성한 문서 [[bioinfo:whole-genome_alignment|Whole-genome alignment]]를 읽어보기 바란다.

===== Varifier =====
박테리아의 유전체가 갖는 고유 특성(mosaic) 때문에 core genome을 기반으로 하는 SNP 분석으로는 만족할 만한 결과를 얻기 어렵다. 2021년 Genome Biology에 발표된 [[https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02473-1|Pandora: nucleotide-resolution bacterial pan-genomics with reference graphs]]라는 논문에서는 시퀀싱된 genome을 reference의 재조합체로 근사하여 새로운 변이를 탐지하며, 여러 샘플을 //pan-genotype//으로 분류한다. 이 논문에서 활용한 [[https://github.com/iqbal-lab-org/varifier|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>