Genome Informatics Laboratory at KRIBB

User Tools

Site Tools

Trace:
bioinfo:roary

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
bioinfo:roary [2022/06/22 16:47] – [cluster ID를 이용하여 GFF 파일에서 지정된 유전자의 염기서열 추출하기] hyjeongbioinfo:roary [2023/06/26 20:02] (current) – [결과물 이해하기] hyjeong
Line 33: Line 33:
  
 ==== 다른 방법으로 마련한 통해 마련한 GFF3을 쓰면 안되나? ==== ==== 다른 방법으로 마련한 통해 마련한 GFF3을 쓰면 안되나? ====
-GFF3 파일은 이를 만들어내는 프로그램에 따라서 형식이 조금씩 다르다. Roary가 사용하는 GFF3 파일은 뒷부분 '##FASTA'라는 선언을 시작으로 전체 염기서열과 각 단백질의 서열이 포함되는 형태이다. Prokka가 만든 GFF 파일은 큰 문제 없이 사용 가능하지만 NCBI에서 직접 내려받은 GFF 파일은 곤란하다. 또한 PGAP(local)이 만든 GenBank 파일을 GFF3로 전환하여 그대로 쓰면 일부 필드가 비어 있어서 에러가 발생할 것이다. 다음은 로컬 머신에서 PGAP을 실행하여 얻은 GenBank 파일의 앞부분 사례이다. Accession과 version 필드가 빈 상태이기 때문에 이를 그대로 GFF3로 전환한 뒤 Roary를 실행하면 제대로 결과가 나오지 않는다.+GFF3 파일은 이를 만들어내는 프로그램에 따라서 형식이 조금씩 다르다. Roary가 사용하는 GFF3 파일은 feature에 대한 정보를 지나쳐서 뒷부분의 '##FASTA'라는 선언을 시작으로 전체 염기서열과 각 단백질의 서열이 포함되는 형태이다. Prokka가 만든 GFF 파일은 큰 문제 없이 사용 가능하지만 NCBI에서 직접 내려받은 GFF 파일은 곤란하다. 또한 PGAP(local)이 만든 GenBank 파일을 GFF3로 전환하여 그대로 쓰면 일부 필드가 비어 있어서 에러가 발생할 것이다. 다음은 로컬 머신에서 PGAP을 실행하여 얻은 GenBank 파일의 앞부분 사례이다. Accession과 version 필드가 빈 상태이기 때문에 이를 그대로 GFF3로 전환한 뒤 Roary를 실행하면 제대로 결과가 나오지 않는다.
   LOCUS       chromosome           3744397 bp    DNA     circular BCT 15-MAR-2022   LOCUS       chromosome           3744397 bp    DNA     circular BCT 15-MAR-2022
   DEFINITION  Ralstonia pseudosolanacearum strain SL1931 chromosome, complete   DEFINITION  Ralstonia pseudosolanacearum strain SL1931 chromosome, complete
