gff 파일에서 locus tag과 product만 추출하기

NCBI GenBank에서 유전체 정보를 보면 자주 참고하는 자료가 있고 계속 반복되어 나오지만 실제로는 보지 않는 자료가 있다

 

gbk 형식은 아래와 같이 생겼는데

 

     gene            complement(1616..1951)
                     /locus_tag="MAE_00030"
     CDS             complement(1616..1951)
                     /locus_tag="MAE_00030"
                     /codon_start=1
                     /transl_table=11
                     /product="ferredoxin"
                     /protein_id="BAF99824.1"
                     /translation="MPRITVYGQTITCDRGENLRRILLKHDISLYNGASKLINCRGIG
                     SCGTCAVAIVGEVSAINWQEKARLSLPPHNPDNNRRLACQVKVFGDIEVTKYDGFWGQ
                     GDSVISDQF"

/codon_start, /transl_table, 은 사용할 일이 없고 분석 목적에 따라

complement, /protein_id, /translation 부분도 선택적으로 사용하는 경우가 많다.

 

아마도 가장 많이 보게 되는 것이

/locus_tag

/product

이 두 가지 정보일 것 같다

 

그래서 보기 편하게 아래와 예시와 같이 locus tag과 product 정보만을 나열하려고 한다.

FXO09_19305	leucine-rich repeat domain-containing protein
FXO09_19300	hypothetical protein
FXO09_19295	hypothetical protein
FXO09_19290	AAA family ATPase

 

NCBI GenBank FTP에서 받은 GFF 파일로부터 위의 정보를 추출해내려고 한다

GFF 파일은 아래와 같이 생겼다.

 

CP012375.1	Genbank	gene	4805	4926	.	+	.	ID=gene-amyaer_4374;Name=amyaer_4374;gbkey=Gene;gene_biotype=rRNA;locus_tag=amyaer_4374
CP012375.1	Genbank	rRNA	4805	4926	.	+	.	ID=rna-amyaer_4374;Parent=gene-amyaer_4374;Dbxref=SEED:fig|1126.22.rna.4;gbkey=rRNA;locus_tag=amyaer_4374;product=5S ribosomal RNA
CP012375.1	Genbank	exon	4805	4926	.	+	.	ID=exon-amyaer_4374-1;Parent=rna-amyaer_4374;Dbxref=SEED:fig|1126.22.rna.4;gbkey=rRNA;locus_tag=amyaer_4374;product=5S ribosomal RNA
CP012375.1	Genbank	gene	4972	7596	.	-	.	ID=gene-amyaer_0001;Name=amyaer_0001;gbkey=Gene;gene_biotype=protein_coding;locus_tag=amyaer_0001
CP012375.1	Genbank	CDS	4972	7596	.	-	0	ID=cds-AOC50756.1;Parent=gene-amyaer_0001;Dbxref=GO:0003887,SEED:fig|1126.22.peg.1,NCBI_GP:AOC50756.1;Name=AOC50756.1;gbkey=CDS;locus_tag=amyaer_0001;product=DNA polymerase III alpha subunit;protein_id=AOC50756.1;transl_table=11
CP012375.1	Genbank	gene	7680	8510	.	-	.	ID=gene-amyaer_0002;Name=amyaer_0002;gbkey=Gene;gene_biotype=protein_coding;locus_tag=amyaer_0002
CP012375.1	Genbank	CDS	7680	8510	.	-	0	ID=cds-AOC50757.1;Parent=gene-amyaer_0002;Dbxref=SEED:fig|1126.22.peg.2,NCBI_GP:AOC50757.1;Name=AOC50757.1;gbkey=CDS;locus_tag=amyaer_0002;product=hypothetical protein;protein_id=AOC50757.1;transl_table=11

 

위 파일로부터 locus tag만 먼저 추출하고 product 를 추출해서 두 정보를 합치려고 한다

locus tag은 있고 product는 없는 경우도 있는데 CDS라고 된 행만을 먼저 추출하여 사용하면 문제 없을 것 같다

 

CDS를 포함하는 행만을 먼저 추출해서 새로운 파일로 저장했다

grep "CDS" LEGE_00239.gff > LEGE_00239.gff.CDS

head -5 LEGE_00236.gff.CDS

