[roary] query_pan_genome 스크립트로 roary 결과 다루기 (1) core gene의 protein id만 구하기

34개의 유전체를 roary로 분석하였다. amino acid identity 50% 기준이다

 

roary --f ./ -e -n -i 50 -p 64 *.gff

분석한 모든 유전체의 core gene을 query_pan_genome 스크립트로 구하였다

그 전에 분석에 사용한 모든 roary gff 파일과

query_pan_genome -a intersection *.gff

결과 파일은 확장자도 없이 'pan_genome_results' 라는 파일로 나왔다.

엑셀에서 열어보면 아래와 같다

자세히 살펴보면

gene name (gene name이 있는 경우에만, 없을 땐 protein id)가 맨 앞에 있고 그 다음에 콜론(:)이 공백(space)로 구분되어 있고 그 뒤로는 탭으로 구분된 protein id가 나열되어 있다.

 

즉, pan_genome_results 파일의 형식은

 

한 줄 (row)는 identity 50% 이상으로 (처음 roary 실행할 때 설정했던 값, default는 아마도 95%였던것 같다) 같은 protein group에 속하는 protein이 나열되어 있고,

첫 번째 컬럼은 그 protein group의 대표이름을 gene name으로 표시한 뒤 공백으로 구분하고,

두 번째 컬럼부터 각 유전체에서 해당 protein group에 속하는 protein id를 탭으로 구분하여 나열해놓은 것으로 보인다.

 

그런데 여기서 말하는 protein id라는게 표현 형식이 이상하다.

roary 분석 시작할 때 fixed input file로 바뀌었던 유전체들은 흔히 보던 NCBI protein accession number 형식이 아니라 이상한 랜덤의 문자열과 (f1607fa491...) 하이픈, 숫자로 바뀌어져있다. (e.g. f1607fa4916f5be75e5b256c462ecccc_4887) 이걸 원래 protein id 또는 locus tag이 무엇인지 찾으려면 어떻게 해야하는지 아직은 모르겠지만 일단 넘어가자.

 

지금 내가 하고 싶은 것은 앞 뒤 필요없는 것은 다 버리고 한 컬럼의 protein id 만 깔끔하게 얻어서 KEGG pathway에 mapping 하고 싶다.

 

(사실 이 작업은 엑셀에서도 할 수 있다. openrefine을 사용하면 좀 더 편리하고 정교하게 할 수 있다. 리눅스 커맨드를 새로 배워서라도 사용하는 것이 당장은 힘들어도 앞으로 활용도가 높을 것 같다.)

 

몰라서 이것저것 찾아보았다.

리눅스의 명령어만으로 이 작업이 가능하다

cat으로 파일을 출력하고 awk로 원하는 컬럼만 다시 출력하고 cut으로 구분자를 하이픈(-)으로 설정하면 될 것 같다.

시행착오 끝에 얻은 커맨드는 아래와 같다. 사용한 유전체가 FACHB-1342 strain이라서 아래와 같은 텍스트 파일로 결과를 저장하였다.

cat pan_genome_results | awk '{print $2}'| cut -d '-' -f 2 > protein_id_core_gene_FACHB-1342.txt

결과물은 아래와 같다

MBD2600753.1
MBD2600724.1
MBD2602888.1
MBD2602306.1
.
.
.

위의 엑셀 파일 캡쳐한 것과 비교하면 어떤 결과를 얻었는지 알 수 있다.

 

다음은 KEGG pathway에 core gene을 mapping 하여 시각적으로 core metabolic pathway를 구성하는 것을 기록해보자

 

roary input file (gff) 준비에 대한 내용은 이전 포스팅에 정리하였다. 스크립트를 이용하여 호로록 준비하는 수준은 아니지만 손으로 하는 작업 치곤 그나마 자동화가 되어있다 (못생긴 애들 중에 제일 잘 생긴 느낌?) NCBI에서 gff와 fna 파일을 따로 다운로드 받아서 roary가 사용할 수 있는 gff 형식으로 합치는 방식이다. 리눅스 안쓰고 탐색기에서 내가 만든 배치파일을 더블클릭하면 된다. 아래 참고.

https://fine-tuned.tistory.com/35

Roary (Pan-genome pipeline) 관련 메모

 

명령어를 잘 몰라서 찾아가면서 했는데 참고한 블로그들은 아래와 같다.

감사한 마음으로 애드블록을 꺼드렸다^^ 앞으로는 광고 한 개씩 클릭해드려야겠다.

이 글을 보시는 분들도 애드블록을 꺼주시면 내 인생에 less than 1원 정도 도움이 된다.

 

-참고-

 

[Linux] 파이프(pipe)에 대한 이해

https://gracefulprograming.tistory.com/92

 

리눅스 awk 명령어 사용법. (Linux awk command) - 리눅스 파일 텍스트 데이터 검사, 조작, 출력.

(awk를 알고 싶으면 이 곳을 먼저 읽어보시길)

https://recipes4dev.tistory.com/171

 

Roary 'query_pan_genome'의 고약한 출력물

(정해영 박사님의 글은 내가 이해하기에 수준이 높지만 그래도 항상 참고하고 있다. 정 박사님 블로그에 애드센스 넣으셔도 될 것 같은데...?)

http://blog.genoglobe.com/2021/11/roary-querypangenome.html