User Tools

Site Tools


to_be_renamed

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
to_be_renamed [2018/10/30 13:36]
hyjeong [Roary]
to_be_renamed [2019/06/27 11:24] (current)
hyjeong [LS-BSR]
Line 45: Line 45:
 {{:​roary.png?​600|Roary의 pan-genome construction. Supplementary Data file에서 발췌.}} {{:​roary.png?​600|Roary의 pan-genome construction. Supplementary Data file에서 발췌.}}
  
-일반적으로 pan genome을 계산하기 위하여 all-against-all BLAST를 거치므로 대단히 많은 시간이 소요되지만,​ Roary는 CD-HIT을 사용하여 partial sequence를 사전에 제거하고 pre-clustering을 실시하여 후속 all-against-all 검색에 사용할 서열집합의 수를 크게 줄이므로 다른 도구에 비하여 훨씬 빠르게 분석을 수행하게 된다. 계산이 끝난 뒤 union/​intersection/​difference 산출을 도와주는 ​a명령어(query_pan_genome) 등 다양한 유틸리티를 제공한다는 점이 특징이다.+일반적으로 pan genome을 계산하기 위하여 all-against-all BLAST를 거치므로 대단히 많은 시간이 소요되지만,​ Roary는 CD-HIT을 사용하여 partial sequence를 사전에 제거하고 pre-clustering을 실시하여 후속 all-against-all 검색에 사용할 서열집합의 수를 크게 줄이므로 다른 도구에 비하여 훨씬 빠르게 분석을 수행하게 된다. 계산이 끝난 뒤 union/​intersection/​difference 산출을 도와주는 명령어(query_pan_genome) 등 다양한 유틸리티를 제공한다는 점이 특징이다.
  
 설치가 약간 까다로운 편이지만([[https://​github.com/​sanger-pathogens/​Roary/​blob/​master/​README.md|설치 가이드]] 요즘은 Virtual Machine과 Docker image를 제공하고 있다. Roary의 기본적인 사용법을 알아보자. 필요한 입력물은 GFF3 파일이 유일하다. 유전체 서열만 준비된 경우라면 [[http://​www.vicbioinformatics.com/​software.prokka.shtml|Prokka]]를 이용하여 genome annotation을 하면 되고, NCBI에서 다운로드한 GenBank 파일이 있따면 bp_genbank2gff3.pl(BioPerl에 포함)를 이용하여 GFF3 file을 만들면 된다. GenBank file에서 /pseudo로 표시된 유전자는 초기 단계에서 아예 단백질 서열이 추출되지 않으므로 계산에서 제외된다. 설치가 약간 까다로운 편이지만([[https://​github.com/​sanger-pathogens/​Roary/​blob/​master/​README.md|설치 가이드]] 요즘은 Virtual Machine과 Docker image를 제공하고 있다. Roary의 기본적인 사용법을 알아보자. 필요한 입력물은 GFF3 파일이 유일하다. 유전체 서열만 준비된 경우라면 [[http://​www.vicbioinformatics.com/​software.prokka.shtml|Prokka]]를 이용하여 genome annotation을 하면 되고, NCBI에서 다운로드한 GenBank 파일이 있따면 bp_genbank2gff3.pl(BioPerl에 포함)를 이용하여 GFF3 file을 만들면 된다. GenBank file에서 /pseudo로 표시된 유전자는 초기 단계에서 아예 단백질 서열이 추출되지 않으므로 계산에서 제외된다.
Line 51: Line 51:
   $ roary *.gff # core gene의 multiple sequence alignment를 하지 않고 pan genome 분석(기본 동작)   $ roary *.gff # core gene의 multiple sequence alignment를 하지 않고 pan genome 분석(기본 동작)
   $ roary -e --mafft -p 8 -i 90 -f output_dir -z *.gff # 본문 설명 참조   $ roary -e --mafft -p 8 -i 90 -f output_dir -z *.gff # 본문 설명 참조
-세 번째 명령에서는 core gene의 alignment를 mafft로 실행하고(-e --mafft) 8 개의 thread를 사용하며(-p 8) blast percent identity cutoff를 90%(기본은 95%)로 한다는 것(-i 90)이다. -o output_dir을 설정하지 않으면 현재 디렉토리에 결과 파일을 작성한다. -e라고만 하면 prank를 사용하여 core gene multiFASTA alignment를 하게 되는데 codon-aware alignment를 하므로 mafft보다는 정확하지만 느리다. -z는 중간 결과 파일(alignment)을 output_dir/​pan_genome_sequences 디렉토리에 남겨둔다는 뜻이다. _로 시작하는 10여 개 중간 결과 파일도 지우지 않고 남기며, 현 디렉토리에는 각 유전체에서 추출한 단백질 서열 파일을 담은 임시 디렉토리도 남기게 된다.+세 번째 명령에서는 core gene의 alignment를 mafft로 실행하고(-e --mafft) 8 개의 thread를 사용하며(-p 8) blast percent identity cutoff를 90%(기본은 95%)로 한다는 것(-i 90)이다. -e 옵션만 주면 PRANK에 의한 core gene alignment를 실시한다(좀 느리다). -o output_dir을 설정하지 않으면 현재 디렉토리에 결과 파일을 작성한다. -e라고만 하면 prank를 사용하여 core gene multiFASTA alignment를 하게 되는데 codon-aware alignment를 하므로 mafft보다는 정확하지만 느리다. -z는 중간 결과 파일(alignment)을 output_dir/​pan_genome_sequences 디렉토리에 남겨둔다는 뜻이다. _로 시작하는 10여 개 중간 결과 파일도 지우지 않고 남기며, 현 디렉토리에는 각 유전체에서 추출한 단백질 서열 파일을 담은 임시 디렉토리도 남기게 된다.
  
 === Input file에 대한 주의사항 === === Input file에 대한 주의사항 ===
Line 57: Line 57:
 {{:​noaccession.png?​500|}} {{:​noaccession.png?​500|}}
  
-그러므로 PGAP annotation file은 그대로 사용하지 말고 [[fillAccession.pl]]으로 처리하여 accession을 채운 뒤 GFF로 전환해야 한다. 또한 RAST server에서 export한 GFF3 파일도 Roary에서 그대로 쓰일 수가 없다. 왜냐하면 염기서열이 뒷부분에 있지 않기 때문이다. 뿐만 아니라 gene 없이 cds feature만 있다는 것도 문제가 된다. 몇 가지를 테스트해 본 경험으로 가장 바람직한 것은, RAST에서 export한 GFF3 file에서 CDS feature만 골라낸 것 + '##​FASTA'​ 라인 + contig sequence 파일을 합쳐서 새로 만든 GFF3 파일을 사용하는 것이 좋다. ​+그러므로 PGAP annotation file은 그대로 사용하지 말고 [[fillAccession.pl]]으로 처리하여 accession을 채운 뒤 GFF로 전환해야 한다. ​ 
 + 
 +또한 RAST server에서 export한 GFF3 파일도 Roary에서 그대로 쓰일 수가 없다. 왜냐하면 염기서열이 뒷부분에 있지 않기 때문이다. 뿐만 아니라 gene 없이 cds feature만 있다는 것도 문제가 된다. 몇 가지를 테스트해 본 경험으로 가장 바람직한 것은, RAST에서 export한 GFF3 file에서 CDS feature만 골라낸 것 + '##​FASTA'​ 라인 + contig sequence 파일(유전자 염기서열 파일이 아님!)을 합쳐서 새로 만든 GFF3 파일을 사용하는 것이 좋다. '​Name='​도 '​Product='​로 바꾸는 것을 강력 권장한다. 왜냐하면 이것이 gene id로 쓰이게 되면 중간에 공백이 들어 있어서 나중에 매우 불편해지기 때문이다.
  
 === Output files === === Output files ===
Line 93: Line 95:
 {{ :​svg.png?​400 |}} {{ :​svg.png?​400 |}}
  
 +=== Strain-specific gene 찾기 ===
 +R에서 gene_presence_absence.Rtab과 gene_presence_absence.csv 두 파일을 다루면 된다. 다음의 사례에서는 Lb_1-46 균주에서 특이적인 유전자의 id를 추출하는 사례를 보여준다. Lb_1-46은 데이터프레임으로 읽어들이면 Lb_1.46으로 바뀌는 것에 유의해라. 구글링을 잘 하면 이를 원래 이름 그대로 유지하는 방법이 있다([[https://​blog.genoglobe.com/​2018/​12/​r-readtable.html|하루에 한 R 관련글]]).
 +  > dat = read.table("​gene_presence_absence.Rtab",​sep="​\t",​header=T,​row.names=1) ​
 +  > dat$Lb_1.46
 +  > dat.s = subset(dat, rowSums(as.matrix(dat))==1)
 +  > Lb.s = dat.s[which(dat.s$Lb_1.46==1),​]
 +  > Lb.s.genes = rownames(Lb.s)
 +  > dat.2 = read.table("​gene_presence_absence.csv",​sep=",",​header=T,​row.names=1)
 +  > dat.2$Lb_1.46
 +  > dat.f = dat.2[which(rownames(dat.2) %in% Lb.s.genes),​]
 +  > dat.f$Lb_1.46
 +  > write.table(dat.f$Lb_1.46,"​out.txt",​sep="​\t",​quote=F,​col.names=F)
 === 기타 해결할 문제 === === 기타 해결할 문제 ===
 결과 파일을 열어보면 일부 단백질의 ID가 변형되어 쓰인 것을 알 수 있다. 즉 원본 annotation file에서 MT_RS20470라는 locus tag을 갖는 유전자가 "​MT_RS20470.p01_16540"​으로 바뀐 것이다. .p01은 유전자에서 단백질로 번역됨을 나타내는 것이겠지만 "​_16540"​은 무엇인가?​ (원래 '​_'​가 3회 반복된 것이나 DokuWiki에서 그렇게 타이핑하면 긴 가로선이 나오기 때문에 부득이하게 하나만 타이핑하였다) 결과 파일을 열어보면 일부 단백질의 ID가 변형되어 쓰인 것을 알 수 있다. 즉 원본 annotation file에서 MT_RS20470라는 locus tag을 갖는 유전자가 "​MT_RS20470.p01_16540"​으로 바뀐 것이다. .p01은 유전자에서 단백질로 번역됨을 나타내는 것이겠지만 "​_16540"​은 무엇인가?​ (원래 '​_'​가 3회 반복된 것이나 DokuWiki에서 그렇게 타이핑하면 긴 가로선이 나오기 때문에 부득이하게 하나만 타이핑하였다)
 +
 +Prokka에서 만든 gff3 파일을 사용하였더니 tRNA gene을 결과물 중에 포함시키는 현상이 가끔 관찰된다.
  
 === Roary 이후 개발된 프로그램 === === Roary 이후 개발된 프로그램 ===
Line 140: Line 156:
 LS-BSR의 간단한 사용법을 알아보자. genome sequence 파일은 확장자가 .fasta가 아니면 작동을 하지 않는다. markers.fasta는 항상 유전자 LS-BSR의 간단한 사용법을 알아보자. genome sequence 파일은 확장자가 .fasta가 아니면 작동을 하지 않는다. markers.fasta는 항상 유전자
  ​염기서열 파일이어야 한다. 검색용 프로그램은 기본이 tblastn이며(단백질로 번역을 먼저 함), 만약 blastn을 쓰고 싶으면 -b blastn을 더해야 한다. -x test라고 지정하면 test라는 디렉토리를 임시로 생성하여 중간 계산 파일을 넣은 뒤, 최종  ​염기서열 파일이어야 한다. 검색용 프로그램은 기본이 tblastn이며(단백질로 번역을 먼저 함), 만약 blastn을 쓰고 싶으면 -b blastn을 더해야 한다. -x test라고 지정하면 test라는 디렉토리를 임시로 생성하여 중간 계산 파일을 넣은 뒤, 최종
- ​결과 파일은 '​test_'​로 시작하는 접두사를 갖고 만들어진다. 분석이 끝나면 test 디렉토리는 없어진다. /home/hyjeong/github/LS-BSR을 사용해라. --- //​[[hyjeong@kribb.re.kr|Haeyoung Jeong]] ​2018/04/26 15:52//+ ​결과 파일은 '​test_'​로 시작하는 접두사를 갖고 만들어진다. 분석이 끝나면 test 디렉토리는 없어진다. /usr/local/apps/LS-BSR/ls_bsr.py을 사용해라. 비록 conda를 이용하여 LS-BSR을 설치했다 하여도 이는 environment만 구성할 뿐, 실제 작동 스크립트 등은 git로 깔았기 때문이다. --- //​[[hyjeong@kribb.re.kr|Haeyoung Jeong]] ​2019/06/27 11:22//
  
   # gene screen method (peptides in interest are ready)   # gene screen method (peptides in interest are ready)
-  # /​data/​apps/​LS-BSR/​ls_bsr.py <= 스크립트 위치 
   $ python path_to_LS-BSR/​ls_bsr.py -d genomes -g markers.fasta -x test   $ python path_to_LS-BSR/​ls_bsr.py -d genomes -g markers.fasta -x test
   # de novo gene prediction method   # de novo gene prediction method
   $ python path_to_LS-BSR/​ls_bsr.py -d genomes -c usearch -x test   $ python path_to_LS-BSR/​ls_bsr.py -d genomes -c usearch -x test
  
-gene prediction method에서는 prodigal로 유전자를 예측한 뒤 단백질로 번역한 다음 usearch로 클러스터링(default:​ -i 0.9)하여 사용한다. 클러스터링 방법은 usearch/​vsearch/​cd-hit이 가능하다($PATH에 있어야 함).+계산이 끝나면 현 디렉토리에 markers.fasta가 번역된 genes.pep 파일이 생긴다. ​gene prediction method에서는 prodigal로 유전자를 예측한 뒤 단백질로 번역한 다음 usearch로 클러스터링(default:​ -i 0.9)하여 사용한다. 클러스터링 방법은 usearch/​vsearch/​cd-hit이 가능하다($PATH에 있어야 함).
  
  
to_be_renamed.1540874197.txt.gz · Last modified: 2018/10/30 13:36 by hyjeong