perl -i -pe 's/defined \(%/\(%/' mummerplot
이유는 잘 모르겠지만 우분투 14.04 를 설치하고 나니 mummerplot 에서 계속 에러가 난다. 찾아보니 이걸로 수정하면 된단다.
2016년 12월 7일 수요일
mummperplot 에서 define 에러가 남.
이유는 잘 모르겠지만 우분투 14.04 를 설치하고 나니 mummerplot 에서 계속 에러가 난다. 찾아보니 이걸로 수정하면 된단다.
2016년 11월 18일 금요일
fasta 를 genbank 로 바꾸기..
일하다 보니 정말 이것저것 포멧 바꿔야 될 일이 많다.
fasta file 은 첫 줄이 '>' 로 시작하고 이 부분은 document 로 처리된다. 그 이후 줄부터는 문자만 허용된다. 흔히 DNA 나 protein sequence 를 저장하기 위해 많이 사용된다.
하나의 file 에 하나의 fasta 만 존재할 경우 .ffn 혹은 .fsa 확장자를 사용하고 여러 개가 하나의 파일에 있을 경우, .fna 확장자를 사용하지만 거의 지켜지지 않는다. 확장자가 .fa 일 경우 protein sequence 인 경우이지만 역시나 거의 지켜지지 않는다.
genbank 파일은 NCBI 에서 주로 보여지는 형태로 확장자는 .gb 혹은 .gbk 를 사용한다. 저자, 논문, accession number 등 매우 많은 정보를 같이 저장할 수 있다.
GFF 포멧도 있다. 위의 genbank 포멧에는 종속 관계(?) 를 표시할 수 없다는 단점이 있어서 새롭게 제안된 포멧이라고 한다. 기존 gff 를 거쳐 gff2, gff3 로 발전 중이다. 흔히 gff2 는 fasta file 과 같이 존재해야 되며 gff3 는 gff2 의 뒤에 '###', '###FASTA' 를 더해준 뒤 sequence 를 덧붙여 준 형태이다.
이외에도 embl 포멧도 있지만 난 사용할 줄 모르는 관계로 패스..
emboss 의 seqret 를 사용해 fasta 를 gbk 로 바꿔보자.
seqret -sequence fasta.file -fopenfile1 fasta -outseq output.file -osformat2 gb
-fopenfile1 과 -osformat2 는 각각 input 과 output 할 포멧을 결정해주는 옵션이다. 없으면 fasta 가 디폴트인 듯하다.
gff3 --> gff
gff2 --> gff2
embl --> em
genbank --> gb, refseq
ddbj --> ddbj
refseqp --> refseqp
pir --> nbrf
swissprot --> swiss, sw
로 표기하는 듯 하다.
fasta file 은 첫 줄이 '>' 로 시작하고 이 부분은 document 로 처리된다. 그 이후 줄부터는 문자만 허용된다. 흔히 DNA 나 protein sequence 를 저장하기 위해 많이 사용된다.
하나의 file 에 하나의 fasta 만 존재할 경우 .ffn 혹은 .fsa 확장자를 사용하고 여러 개가 하나의 파일에 있을 경우, .fna 확장자를 사용하지만 거의 지켜지지 않는다. 확장자가 .fa 일 경우 protein sequence 인 경우이지만 역시나 거의 지켜지지 않는다.
genbank 파일은 NCBI 에서 주로 보여지는 형태로 확장자는 .gb 혹은 .gbk 를 사용한다. 저자, 논문, accession number 등 매우 많은 정보를 같이 저장할 수 있다.
GFF 포멧도 있다. 위의 genbank 포멧에는 종속 관계(?) 를 표시할 수 없다는 단점이 있어서 새롭게 제안된 포멧이라고 한다. 기존 gff 를 거쳐 gff2, gff3 로 발전 중이다. 흔히 gff2 는 fasta file 과 같이 존재해야 되며 gff3 는 gff2 의 뒤에 '###', '###FASTA' 를 더해준 뒤 sequence 를 덧붙여 준 형태이다.
이외에도 embl 포멧도 있지만 난 사용할 줄 모르는 관계로 패스..
emboss 의 seqret 를 사용해 fasta 를 gbk 로 바꿔보자.
seqret -sequence fasta.file -fopenfile1 fasta -outseq output.file -osformat2 gb
-fopenfile1 과 -osformat2 는 각각 input 과 output 할 포멧을 결정해주는 옵션이다. 없으면 fasta 가 디폴트인 듯하다.
gff3 --> gff
gff2 --> gff2
embl --> em
genbank --> gb, refseq
ddbj --> ddbj
refseqp --> refseqp
pir --> nbrf
swissprot --> swiss, sw
로 표기하는 듯 하다.
2016년 10월 5일 수요일
GATK 로 variant calling 하기..
설명해 놓은 홈페이지..
https://software.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_tools_walkers_haplotypecaller_HaplotypeCaller.php
GVCF 와 variant-only 인데 GVCF 가 자기들이 새로 개발했는데 좋다는 듯.. ㅡ_ㅡ..
java -jar GenomeAnalysisTK.jar -R refence.fasta -T HaplotypeCaller -I input.bam --emitRefConfidence GVCF --dbsnp known.vcf -o output.snps.indels.g.vcf -variant_index_type LINEAR -variant_index_parameter 128000
* -variant_index_type LINEAR -variant_index_parameter 128000 : 위의 웹에는 설명되어 있지 않지만 5.6 버전 이후라서 그런가.. 이거 없다고 에러 뜬다.
* .g.vcf : 이게 확장자란다. 없다고 큰일 나는 건 아니지만 GATKVCFUtils 가 안된다고 warning 뜬다.
HaplotypeCaller 말고도 UnifiedGenotyper 라고 있는데 HaplotypeCaller 대비 성능이 좀 떨어진다는 듯..
variant-only calling on DNAseq 방법..
java -jar GenomeAnalysisTK.jar -R refence.fasta -T HaplotypeCaller -I input.bam --dbsnp known.vcf -stand_call_conf 30 -stand_emit_conf 10 -o output.raw.snps.indels.vcf
snpEffector 를 사용하기 위해서는 variant-only 방법을 사용해야 한다.
추가1) GATK 3.8 버전에서는 -stand_emit_conf 옵션을 빼라는 에러가 뜬다.
https://software.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_tools_walkers_haplotypecaller_HaplotypeCaller.php
GVCF 와 variant-only 인데 GVCF 가 자기들이 새로 개발했는데 좋다는 듯.. ㅡ_ㅡ..
java -jar GenomeAnalysisTK.jar -R refence.fasta -T HaplotypeCaller -I input.bam --emitRefConfidence GVCF --dbsnp known.vcf -o output.snps.indels.g.vcf -variant_index_type LINEAR -variant_index_parameter 128000
* -variant_index_type LINEAR -variant_index_parameter 128000 : 위의 웹에는 설명되어 있지 않지만 5.6 버전 이후라서 그런가.. 이거 없다고 에러 뜬다.
* .g.vcf : 이게 확장자란다. 없다고 큰일 나는 건 아니지만 GATKVCFUtils 가 안된다고 warning 뜬다.
HaplotypeCaller 말고도 UnifiedGenotyper 라고 있는데 HaplotypeCaller 대비 성능이 좀 떨어진다는 듯..
variant-only calling on DNAseq 방법..
java -jar GenomeAnalysisTK.jar -R refence.fasta -T HaplotypeCaller -I input.bam --dbsnp known.vcf -stand_call_conf 30 -stand_emit_conf 10 -o output.raw.snps.indels.vcf
snpEffector 를 사용하기 위해서는 variant-only 방법을 사용해야 한다.
추가1) GATK 3.8 버전에서는 -stand_emit_conf 옵션을 빼라는 에러가 뜬다.
2016년 9월 19일 월요일
LTS 에 커널 버전 업데이트는 따로 있다고??
우분투 14.04.3 을 서버 시스템에 설치하여 사용 중이었다. 그런데 이게 알고 보니..
14.04.1 --> 3.13
14.04.2 --> 3.16
14.04.3 --> 3.19
14.04.4 --> 4.2
14.04.5 --> 4.4
왼쪽이 14.04 의 버전이고 오른쪽이 커널 버전 업데이트 현황이다. 충격적인건 LTS 로 지원하는 건 3.13 과 4.4 라는 것.. 즉, 3.16, 3.19, 4.2 는 커널 업데이트를 따로 해야만 된다.
못 믿겠다고? 밑에 링크에 들어가봐라.
https://wiki.ubuntu.com/Kernel/Support
https://wiki.ubuntu.com/Kernel/LTSEnablementStack
이게 HWE 뭐시기랑 연관이 되어 있단다. 아, 진심 귀찮음.. 커널만 업데이트 할 수 없을까 찾아봤지만 그냥 통째로 16.04 로 업그레이드 하란다. 이거 또 엉뚱하게 4.6 이나 4.8 로 업그레이도 되면 또 자주 업그레이드 하게 되므로 4.4 로 업그레이드 해야겠다.
방법은 나중에 찾자.
14.04.1 --> 3.13
14.04.2 --> 3.16
14.04.3 --> 3.19
14.04.4 --> 4.2
14.04.5 --> 4.4
왼쪽이 14.04 의 버전이고 오른쪽이 커널 버전 업데이트 현황이다. 충격적인건 LTS 로 지원하는 건 3.13 과 4.4 라는 것.. 즉, 3.16, 3.19, 4.2 는 커널 업데이트를 따로 해야만 된다.
못 믿겠다고? 밑에 링크에 들어가봐라.
https://wiki.ubuntu.com/Kernel/Support
https://wiki.ubuntu.com/Kernel/LTSEnablementStack
이게 HWE 뭐시기랑 연관이 되어 있단다. 아, 진심 귀찮음.. 커널만 업데이트 할 수 없을까 찾아봤지만 그냥 통째로 16.04 로 업그레이드 하란다. 이거 또 엉뚱하게 4.6 이나 4.8 로 업그레이도 되면 또 자주 업그레이드 하게 되므로 4.4 로 업그레이드 해야겠다.
방법은 나중에 찾자.
GATK 사용 중 에러..
GenomeAnalysisTK 에서 RealignerTargetCreator 를 사용하던 중에 뭔지 모를 에러가 났다. 에러 메세지 중에 "unknown index" 라는 말이 있어서 bai files 도 다시 만들어 보고 reference index file 도 다시 만들어 봤지만 도통 해결될 기미가 보이지 않는다..
그러다가 결국 picard 를 2.0.1 --> 2.6.1 로 업그레이드 하고 GATK 도 3.5 --> 3.6 으로 업그레이드 하였다. 덕택에 Oracle JAVA 7 과 8을 혼용하다가 8만 사용 가능하게 되었다.
본래 picard 는 JAVA 8 을 일찌감치 사용하였는데 GATK 가 JAVA 7 을 사용하다가 3.6 버전에서 JAVA 8 로 바뀐 것.. 그런데.. 난 본래 8 이었을텐데.. ;;
여하튼.. 그래도 해결이 안된다.. 사흘 가까이 머리를 쥐어 뜯던 중..
'-known /.../.vcf' 옵션이 보였고 여기에 가보니 idx 파일이 보인다. 이걸 지우고 다시 해보니 잘 된다. 엉엉.... ㅠ_ㅠ
vcf index 파일의 일종인데 GATK 가 실행될 때 자동으로 생성되는 파일이라 그동안 몰랐던 것. 젠장.. 에러 메시지면 에러 메시지 답게 뭐가 잘못된 건지 잘 알려줘야 될 것 아닌가.. ;;
다시 며칠 간 지켜보며 진행하다 보니 어딘가에서 에러가 또 있다. 대충 에러 메시지를 살펴보니 Quality 가 맞지 않는단다.. =_=
예전에 시스템이 두번 뒤집힌 적이 있었는데 그때 원본 파일에 손상이 갔는지 일부 데이터가 깨져 있어서 다시 원본으로 덮어 쓴 적이 있었다. 그럼에도 복구가 안된 듯..
아니면 본래 데이터가 phd 64 와 phd 33 가 혼용이 되어 있어서 phd 64 를 phd 33 으로 바꿔놨는데 이게 잘 안됐나보다. 여기서 쓸 수 있는 방법을 찾아보니...
--fix_misencoded_quality_scores / -fixMisencodedQuals
-allowPotentiallyMisencodedQuals / --allow_potentially_misencoded_quality_scores
두 가지 방법을 안내해주더만.. 그중에서 오른쪽 밑에 옵션을 사용했더니 진행이 되었다.
그러다가 결국 picard 를 2.0.1 --> 2.6.1 로 업그레이드 하고 GATK 도 3.5 --> 3.6 으로 업그레이드 하였다. 덕택에 Oracle JAVA 7 과 8을 혼용하다가 8만 사용 가능하게 되었다.
본래 picard 는 JAVA 8 을 일찌감치 사용하였는데 GATK 가 JAVA 7 을 사용하다가 3.6 버전에서 JAVA 8 로 바뀐 것.. 그런데.. 난 본래 8 이었을텐데.. ;;
여하튼.. 그래도 해결이 안된다.. 사흘 가까이 머리를 쥐어 뜯던 중..
'-known /.../.vcf' 옵션이 보였고 여기에 가보니 idx 파일이 보인다. 이걸 지우고 다시 해보니 잘 된다. 엉엉.... ㅠ_ㅠ
vcf index 파일의 일종인데 GATK 가 실행될 때 자동으로 생성되는 파일이라 그동안 몰랐던 것. 젠장.. 에러 메시지면 에러 메시지 답게 뭐가 잘못된 건지 잘 알려줘야 될 것 아닌가.. ;;
다시 며칠 간 지켜보며 진행하다 보니 어딘가에서 에러가 또 있다. 대충 에러 메시지를 살펴보니 Quality 가 맞지 않는단다.. =_=
예전에 시스템이 두번 뒤집힌 적이 있었는데 그때 원본 파일에 손상이 갔는지 일부 데이터가 깨져 있어서 다시 원본으로 덮어 쓴 적이 있었다. 그럼에도 복구가 안된 듯..
아니면 본래 데이터가 phd 64 와 phd 33 가 혼용이 되어 있어서 phd 64 를 phd 33 으로 바꿔놨는데 이게 잘 안됐나보다. 여기서 쓸 수 있는 방법을 찾아보니...
--fix_misencoded_quality_scores / -fixMisencodedQuals
-allowPotentiallyMisencodedQuals / --allow_potentially_misencoded_quality_scores
두 가지 방법을 안내해주더만.. 그중에서 오른쪽 밑에 옵션을 사용했더니 진행이 되었다.
2016년 2월 5일 금요일
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
지금부터 할 일은 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 를 생성하는거다.
먼저 만약 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 를 생성하는거다.
2016년 2월 3일 수요일
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
알고 보니 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
피드 구독하기:
글 (Atom)