Genome Informatics Laboratory at KRIBB

User Tools

Site Tools

Trace:
to_be_renamed

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
to_be_renamed [2018/10/30 14:10] – [Roary] hyjeongto_be_renamed [2022/06/18 13:12] (current) – [Roary] hyjeong
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 61: Line 61:
 또한 RAST server에서 export한 GFF3 파일도 Roary에서 그대로 쓰일 수가 없다. 왜냐하면 염기서열이 뒷부분에 있지 않기 때문이다. 뿐만 아니라 gene 없이 cds feature만 있다는 것도 문제가 된다. 몇 가지를 테스트해 본 경험으로 가장 바람직한 것은, RAST에서 export한 GFF3 file에서 CDS feature만 골라낸 것 + '##FASTA' 라인 + contig sequence 파일(유전자 염기서열 파일이 아님!)을 합쳐서 새로 만든 GFF3 파일을 사용하는 것이 좋다. 'Name='도 'Product='로 바꾸는 것을 강력 권장한다. 왜냐하면 이것이 gene id로 쓰이게 되면 중간에 공백이 들어 있어서 나중에 매우 불편해지기 때문이다. 또한 RAST server에서 export한 GFF3 파일도 Roary에서 그대로 쓰일 수가 없다. 왜냐하면 염기서열이 뒷부분에 있지 않기 때문이다. 뿐만 아니라 gene 없이 cds feature만 있다는 것도 문제가 된다. 몇 가지를 테스트해 본 경험으로 가장 바람직한 것은, RAST에서 export한 GFF3 file에서 CDS feature만 골라낸 것 + '##FASTA' 라인 + contig sequence 파일(유전자 염기서열 파일이 아님!)을 합쳐서 새로 만든 GFF3 파일을 사용하는 것이 좋다. 'Name='도 'Product='로 바꾸는 것을 강력 권장한다. 왜냐하면 이것이 gene id로 쓰이게 되면 중간에 공백이 들어 있어서 나중에 매우 불편해지기 때문이다.
  
 +Roary를 실행하면 십중팔구는 다음과 같은 메시지와 함께 GFF 파일을 수정하여 fixed_input_files 디렉토리에 복사하게 된다. 나중에 query_pan_genome 스크립트로 GFF 파일을 대상으로 하는 작업을 할 때에는 수정된 것을 써야 한다.
 +2021/06/30 15:53:28 Input file contains duplicate gene IDs, attempting to fix by adding a unique suffix, new GFF in the fixed_input_files directory: 
 === Output files === === Output files ===
 {{:output.png?600|test}} {{:output.png?600|test}}
Line 80: Line 82:
  
 === Command line tool의 사용법 === === Command line tool의 사용법 ===
-테스트를 해 보았는데 그 동작이 완벽한 것 같지는 않다.+원본 GFF 파일이 아니라 roary가 수정한 것(fixed_input_files/ 디렉토리에 있는 것)을 사용해야 한다. 인수로 주어질 수 있는 것은 *gff 또는 단일 GFF 파일이다.  
   $ query_pan_genome -h # for help   $ query_pan_genome -h # for help
   $ query_pan_genome -a union *.gff # 결과물: pan_genome_results   $ query_pan_genome -a union *.gff # 결과물: pan_genome_results
-  $ query_pan_genome -a intersection *.gff # 결과물: pan_genome_results (텅 비었음. why?)+  $ query_pan_genome -a intersection *.gff # 결과물: pan_genome_results
   # difference   # difference
   # complement   # complement
Line 95: Line 97:
 {{ :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 142: Line 158:
 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.1540876246.txt.gz · Last modified: (external edit)