• Reference로는 E. coli B str. REL606의 유전체 서열을 NCBI에서 다운로드하여 사용한다.
• Short read mapper로는 bowtie2 사용(PATH 환경변수에 설정)
• samtools, bctools 및 vcfutils.pl(링크)는 편의상 A5-miseq에 포함된 것을 사용하고, 역시 PATH 환경변수에 포함시킨다.

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

• bowtie2 manual page

$ 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

$ bowtie2 -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

• IGV나 tablet을 사용하여 assembly 및 alignment의 시각화 가능
• [IGV] Genomes → Create .genome File에서 fasta file을 로드한 다음 File → Load from file에서 sorting 및 indexing이 완료된 bam file을 로드

$ 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의 측정)

$ samtools faidx GCA_000017985.1_ASM1798v1_genomic.fna
$ samtools sort BL21.bam BL21.sorted
$ samtools sort BL21.bam BL21.sorted
$ 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)

$ 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에 관해서는 Understanding and manipulating SAM/BAM files을 참조할 것
• [SAMtools 공식 문서] Calling SNPs/InDels with SAMtools/BCFtools
• SAM and BAM filtering onliners
• getinsertsize.py 다운로드 링크 Estimating paired-end read insert size from BAM/SAM files