awk와 sed를 사용하여 muti-fasta 파일의 서열 ID를 일괄적으로 바꾸기

수십 개의 대장균 유래 MiSeq 데이터를 Unicycler(GitHub 링크)로 조립하였다. 실행 결과물은 전부 샘플별 디렉토리에 나뉘어 저장된 상태이다. 디렉토리 이름은 전부 다르지만, 그 내용물의 파일 이름은 전부 같다. 예를 들자면 최종 조립 결과 파일의 이름은 모두 assembly.fasta이다. 이를 한데 모아서 저장하자면 파일 이름에 샘플의 이름이 들어가도록 바꾸어야만 한다.

이 작업은 파일의 내용을 직접 건드리는 것이 아니라서 비교적 쉽다.

1 2 3 4 5

$ mkdir final $ find . -name assembly.fasta | while read f > do > cp $f `echo $f | awk -F/ '{sub(/assembly/, ""); print "final/"$2"_plasmid"$3}'` > done

이상의 명령을 실행하면 ./123/assembly.fasta 파일은 ./final/123_plasmid.fasta로 바뀐다. 세 번째 줄을 참고하면 awk의 sub() 함수를 사용하여 assembly라는 문자열을 제거하고, 파일 경로로부터 샘플의 이름(원래는 디렉토리 이름이었음)을 추출하는 것이 핵심이다.

다음으로는 각 fasta 파일로 들어가서 서열 ID를 바꾸는 방법을 알아보자.

1 2 3 4 5 6 7

$ ls *fasta 1039_plasmid.fasta 1536_plasmid.fasta ... $ ls *.fasta | while read f > do > h=${f%%_*} > sed -i "s/>/>s${h}p/" $f > done

현 디렉토리에 존재하는 .fasta 파일의 이름을 변수 f에 넣은 다음 밑줄('_')부터 끝까지를 제거하는 명령(다섯 번째 줄)을 먼저 눈여겨 보자. Shell에서 문자열을 조작하는 매우 유용한 방법이니 잘 기억해 두어야 한다. %%, %, ##, #에 따라서 삭제의 범위와 위치 기준이 달라지고, 치환도 가능하다. 상세한 사항은 Advances Bash-scripting Guide: Manipulating Variables를 참조하자. 다음에는 서열 ID를 sed 명령으로 치환한다. sed의 -i 옵션은 파일을 직접 고치는 것을 뜻한다. 따라서 실수를 하면 곤란하다. 그러나 -i.bak라고 입력하면 원본파일.bak이라는 백업본이 남는다.  흥미롭게도 -i와 SUFFIX 문자열 사이에는 공백이 없어야 한다. Shell variable을 sed의 치환 명령어 내에서 그대로 쓸 수 있다는 것이 매우 편리하다. 단, 명령어를 겹 따옴표로 둘러싸야 한다. 이상을 실행하면 서열의 ID는 다음과 같이 바뀐다.

(변경 전) >1 length=25428 depth=1.01x
(변경 후) >s1039p1 length=25428 depth=1.01x

awk를 사용하면 똑같은 일을 좀 더 복잡하게(?)할 수 있다. Shell variable을 awk 명령어 내에서 쓰려면 -v 옵션을 이용해야 하고, 서열 ID가 있는 라인과 서열만 포함한 라인을 모두 출력하려면 next 명령을 써야 한다. 그러나 printf() 함수를 이용하여 좀 더 세련된 결과를 얻을 수 있다. 다음의 사례에서 세번째 awk 명령어에 주목해 보자. awk의 sub( ) 함수 용법을 먼저 이해해야 하고, next의 의미는 곧 이어서 설명하겠다.

1 2 3 4

$ ls *fasta | while read f > do > awk -v s="${f%%_*}" '$1~/^>/{sub (/>/, ""); printf ">s%d-p%d %s %s\n", s, $1, $2, $3; next}{print}' ${f} > ./edit/${f}; done > done

세 번째 줄의 awk 명령어 구조가 좀 난해하다. 이를 이해하려면 다음의 두 명령어가 무엇이 다른지를 알아야 한다.

1. awk '조건{명령어A}{명령어B}' 파일
2. awk '조건{명령어A; next}{명령어B}' 파일

작업 중이 라인이 조건에 맞으면 {명령어A} 블록을 실행하고, 맞지 않으면 {명령어B} 블록을 실행하는 것이 아니다! 첫 번째 명령어는 조건에 맞으면 A를 실행하고, 같은 라인에 대해서 그대로 B를 수행함을 뜻한다. B를 수행함에 있어서는 앞에서 선언한 조건과는 상관이 없다. 따라서 조건에 맞는 줄에 대해서 A를 실행하고, 다음 줄로 넘어가려면 그 블록 안에서 next 명령을 선언해야 한다. 이렇게 해야 조건에 맞지 않는 줄에 한하여 명령어 B 블록이 수행된다. 이렇게 표현하는 것이 마음에 들지 않으면 if ~ else ~ 구문을 써야 한다.

아무 fasta 파일에 대해서 다음의 명령을 실행시켜 보라. NR은 현재 작업 중인 라인의 번호를 의미한다. 그 다음에는 next를 제거한 뒤 다시 똑같은 명령을 실행해 보라. 어떤 차이가 있는지를 쉽게 이해할 수 있을 것이다.

$ awk '$1~/^>/{print NR, ": header"; next}{print NR, ": seq"}' test.fasta

오늘 소개한 방법이 최선이 아닐 수도 있다. 다음 웹페이지에서 소개하는 짤막한 awk 코드에도 공부할 것이 많이 숨어 있다.

one liner to split a multifasta into separate single file

1 2 3 4

print 'hello world!'cat hg18.fa | awk '{ if (substr($0, 1, 1)==">") {filename=(substr($0,2) ".fa")} print $0 > filename }'