====== Reference mapping 및 variant calling ======
===== 개요 =====
  * Reference로는 //E. coli// B str. REL606의 유전체 서열을 NCBI에서 다운로드하여 사용한다.
  * Short read mapper로는 bowtie2 사용(PATH 환경변수에 설정)
  * samtools, bctools 및 vcfutils.pl([[http://www.htslib.org/|링크]])는 편의상 A5-miseq에 포함된 것을 사용하고, 역시 PATH 환경변수에 포함시킨다.

  $ export PATH=$PATH:/usr/local/Bio/bowtie2-2.2.6
  $ export PATH=$PATH:/usr/local/Bio/a5_miseq_linux_20140604/bin

===== 1. Bowtie2를 이용한 mapping =====
  * [[http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml|bowtie2 manual page]]
==== Reference sequence의 다운로드 및 인덱스 생성 ====
  $ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000017985.1_ASM1798v1/GCA_000017985.1_ASM1798v1_genomic.fna.gz
  $ gzip -d GCA_000017985.1_ASM1798v1_genomic.fna.gz
  $ mkdir bowtie2 (반드시 별도의 디렉토리를 만들지 않아도 됨)
  $ bowtie2-build GCA_000017985.1_ASM1798v1_genomic.fna bowtie2/REL606
  
==== Mapping ====
만약 paired read의 방향이 일반적인 fr(forward-reverse)가 아니고 rf(from mate-pair library) 혹은 ff라면 어떻게 할 것인가? 또한 library insert size를 mapping 단계에서 지정하려면 어떻게 해야 할 것인지 생각해 보자. SAM 및 BAM 파일에 대한 지식이 필요하면 [[https://github.com/davetang/learning_bam_file|Learning BAM files 사이트]]를 둘러볼 것을 권장한다.

  $ bowtie2 -p 16 -x bowtie2/REL606 -1 BL21-20x_1.fastq -2 BL21-20x_2.fastq -S BL21.sam
  $ samtools view -b -S -o BL21.bam BL21.sam
  
===== (optional) Read alignment의 시각화 =====
  * [[http://www.broadinstitute.org/igv/|IGV]]나 [[https://ics.hutton.ac.uk/tablet/|tablet]]을 사용하여 assembly 및 alignment의 시각화 가능
  * [IGV] Genomes -> Create .genome File에서 fasta file을 로드한 다음 File -> Load from file에서 sorting 및 indexing이 완료된 bam file을 로드

===== (optional) Mapping 관련 수치의 출력 =====
  $ samtools flagstat BL21.bam
  $ samtools view –f 0x04 –c BL21.bam | wc -l (unmapped read의 수)
  $ samtools view –F 0x04 –c BL21.bam | wc -l (mapped 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의 측정)

===== (optional) 각 contig의 average read depth 산출 =====
bam 파일로 소팅한 다음 samtools depth 명령을 실행하면 reference 서열의 모든 염기 위치에 대한 read depth가 출력된다. 더욱 상세한 정보를 원한다면 samtools mpileup 명령을 쓰면 되나, 이번 작업에서는 그렇게까지는 할 필요가 없다. 다만 reference sequence 단위의 집계를 위해서 간단한 스크립트 [[read_depth.pl]]이 필요하다. samtools 1.4.1에서는 sort와 index에 관련한 명령어 문법이 약간 바뀌었으므로 [[http://www.htslib.org/doc/samtools-1.4.1.html|매뉴얼 페이지]]를 참조하기 바란다.

  $ samtools sort -o BL21_sorted.bam BL21.bam
  $ samtools depth BL21_sorted.bam (출력물이 어떠한지 관찰해 보기)
  $ samtools depth BL21_sorted.bam | read_depth.pl > depth.csv

===== 2. Variant calling =====
==== [1] 일반적인 과정 ====
  $ samtools faidx GCA_000017985.1_ASM1798v1_genomic.fna
  $ samtools sort -o BL21.sorted.bam BL21.bam
  $ samtools index BL21.sorted.bam
  $ samtools mpileup -u -f GCA_000017985.1_ASM1798v1_genomic.fna BL21.sorted.bam > BL21.bcf
  $ bcftools view -v -c -g BL21.bcf > BL21.vcf
  $ vcfutils.pl varFilter -D100 BL21.vcf > BL21.vcf.filtered (filter: max read depth 100)
==== [2] Pipe를 사용한 간략화 ====
  $ samtools mpileup -uf GCA_000017985.1_ASM1798v1_genomic.fna BL21.sorted.bam | bcftools view -bvcg - > BL21.bcf
  $ bcftools view BL21.bcf | vcfutils.pl varFilter -D100 > BL21.vcf.filtered
====== 참고자료 ======
  * SAM/BAM 파일의 이해와 sorting/indexing에 관해서는 [[http://homer.salk.edu/homer/basicTutorial/samfiles.html|Understanding and manipulating SAM/BAM files]]을 참조할 것
  * [SAMtools 공식 문서] [[http://samtools.sourceforge.net/mpileup.shtml|Calling SNPs/InDels with SAMtools/BCFtools]]
  * [[https://gist.github.com/davfre/8596159|SAM and BAM filtering onliners]] 
  * getinsertsize.py [[https://gist.github.com/davidliwei/2323462|다운로드 링크]] [[http://allaboutbioinfo.blogspot.kr/2012/04/estimating-paired-end-read-insert.html|Estimating paired-end read insert size from BAM/SAM files]]