GATK 사용기 2

이제 GATK 로 넘어간다.
지금부터 할 일은 snp 라고 추측되는 부분을 중심으로 reads를 재분석하여 다시 realignment 해야 된다. 왜냐하면 이 부분에서 mapping 이 잘못되는 경우가 많기 때문. 그래서 이미 알려진 snp db 를 가지고 다시 분석한다.

GATK의 help 로 보면 입이 떡 벌어질 길고도 복잡한 페이지가 보인다. 읽는 것을 포기한다.

java -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R reference.fa -I dedup_reads.bam -L 20 -known gold_indels.vcf -o realignment_targets.list

홈페이지에 나온 간략본이다.

-T - GATK 의 수많은 기능 중 이번에 사용할 기능 RealignerTargetCreator
-R - reference sequence
-I - picard 로 duplicated reads 를 marking 하고 난 뒤의 bam file
-L - 특정 contig 를 대상으로 분석하고 싶을 때 쓰는 옵션.
-known - 이미 알려진 snp 가 있다면 찾아서 넣어주면 좀 더 잘 찾는다고함. reference sequence 와 같은 contig name 을 가지고 있어야 한다. 좀 더 구분하기 쉬우라고 'chr01' 이런 식으로 고쳤다가 피봤다.
-o - output file

java -jar GenomeAnalysisTK.jar -T IndelRealigner -R reference.fa -I dedup_reads.bam -targetIntervals realignment_targets.list -known gold_indels.vcf -o realigned_reads.bam

대부분의 옵션은 위와 겹치니까 설명은 안 하겠다. (귀찮음.. ;o;)
위에서 구한 list 와 SNPdb 를 가지고 realign 을 수행하여 bam file 을 다시 만든다.


계속해서 GATK 를 사용해보자.

java -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R referencefa -I realigned_reads.bam -L 20 -knownSites dbsnp.vcf -knownSites gold_indels.vcf -o recal_data.table


java -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R reference.fa -I realigned_reads.bam -L 20 -knownSites dbsnp.vcf -knownSites gold_indels.vcf -BQSR recal_data.table -o post_recal_data.table

java -jar GenomeAnalysisTK.jar -T AnalyzeCovariates -R reference.fa -L 20 -before recal_data.table -after post_recal_data.table -plots recalibration_plots.pdf

java -jar GenomeAnalysisTK.jar -T PrintReads -R reference.fa -I realigned_reads.bam -L 20 -BQSR recal_data.table -o recal_reads.bam

picard 와 GATK 사용하기

GATK 를 사용해서 snp 분석을 하기 전에 picard 를 써서 sam files 를 정리해야 된다. 이와중에 시행착오가 참 많았다..

먼저 만약 1개의 sample 로 여러 개의 NGS data를 만들었을 경우, picard 의 MergeSamFiles 명령어로 여러 sam files 합치는 것과 동시에 coordinate sorting 하자.
어차피 sam file 은 reads 순서로 정렬이 되어 있으므로 reference mapping 결과대로 다시 sorting 해야 되니까 이참에 하자.

java -jar -Xmx128g picard.jar MergeSamFiles I=1.sam I=2.sam O=filename.merge.sam SO=coordinate USE_THREADING=true TMP_DIR=/.../ R=reference.fasta

합치고 싶은 sam files 각각을 I 옵션에 넣어주면 된다.
USE_THREADING - multi thread 를 활성화 시켜줘서 20% 정도 속도 향상을 시켜준다고 하는데 딱히 체감되지는 않지만...
TMP_DIR - sorting 할 때 디스크를 많이 쓴다. 여분의 디스크를 따로 설치하고 그곳을 지정해주는게 속도 향상에 더 큰 영향을 주는 듯. samtools 는 메모리 영역을 크게 할당할 수 있었는데 picard 는 그런 옵션이 없는 듯.
R - reference 파일인데 딱히 필요하지는 않는다.
* -Xms32g, -Xmx128g - java 옵션이다. Xms 는 최소 메모리, Xmx 는 최대 메모리를 지정한다. 이걸 크게 설정해 놓으면 disk I/O 작업이 빈번하게 일어나는 것을 막아 disk drive 가 고장날 확률을 줄여주고 작업 속도도 빨라진다.... 라고 기대했는데 아닌가 보다. TMP 디렉토리에 잔뜩 만들어지는 걸 보아하니..

