53GB WGS 데이터를 처리하는 29단계 파이프라인 만들기

이어드림스쿨 기업 연계 프로젝트로 반려동물 유전자검사 기업에서 데이터를 받았다. 반려견 한 마리의 전장유전체(WGS) 원시 데이터였는데 압축 파일만 53GB, 시퀀싱 리드로는 약 8억 개였다. 여기서 혈액질환과 관련된 후보 변이를 찾아내는 게 과제였다.

바이오인포매틱스는 이번이 처음이었다. 그래서 1주차에는 품질 확인부터 변이 주석까지 아홉 단계를 전부 손으로 실행했다. 명령 하나를 돌리고 끝나면 결과를 확인한 뒤 다음 명령을 입력하는 식이었다. 분석이 어떻게 흘러가는지 이해하는 데는 이만한 방법이 없었지만 20시간 넘게 걸리는 작업을 계속 이렇게 지켜보고 있을 수는 없었다.

손으로 시작했던 파이프라인에는 나중에 29개 단계가 등록됐다. 목표 단계만 고르면 필요한 선행 작업을 알아서 계산해 실행한다. 아래 영상은 새 시료를 올리고 분석을 시작하는 과정이다.

PawCode 파이프라인 29단계 화면

실수 한 번의 비용

각 단계를 돌려보니 예상보다 훨씬 오래 걸렸다.

단계	실제 실행 시간
정렬	2시간 13분
중복 마킹	1시간 28분
전장 변이 검출	20시간 36분

이 정도 시간이 걸리면 실수 한 번이 비싸진다. 앞 단계 설정이 틀리면 뒤에 이어지는 작업을 전부 다시 돌려야 하는데, 설정이 잘못됐다는 사실은 대개 작업이 다 끝난 뒤에야 드러났다.

한번은 변이 파일을 새로 만들면서 색인 하나만 갱신하고 다른 형식의 낡은 색인을 그대로 남겨둔 적이 있다. 조회 도구는 하필 낡은 쪽을 먼저 읽었고 오류 한 줄 없이 빈 결과만 돌아왔다. 구조 변이 검출 때는 1시간 27분을 돌고 나서 변환 단계에서 실패했다. 당시에는 변환 작업을 별도 단계로 등록해두지 않아서 산출물이 없다는 사실조차 바로 알지 못했다.

이런 일을 몇 번 겪고 나니 서버에 들어가 파일을 뒤지는 식으로는 오래 못 버티겠다 싶었다. 나중에 실행 구조를 잡으면서 자원 관리, 중단 후 재개, 산출물 출처 기록을 유난히 챙긴 것도 이때 덴 기억 때문이다.

파이프라인이 도는 동안에는 현재 단계와 완료된 단계, 남은 시간을 한 화면에서 보게 했다. 실행 기록이 없어서 예상 시간을 계산할 수 없는 단계는 숫자를 지어내는 대신 확인되지 않음으로 표시했다.

실행 중인 파이프라인의 단계별 상태

도구 선택 기준

일곱 가지 정렬 도구를 같은 데이터와 같은 서버 조건에서 비교했다. 처음엔 실행 시간과 정렬률만 보면 되겠지 싶었는데 막상 돌려보니 몇몇 도구는 정렬은커녕 인덱스를 만드는 단계에서 메모리가 부족해 멈췄다. 정렬에 필요한 메모리와 인덱스 생성에 필요한 메모리가 다르다는 걸 이때 알았다.

속도와 정렬 결과, 메모리 사용량을 같이 놓고 minibwa를 골랐다. 전체 데이터를 정렬한 결과 정렬률은 99.80%, 중복 마킹 이후 평균 깊이는 41x였다. 품질 확인은 요약값만 따로 만들지 않고 fastp가 생성한 원본 리포트를 그대로 볼 수 있게 열어뒀다.

필터 통계와 품질 곡선을 확인하는 fastp QC 리포트

변이 검출기는 GATK HaplotypeCaller와 DeepVariant 중에 정하기가 어려워서 1번 염색체만 잘라 둘 다 돌려봤다.

두 도구가 모두 찾은 변이   296,468건 (87.1%)
GATK만 찾은 변이            23,196건
DeepVariant만 찾은 변이     20,864건

1번 염색체에서 두 도구의 결과는 크게 다르지 않았다. GATK HaplotypeCaller는 후보 구간을 국소 조립한 뒤 확률 모델로 유전형을 계산하고, DeepVariant는 정렬된 리드를 이미지 형태로 바꿔 딥러닝 모델로 유전형을 판단한다.

DeepVariant가 딥러닝 모델을 사용한다는 점도 선택에 영향을 줬다. 여기에 DeepVariant 논문의 SNP·소규모 삽입·결실 벤치마크와 멘토 피드백을 참고해 전장 분석의 기본값을 DeepVariant로 정했다. 다만 이번 데이터에는 정답셋이 없어 GATK와 DeepVariant 중 어느 쪽이 더 정확한지는 직접 확인하지 못했다. GATK는 결과를 대조할 때 선택할 수 있도록 남겼고 이 한계도 결과에 같이 적어뒀다.