JADEVY010000100.1       GeneMarkS-2+    CDS     194     472     .       -       0       ID=cds-MBE9245916.1;Parent=gene-IQ223_15790;Dbxref=NCBI_GP:MBE9245916.1;Name=MBE9245916.1;gbkey=CDS;inference=COORDINATES: ab initio prediction:GeneMarkS-2+;locus_tag=IQ223_15790;product=hypothetical protein;protein_id=MBE9245916.1;transl_table=11
JADEVY010000100.1       Protein Homology        CDS     850     1236    .       -       0       ID=cds-MBE9245917.1;Parent=gene-IQ223_15795;Dbxref=NCBI_GP:MBE9245917.1;Name=MBE9245917.1;gbkey=CDS;inference=COORDINATES: similar to AA sequence:RefSeq:WP_002794531.1;locus_tag=IQ223_15795;product=hypothetical protein;protein_id=MBE9245917.1;transl_table=11
JADEVY010000100.1       Protein Homology        CDS     1463    1693    .       -       0       ID=cds-MBE9245918.1;Parent=gene-IQ223_15800;Dbxref=NCBI_GP:MBE9245918.1;Name=MBE9245918.1;gbkey=CDS;inference=COORDINATES: similar to AA sequence:RefSeq:WP_008205545.1;locus_tag=IQ223_15800;product=type II toxin-antitoxin system CcdA family antitoxin;protein_id=MBE9245918.1;transl_table=11
JADEVY010000100.1       Protein Homology        CDS     2027    3838    .       +       0       ID=cds-MBE9245919.1;Parent=gene-IQ223_15805;Dbxref=NCBI_GP:MBE9245919.1;Name=MBE9245919.1;gbkey=CDS;gene=lepA;inference=COORDINATES: similar to AA sequence:RefSeq:WP_008277316.1;locus_tag=IQ223_15805;product=elongation factor 4;protein_id=MBE9245919.1;transl_table=11
JADEVY010000100.1       Protein Homology        CDS     3894    5072    .       +       0       ID=cds-MBE9245920.1;Parent=gene-IQ223_15810;Dbxref=NCBI_GP:MBE9245920.1;Name=MBE9245920.1;gbkey=CDS;inference=COORDINATES: similar to AA sequence:RefSeq:WP_004158888.1;locus_tag=IQ223_15810;product=AAA family ATPase;protein_id=MBE9245920.1;transl_table=11

위 파일에서는 IQ223_..... 이 부분이 필요하다

 

";locus_tag=" 과 ";product=" 두 문자열을 @로 치환한다

 

그리고 어떤 유전자가 internal stop codon 등의 이유로 fragment로 존재할 결우

locus_tag 뒤에 "partial=true" 라는  표시가 붙는데 이 문자열도 @로 치환한다

 

그리고 2번째 column을 출력하면 locus tag만을 얻을 수 있다

locus tag만 추출하여 LEGE_00239.gff.CDS.locustag 에 저장한다

 

sed -n -e 's/;locus_tag=/@/g' -e 's/;product=/@/g' -e 's/;partial=true/@/g' -e '/$/p' LEGE_00239.gff.CDS | cut -d '@' -f 2 >LEGE_00239.gff.CDS.locustag

 

비슷한 요령으로 product도 추출하여 LEGE_00239.gff.CDS.product 에 저장한다

 

sed -n -e 's/;product=/@/g' -e 's/;protein_id=/@/g' -e 's/;pseudo=true/@/g' -e '/$/p' LEGE_00239.gff.CDS | cut -d '@' -f 2 >LEGE_00239.gff.CDS.product

 

각 파일은 아래와 같이 생겼다

 

head -5 LEGE_00239.gff.CDS.locustag
IQ223_15790
IQ223_15795
IQ223_15800
IQ223_15805
IQ223_15810

head -5 LEGE_00239.gff.CDS.product
hypothetical protein
hypothetical protein
type II toxin-antitoxin system CcdA family antitoxin
elongation factor 4
AAA family ATPase

 

두 파일을 나란히 붙인다

paste LEGE_00239.gff.CDS.locustag LEGE_00239.gff.CDS.product > LEGE_00239.gff.locustag_product

head LEGE_00239.gff.locustag_product

IQ223_15790     hypothetical protein
IQ223_15795     hypothetical protein
IQ223_15800     type II toxin-antitoxin system CcdA family antitoxin
IQ223_15805     elongation factor 4
IQ223_15810     AAA family ATPase
IQ223_15815     hypothetical protein
IQ223_15820     ABC transporter ATP-binding protein
IQ223_15825     DUF4164 domain-containing protein
IQ223_15830     potassium-transporting ATPase subunit A
IQ223_15835     potassium-transporting ATPase subunit KdpB

 

이제 .locustag_product 파일을 사용하면 된다

 

유사한 방식으로 원하는 정보를 추출하면 된다

 

grep "CDS" LEGE_00239.gff > LEGE_00239.gff.CDS
sed -n -e 's/;locus_tag=/@/g' -e 's/;product=/@/g' -e 's/;partial=true/@/g' -e '/$/p' LEGE_00239.gff.CDS | cut -d '@' -f 2 >LEGE_00239.gff.CDS.locustag
sed -n -e 's/;product=/@/g' -e 's/;protein_id=/@/g' -e 's/;pseudo=true/@/g' -e '/$/p' LEGE_00239.gff.CDS | cut -d '@' -f 2 >LEGE_00239.gff.CDS.product
paste LEGE_00239.gff.CDS.locustag LEGE_00239.gff.CDS.product > LEGE_00239.gff.locustag_product