굳이 합칠 일이 없을 경우 sorting 하는 방법이다.

java -jar /share/picard-2.0.1/dist/picard.jar SortSam I=filename.sam O=filename_sorted.sam SO=coordinate TMP_DIR=/.../tmp/ R=reference.fasta

Sorting 이 끝나면 이제 duplication 제거를 한다. 이것은 sequencing 될 때 적당한 강도의 signal 확보를 위해 PCR 을 하게 되고 이 과정에서 특정 DNA 만 유독 강하게 나올 수 있다. 이것을 찾아 가장 좋은 것 하나로 만드는 작업이다. 그렇지 않으면 depth 계산이 엉뚱하게 나올 수 있다.

java -jar picard.jar MarkDuplicatesWithMateCigar I=filename.sorted.sam O=filename.sorted.dedup.bam M=filename.sorted.dedup.bam.metrics.txt TMP_DIR=/.../tmp/ R=/reference.fasta

다음 글에서 GATK 로 넘어간다..... 그냥 넘어가면 에러 뜬다.. ㅡ_ㅡ...
bam file 의 index 를 만들어야 된다.

java -jar picard.jar BuildBamIndex I=filename.sorted.dedup.bam O=filename.sorted.dedup.bai TMP_DIR=/.../tmp/ R=/reference.fasta

굳이 output 옵션을 쓸 필요는 없다.. 알아서 현재 디렉토리에 .bai 로 나온다.

하나 더 해야 된다.. 아.. 짜증.. 진작 말하지..

java -jar picard.jar CreateSequenceDictionary R=reference.fasta O=refrence.dict

쉽게 말해 GATK 에서 사용할 reference fasta 의 index 를 생성하는거다.

bowtie2 와 bwa mem 명령어

GATK 를 사용해 snp 분석을 하려고 했더니 첩첩산중... 너무 어렵다.
알고 보니 sam 파일을 만들 때부터 뭔가 잔뜩 들어간다. 이제부터 그 이야기를 써보자.

GATK 에서 error 가 발생하는 원인 중 하나는 multiple match..
정확한 분석을 위해서는 하나의 match 만이 인정되어야 한다. bwa 의 경우 -M 옵션을 붙이면 해결이 된다.
bowtie2 의 경우는 잘 모르겠다.

또다른 error 의 원인 중 하나는 sam header 를 잘 써줘야 한다. bwa 의 경우 -R '@RG\tID:text\tSM:text\tPL:text'

@RG - sam header 중 RG 를 쓰겠다는 뜻..
\t - tab seperate
ID:text - ID section 에 text 를 집어넣겠다. RG header 를 쓸 때 필수.
DS:text - 설명문 넣겠음.
LB:text - library 이름을 넣겠음.
PL:text - Platform/technology 넣는 곳으로 GATK 에서 필수. CAPILLARY, LS454, ILLUMINA, SOLID, HELICOS, IONTORRENT, ONT, PACBIO 로 선점되어 있으니 맞는 거 골라써라.
PM:text - PL 과 비슷하나 좀더 상세 모델명을 쓸 수 있다. 자유롭게 쓰면 된다.
SM:text - sample name

bwa mem 명령을 이용해 paired end fastq 파일을 sam 으로 만들어 보자.
bwa mem -t 4 -M -R '@RG\tID:rice\tLB:PE1\tPL:ILLUMINA\tPM:HISEQ2000\tSM:BJJNo1' reference.fasta paired1.fq paired2.fq >filename.sam

bowtie2 를 사용해보자.
bowtie2 -q --phred33 --rg-id rice --rg "LB:PE1" --rg "PL:ILLUMINA" --rg "PM:HISEQ2000" --rg "SM:BJJNo1" --fr -p 4 -x reference.fasta -1 paired1.fq -2 paired2.fq -S filename.sam