$ export PATH=$PATH:/usr/local/Bio/bowtie2-2.2.6 $ export PATH=$PATH:/usr/local/Bio/a5_miseq_linux_20140604/bin
$ 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
$ 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