최종 실행에서 DeepVariant가 검출한 전장 변이는 7,661,621건이었다.

서버 한 대를 나눠 쓰기

프로젝트를 진행하며 제공받은 컴퓨팅 환경은 8코어 서버 한 대였다. 이 서버에서 여러 작업이 같이 돌았기 때문에 무거운 단계 하나가 CPU를 다 쓰면 API 응답이 밀리고 다른 분석도 느려졌다.

각 단계에는 설정한 스레드 수가 아니라 실제로 동시에 사용하는 코어 수를 계산해 등록했다. 정렬처럼 여러 도구가 파이프로 같이 도는 단계는 입력한 스레드 수보다 실제 사용량이 크기 때문이다.

전체 코어에서 API와 DB가 사용할 2코어를 먼저 빼고 나머지를 분석 작업의 예산으로 썼다. 무거운 단계만 요청 순서대로 실행했고 장부에 등록된 사용량과 실제 CPU 사용량 중 큰 값을 기준으로 자리가 있는지 판단했다. 예산이 부족하면 병렬도를 낮췄고 무거운 분석 작업의 CPU 우선순위를 낮춰 API 요청이 들어왔을 때 먼저 처리되도록 했다.

다만 코어 예산을 코드에 박아두지는 않았다. 서버의 코어 수에 따라 예산을 다시 계산하므로 더 큰 서버에서는 단계별 병렬도를 높이거나 여러 작업을 함께 돌릴 수 있다.

중단 후 재개

각 단계가 끝나면 완료 여부와 입력 파일의 지문을 남겼다. 입력 파일과 파라미터가 같고 산출물이 정상이면 건너뛴다. 입력이나 설정이 달라졌을 때만 그 단계부터 다시 돈다.

한 시료에 실행 이력이 27회 쌓였는데 전부 처음부터 돈 건 아니다. 설정을 고치거나 실패한 작업을 다시 시작할 때 이미 끝난 단계는 그대로 썼다. 실행 이력에서는 성공과 실패만 나열하지 않고 어떤 입력과 설정으로 어디까지 갔는지 다시 확인할 수 있게 했다.

27회 누적된 실행 이력

산출물의 출처 기록

참조 자료, 도구 버전, 파라미터, 입출력 건수, 실행 환경까지 11개 항목을 기록했다. 리포트는 숫자를 따로 만들지 않고 이 기록을 읽어 화면과 문서를 생성한다. 다시 실행하면 이전 리포트를 덮어쓰지 않고 새 버전으로 생성된다. 결과가 달라졌을 때 이전 문서와 새 문서를 나란히 놓고 어느 실행에서 나온 값인지 따져보고 싶었다.

실행할 때마다 새 버전으로 생성된 리포트

산출물의 상태는 다섯 가지로 구분했다.

이번 실행에서 만든 것
이전 실행에서 변환한 것
이전 산출물을 그대로 사용한 것
만들지 못한 것
실행하지 않은 것

자주 헷갈린 건 마지막 두 상태였다. 실패한 작업을 0건으로 적으면 정말 결과가 없었던 것처럼 보인다. 그렇다고 실행하지 않은 작업을 실패로 적으면 그것대로 사실과 다르다. 구조 변이 변환 실패처럼 만드는 단계와 읽는 단계가 끊어진 일을 겪은 뒤로 이 구분을 실행 기록과 화면에 같이 적용했다.

리드 수, 변이 수, 행 수

어느 단계에서 무엇이 어떤 기준으로 제외됐는지 파이프라인 화면 한곳에 모았다.

단계별 필터링 현황

이 화면을 만들면서 은근히 애먹은 게 단위였다. 한 변이가 여러 전사체에 걸리면 후보표에는 여러 행으로 나타나서 변이 수와 행 수가 같지 않다. 리드 수, 변이 수, 전사체별 행 수처럼 단위가 다른 숫자를 깔때기 하나로 이어 그리면 실제보다 매끈한 흐름처럼 보인다. 그래서 화면에서는 각 단계가 무엇을 세는지부터 구분했다.

합계도 뺄셈으로 맞추지 않았다. 근거에 걸린 후보와 걸리지 않은 후보를 더하면 18,902건인데 고유 변이는 18,853건이다. 같은 변이라도 어떤 행은 근거에 걸리고 다른 행은 안 걸려서 49건이 양쪽에 들어갔기 때문이다. 화면은 계산으로 맞춘 숫자 대신 리포트가 기록한 합집합을 그대로 쓴다.

PawCode에서는 직접 검증한 실행 구성 하나를 기본값으로 제공한다. 목적별 프로필을 고르거나 대화로 단계를 조합하는 기능까지는 만들지 못했다. 한 번 돌리고 끝나던 명령은 다시 시작할 수 있는 작업 단위가 됐다. 각 결과가 어느 실행에서 나왔는지도 따라갈 수 있다.

처음에는 도구를 순서대로 실행하면 되는 줄 알았다. 직접 해보니 단계마다 확인할 조건이 많았고 한 번의 실수가 긴 재실행으로 이어졌다. 유전체 분석 결과 하나를 얻기까지 이렇게 많은 작업과 확인이 필요하다는 걸 이번에 알게 됐다.