Line 115: Line 115:
 [[https://sanger-pathogens.github.io/Roary/|Roary 공식 웹사이트]]의 <output files> 항목에 꽤 상세한 설명이 나온다. 여기에서는 공식 문서에서도 해소되지 않는 궁금증을 풀어 나가도록 한다.  [[https://sanger-pathogens.github.io/Roary/|Roary 공식 웹사이트]]의 <output files> 항목에 꽤 상세한 설명이 나온다. 여기에서는 공식 문서에서도 해소되지 않는 궁금증을 풀어 나가도록 한다. 
  
-clustered_proteins는 가장 이해하기 쉬운 파일이다. summary_statistics.txt에 나타난 total genes는 결국 gene cluster(singleton 포함)이며, 각 클러스터를 구성하는 균주별 유전자의 목록이 clustered_proteins 파일에 라인 단위로 보여진다. 클러스터의 대표 서열은 pan_genome_reference.fa 파일에 수록된다. 10개의 genome 중에서 어떤 기준으로 다음의 것이 선정되었는지는 나도 잘 모른다. 아마 Roary 논문에 이에 대한 설명이 있을 것이다. 다음의 사례에서 서열 ID는 실제로 선정된 유전자이고, description 항목의 'dna'는 cluster ID에 해당한다.+**clustered_proteins**는 가장 이해하기 쉬운 파일이다. **summary_statistics.txt**에 나타난 total genes는 결국 gene cluster(singleton 포함)이며, 각 클러스터를 구성하는 균주별 유전자의 목록이 clustered_proteins 파일에 라인 단위로 보여진다. 클러스터의 대표 서열은 **pan_genome_reference.fa** 파일에 수록된다. 10개의 genome 중에서 어떤 기준으로 다음의 것이 선정되었는지는 나도 잘 모른다. 아마 Roary 논문에 이에 대한 설명이 있을 것이다. 다음의 사례에서 서열 ID는 실제로 선정된 유전자이고, description 항목의 'dna'는 cluster ID에 해당한다.
   >f62c63acf1630b40d04ccf54718a60ec_3 dnaA   >f62c63acf1630b40d04ccf54718a60ec_3 dnaA
   GTGGACAGCCATACCTCTGAACTATGGCAGCAAATTCTATCCATTATACAAACCAAGCTG   GTGGACAGCCATACCTCTGAACTATGGCAGCAAATTCTATCCATTATACAAACCAAGCTG
Line 314: Line 314:
  
 === 종합 실습 === === 종합 실습 ===
-Core gene 목록을 뽑은 뒤 이에 해당하는 유전자 염기서열을 GCF_000520775.1.gff에서 뽑아내어 FASTA file로 만든다. 작업 디렉토리는 fixed_input_files로 한다. 가장 먼저 할 일은 뽑아낼 유전자(core gene)의 cluster ID를 추출하는 것이다.+Core gene 목록을 뽑은 뒤 이에 해당하는 유전자 염기서열을 **GCF_000520775.1.gff**에서 뽑아내어 FASTA file로 만든다. 작업 디렉토리는 fixed_input_files로 한다. 가장 먼저 할 일은 뽑아낼 유전자(core gene)의 cluster ID를 추출하는 것이다.
   $ query_pan_genome -a intersection -g ../clustered_proteins *gff   $ query_pan_genome -a intersection -g ../clustered_proteins *gff
   $ awk '{print $1}' pan_genome_results | sed 's/://' | sort > core_genes.txt   $ awk '{print $1}' pan_genome_results | sed 's/://' | sort > core_genes.txt
Line 323: Line 323:
   - gene_presence_absence.csv는 콤마로 각 컬럼을 구분하고 있으며, 컬럼은 따옴표로 둘러싸여 있다. Annotation 컬럼은 내부적으로 콤마를 포함할 수도 있다. 많은 유전자의 product 이름에는 콤마가 들어있는 경우가 많다.   - gene_presence_absence.csv는 콤마로 각 컬럼을 구분하고 있으며, 컬럼은 따옴표로 둘러싸여 있다. Annotation 컬럼은 내부적으로 콤마를 포함할 수도 있다. 많은 유전자의 product 이름에는 콤마가 들어있는 경우가 많다.
   - Paralog가 출현하면 하나의 컬럼에 두 개 이상의 유전자가 있을 수 있다.    - Paralog가 출현하면 하나의 컬럼에 두 개 이상의 유전자가 있을 수 있다. 
-1, 즉 따옴표로 둘러싼 컬럼을 콤마로 구분하는 csv 파일을 파싱하는 것은 매우 흔한 일이면서도 의외로 까다롭다. [[https://github.com/BurntSushi/xsv|xsv]], [[http://lorance.freeshell.org/csv/|AWK CSV parser]], [[https://csvkit.readthedocs.io/en/latest/index.html|csvkit]] 등 다양한 솔루션이 있는데, 나는 가장 마지막의 csvkit을 써 보기로 하였다. gene_presence_absence.csv의 1번 컬럼('Gene')이 core_genes.txt와 겹치는 라인의 17번 컬럼('GCF_000520775.1')을 추출하되 하나의 core gene에 대하여 유전자가 1개씩 정확히 있는지를 특히 눈여겨 보아야 한다.+1, 즉 따옴표로 둘러싼 컬럼(특히 그 내부에 콤마를 포함하는 경우)을 콤마로 구분하는 csv 파일을 파싱하는 것은 매우 흔한 일이면서도 의외로 까다롭다(([[https://stackoverflow.com/questions/71543341/how-can-i-parse-csv-files-with-quoted-fields-containing-commas-in-awk|How can I parse CSV files with quoted fields containing commas, in awk?]]에 의하면 awk -v FPAT='([^,]+)|(\"[^\"]+\")' file.csv 명령을 써서 필드 구분자가 아니라 필드 자체를 매치하도록 만들면 된다.)). [[https://github.com/BurntSushi/xsv|xsv]], [[http://lorance.freeshell.org/csv/|AWK CSV parser]], [[https://csvkit.readthedocs.io/en/latest/index.html|csvkit]] 등 다양한 솔루션이 있는데, 나는 가장 마지막의 csvkit을 써 보기로 하였다. gene_presence_absence.csv의 1번 컬럼('Gene')이 core_genes.txt와 겹치는 라인의 17번 컬럼('GCF_000520775.1')을 추출하되 하나의 core gene에 대하여 유전자가 1개씩 정확히 있는지를 특히 눈여겨 보아야 한다.
   $ csvcut -c 1,17 ../gene_presence_absence.csv | sort -t, -k1,1 > gene_presence_absence_GCF_000520775.1.csv   $ csvcut -c 1,17 ../gene_presence_absence.csv | sort -t, -k1,1 > gene_presence_absence_GCF_000520775.1.csv
   $ join -t, -1 1 -2 1 core_genes.txt gene_presence_absence_GCF_000520775.1.csv > core_gene_info_GCF_000520775.1.txt   $ join -t, -1 1 -2 1 core_genes.txt gene_presence_absence_GCF_000520775.1.csv > core_gene_info_GCF_000520775.1.txt
-GCF_000520775.1_info.txt 파일을 콤마로 구분했을 때 두 번째 컬럼에 유전자가 정확히 1개씩 있는지, 혹은 탭을 구분자로 하여 2개 이상이 있는지를 확인해 둔다.  +**core_gene_info_GCF_000520775.1.txt** 파일을 콤마로 구분했을 때 두 번째 컬럼에 유전자가 정확히 1개씩 있는지, 혹은 탭을 구분자로 하여 2개 이상이 있는지를 확인해 둔다.  
-== Headline ==+== [1] pan_genome_sequences 디렉토리의 core gene alignment에서 염기서열을 뽑아내기 == 
 +  # 현 위치는 fixed_input_files
   $ cp core_gene_info_GCF_000520775.1.txt ../pan_genome_sequences/   $ cp core_gene_info_GCF_000520775.1.txt ../pan_genome_sequences/
   $ cd ../pan_genome_sequences/    $ cd ../pan_genome_sequences/ 
Line 335: Line 336:
   > seqkit grep -r -p $b $a.fa.aln | seqret fasta::stdin fasta:stdout | sed "s/>/>$a /" >> core_genes_from_GCF_000520775.1.fa   > seqkit grep -r -p $b $a.fa.aln | seqret fasta::stdin fasta:stdout | sed "s/>/>$a /" >> core_genes_from_GCF_000520775.1.fa
   > done   > done
-== Headline ==+염기서열이 추출된 유전자 수는 하나가 적다. core_gene_info_GCF_000520775.1.txt에는 abc-f라는 core gene이 있지만 정작 pan_genome_sequences 디렉토리에 있는 파일은 abc_f.fa.aln이다. 파일명에는  
 +-'를 쓰지 않게 만드는 것 같다. 다소 불편하지만 실행 전후의 유전자 수를 점검해 보는 것이 좋을 것이다. core_genes_from_GCF_000520775.1.fa 내의 gap을 없내는 것은 EMBOSS의 degapseq 유틸리티를 이용한다. 
 +  $ degapseq core_genes_from_GCF_000520775.1.fa core_genes_from_GCF_000520775.1.fa.nogaps 
 +== [2] GFF 파일과 이로부터 추출한 FASTA file을 이용하는 방법 ==
  
 +  # 현 위치는 fixed_input_files
 +  $ perl test_20220622.pl core_gene_info_GCF_000520775.1.txt GCF_000520775.1.gff > core_genes_from_GCF_000520775.1.fa
 +Perl 스크립트 'test_20220622.pl'는 그렇게 효율적이지는 못하다. EMBOSS의 seqret 명령어를 추출할 core gene 수만큼 실행하기 때문이다. GCF_000520775.1.fna 파일을 스크립트 실행 개시 부분에서 Bio::Seq 개체로 한 번에 읽어들인 뒤, core gene이 위치한 서열 개체 정보를 이용하여 출력하는 것이 더 빠를 수도 있다. 별로 신통치 않은 스크립트이지만 소개하여 본다.
  
 +  #!/usr/bin/perl
 +  #
 +  
 +  open LIST, $ARGV[0];
 +  while ( <LIST> ) {
 +          chomp;
 +          @data = split /,/, $_;
 +          $seen{$data[1]} = $data[0];
 +  }
 +  
 +  ( $strain = $ARGV[1] ) =~ s/\.gff$//;
 +  open GFF, $ARGV[1];
 +  while ( <GFF> ) {
 +          chomp;
 +          @temp = split /\t/, $_;
 +          next unless $temp[2] eq 'CDS';
 +          $temp[8] =~ m/ID=([^;]+);/;
 +          $op = '';
 +          if ( exists $seen{$1} ) {
 +                  $seq = $strain . '.fna:' . $temp[0];
 +                  $op = '-sreverse' if $temp[6] eq '-';
 +                  $idndesc = '>' . $seen{$1} . ' ' . $1;
 +                  system( "seqret -sbegin $temp[3] -send $temp[4] $op $seq fasta::stdout | sed \"s/>.*/$idndesc/\" " );
 +          }
 +  }
 +두 번째 방법으로 추출한 FASTA 파일에는 gap이 존재하지 않으며, 오로지 위치를 이용하여 염기서열을 추출하게 되므로 core gene의 group ID가 약간 틀려도 문제를 일으키지 않는다는 장점이 있다. 
  
 ===== Roary에게 경의를! Scory ===== ===== Roary에게 경의를! Scory =====
bioinfo/roary.1655884042.txt.gz · Last modified: 2022/06/22 16:47 by hyjeong