Differences

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

--- to_be_renamed [2018/10/30 08:29] – [Roary] hyjeong
+++ to_be_renamed [2022/06/18 13:12] (current) – [Roary] hyjeong
@@ Line 45: / Line 45: @@
 {{: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로 표시된 유전자는 초기 단계에서 아예 단백질 서열이 추출되지 않으므로 계산에서 제외된다.
@@ Line 51: / Line 51: @@
   $ roary *.gff # core gene의 multiple sequence alignment를 하지 않고 pan genome 분석(기본 동작)
   $ 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에 대한 주의사항 ===
@@ Line 57: / Line 57: @@
 {{: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 파일을 사용하는 것이 좋다. 'Name='도 'Product='로 바꾸는 것을 강력 권장한다. 왜냐하면 이것이 gene id로 쓰이게 되면 중간에 공백이 들어 있어서 나중에 매우 불편해지기 때문이다.
+Roary를 실행하면 십중팔구는 다음과 같은 메시지와 함께 GFF 파일을 수정하여 fixed_input_files 디렉토리에 복사하게 된다. 나중에 query_pan_genome 스크립트로 GFF 파일을 대상으로 하는 작업을 할 때에는 수정된 것을 써야 한다.
+/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.png?600|test}}
@@ Line 63: / Line 68: @@
 == 결과 파일에 대한 설명 ==
   * **gene_presence_absence.csv:** Spreadsheet containing statistics on each group and presence and absence of genes in each isolate
-  * core_gene_alignment.aln: MultiFASTA alignment of the core conserved genes
+  * **core_gene_alignment.aln:** MultiFASTA alignment of the core conserved genes
-  * clustered_proteins: File with 1 cluster per line and sequence identifiers
+  * **clustered_proteins:** File with 1 cluster per line and sequence identifiers. gene_presence_absence.*의 첫번째 컬럼(Gene)과 이 파일의 첫번째 컬럼에 쓰인 식별자가 동일하다. 또한 이 파일에 수록된 모든 cluster의 수(=라인 수)는 summary_statistics.txt 파일 맨 아랫줄의 total genes와 같다.
   * *.Rtab: Tab Files for producing graphs in R
   * *embl: EMBL files for visualizing presence and absence of genes
   * accessory_binary_genes.fa.newick: Newick tree based on presence and absence of genes in accessory genome
-  * pan_genome_reference.fa: A FASTA file with a single nucleotide sequence for every cluster in the pan genome
+  * **pan_genome_reference.fa:** A FASTA file with a single REPRESENTATIVE nucleotide sequence for every cluster in the pan genome (core and accessory)
-  * pangenome_sequences (directory): MultiFASTA alignment files for each cluster
+  * **pangenome_sequences** (directory): MultiFASTA alignment files for each cluster
 Roary에서는 accessory gene의 여부를 이용한 binary tree만을 만들어 주므로, core gene alignment에서 phylogenetic tree를 작성하려면 FastTree 등의 프로그램을 직접 돌려야 한다.계산 시간을 줄이려면 [[https://github.com/sanger-pathogens/snp-sites|snp_site]]([[http://mgen.microbiologyresearch.org/content/journal/mgen/10.1099/mgen.0.000056|논문]])로 필터링을 하여 사용하라. 이 과정에서 recombination은 고려하지 않음에 유의한다.
   $ FastTree –nt –gtr core_gene_alignment.aln > my_tree.newick
@@ Line 77: / Line 82: @@
 === Command line tool의 사용법 ===
-테스트를 해 보았는데 그 동작이 완벽한 것 같지는 않다.
+원본 GFF 파일이 아니라 roary가 수정한 것(fixed_input_files/ 디렉토리에 있는 것)을 사용해야 한다. 인수로 주어질 수 있는 것은 *gff 또는 단일 GFF 파일이다.
   $ query_pan_genome -h # for help
   $ 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
   # complement
@@ Line 92: / Line 97: @@
 {{ :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에서 그렇게 타이핑하면 긴 가로선이 나오기 때문에 부득이하게 하나만 타이핑하였다)
+Prokka에서 만든 gff3 파일을 사용하였더니 tRNA gene을 결과물 중에 포함시키는 현상이 가끔 관찰된다.
 === Roary 이후 개발된 프로그램 ===
@@ Line 139: / Line 158: @@
 LS-BSR의 간단한 사용법을 알아보자. genome sequence 파일은 확장자가 .fasta가 아니면 작동을 하지 않는다. markers.fasta는 항상 유전자
  염기서열 파일이어야 한다. 검색용 프로그램은 기본이 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)
-  # /data/apps/LS-BSR/ls_bsr.py <= 스크립트 위치
   $ python path_to_LS-BSR/ls_bsr.py -d genomes -g markers.fasta -x test
   # de novo gene prediction method
   $ 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에 있어야 함).