• 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

만약 paired read의 방향이 일반적인 fr(forward-reverse)가 아니고 rf(from mate-pair library) 혹은 ff라면 어떻게 할 것인가? 또한 library insert size를 mapping 단계에서 지정하려면 어떻게 해야 할 것인지 생각해 보자. SAM 및 BAM 파일에 대한 지식이 필요하면 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

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

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

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

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

$ 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