expand를 사용하여 파일 이름 생성하기

Snakemake 와일드카드는 많은 파일에 규칙을 쉽게 적용할 수 있게 해주지만, 동시에 새로운 과제를 안겨줍니다. 바로 여러분이 원하는 모든 파일 이름을 어떻게 생성할 것인가 하는 점입니다.

이러한 과제의 예로, 챕터 8의 compare_genomes 규칙에 필요한 게놈 목록을 생각해 보십시오.

rule compare_genomes:
    input:
        "GCF_000017325.1.fna.gz.sig",
        "GCF_000020225.1.fna.gz.sig",
        "GCF_000021665.1.fna.gz.sig",
        "GCF_008423265.1.fna.gz.sig",

이 목록은 와일드카드 규칙에 의해 생성될 스케치들을 지정하기 때문에 매우 중요합니다. 하지만 파일 이름의 일부가 동일하고 반복되기 때문에 이 목록을 일일이 작성하는 것은 번거롭고 실수가 발생하기 쉽습니다.

더 나쁜 것은, 이 목록을 여러 곳에서 사용해야 하거나 동일한 접근 번호로 약간 다른 파일 이름을 만들어야 하는 경우입니다. 목록의 요소를 추가, 삭제 또는 편집하고 싶을 때 여러 곳을 변경해야 하므로 실수가 생길 가능성이 큽니다.

챕터 9에서는 이를 Snakefile 상단의 접근 번호 리스트로 바꾸고 expand라는 함수를 사용하여 목록을 생성하는 방식을 보여드렸습니다.

ACCESSIONS = ["GCF_000017325.1",
              "GCF_000020225.1",
              "GCF_000021665.1",
              "GCF_008423265.1"]

#...

rule compare_genomes:
    input:
        expand("{acc}.fna.gz.sig", acc=ACCESSIONS),

expand를 사용하여 파일 이름 목록을 생성하는 것은 Snakefile에서 흔히 볼 수 있는 패턴입니다. 이번 장에서 좀 더 자세히 살펴보겠습니다!

단일 패턴과 하나의 값 리스트로 expand 사용하기

위의 예제에서 우리는 {acc}.fna.gz.sig라는 단일 패턴을 제공하고, ACCESSIONS의 각 요소에서 acc라는 필드 이름을 채워 여러 파일 이름으로 풀도록 expand에 요청했습니다. (입력 및 출력 블록에서 값을 지정하기 위한 키워드 구문인 acc=ACCESSIONS를 보셨을 것입니다.)

여기서 expand('{acc}.fna.gz.sig', acc=...)의 결과는 파일 이름 네 개를 긴 형식으로 적는 것과 완전히 동일합니다.

"GCF_000017325.1.fna.gz.sig",
"GCF_000020225.1.fna.gz.sig",
"GCF_000021665.1.fna.gz.sig",
"GCF_008423265.1.fna.gz.sig"

즉, expand는 특별한 와일드카드 매칭이나 패턴 추론을 수행하지 않습니다. 단지 값을 채워 넣고 결과 리스트를 반환할 뿐입니다.

여기서 ACCESSIONS는 리스트, 튜플 또는 딕셔너리와 같은 모든 Python 반복 가능한 객체(iterable)가 될 수 있습니다. 자세한 내용은 Python 부록을 참고하세요.

여러 값 리스트로 expand 사용하기

여러 필드 이름과 함께 expand를 사용할 수도 있습니다. 다음을 고려해 보세요.

expand('{acc}.fna.{extension}', acc=ACCESSIONS, extension=['.gz.sig', '.gz'])

이것은 제공된 패턴에 acc와 extension의 가능한 모든 조합을 대입하여 다음과 같은 8개의 파일 이름을 생성합니다.

"GCF_000017325.1.fna.gz.sig",
"GCF_000017325.1.fna.gz",
"GCF_000020225.1.fna.gz.sig",
"GCF_000020225.1.fna.gz",
"GCF_000021665.1.fna.gz.sig",
"GCF_000021665.1.fna.gz",
"GCF_008423265.1.fna.gz.sig",
"GCF_008423265.1.fna.gz"

모든 조합 생성 vs 쌍별(pairwise) 조합 생성

위에서 보았듯이, 여러 패턴이 있는 경우 expand는 가능한 모든 조합을 생성합니다. 즉:

# 대상: -n all_product all_zip all_zip_short

# ANCHOR: 조합
X = [1, 2, 3]
Y = ['a', 'b', 'c']

rule all_product:
   input:
      expand('{x}.by.{y}', x=X, y=Y)
# ANCHOR_END: 조합

assert len(expand('{x}.by.{y}', x=X, y=Y)) == 9

# ANCHOR: zip
X = [1, 2, 3]
Y = ['a', 'b', 'c']

rule all_zip:
   input:
      expand('{x}.by.{y}', zip, x=X, y=Y)
# ANCHOR_END: zip

assert len(expand('{x}.by.{y}', zip, x=X, y=Y)) == 3

# ANCHOR: zip_short
X = [1, 2, 3]
Y = ['a', 'b']

rule all_zip_short:
   input:
      expand('{x}.by.{y}', zip, x=X, y=Y)
# ANCHOR_END: zip_short

assert len(expand('{x}.by.{y}', zip, x=X, y=Y)) == 2

rule create:
    output:
        touch("{x}.by.{y}")

이는 1.by.a, 1.by.b, 1.by.c, 2.by.a 등 9개의 파일 이름을 생성합니다. 만약 expand 문자열에 세 번째 패턴을 추가한다면, expand는 그것 역시 조합에 추가할 것입니다!

어떻게 이런 일이 일어나는 걸까요?

기본적으로 expand는 가능한 모든 조합을 포함하는 “all-by-all” 확장을 수행합니다. (이를 데카르트 곱(Cartesian product), 교차 곱(cross-product) 또는 외부 조인(outer join)이라고도 부릅니다.)

하지만 항상 그런 방식을 원하는 것은 아닙니다. 이 동작을 어떻게 바꿀 수 있을까요?

expand 함수는 선택적으로 두 번째 인자인 결합자(combinator)를 받는데, 이는 뒤따르는 값 리스트들을 어떻게 결합할지 expand에 알려줍니다. 기본적으로 expand는 가능한 모든 조합을 만드는 Python 함수인 itertools.product를 사용하지만, 다른 함수를 지정할 수 있습니다.

특히, product 대신 zip을 사용하여 쌍별(pairwise) 조합을 생성하도록 expand에 지시할 수 있습니다. 이는 와일드카드 예시 중 하나에서 이미 수행했던 작업입니다.

다음은 그 예시입니다:

X = [1, 2, 3]
Y = ['a', 'b', 'c']

rule all_zip:
   input:
      expand('{x}.by.{y}', zip, x=X, y=Y)

이제 1.by.a, 2.by.b, 3.by.c 단 세 개의 파일 이름만 생성됩니다.

여기서 주의할 점은 zip이 가장 짧은 입력 리스트의 길이에 맞춰 출력 리스트를 만든다는 것입니다. 따라서 요소가 세 개인 리스트 하나와 요소가 두 개인 리스트 하나를 주면, 첫 번째 리스트에서 두 개의 요소만 사용하게 됩니다.

예를 들어, 이 Snakefile의 expand에서:

X = [1, 2, 3]
Y = ['a', 'b']

rule all_zip_short:
   input:
      expand('{x}.by.{y}', zip, x=X, y=Y)

두 번째 리스트에 3에 대응하는 짝이 없기 때문에 1.by.a와 2.by.b만 생성됩니다.

자세한 내용은 product 대신 zip을 사용하는 방법에 관한 snakemake 문서를 참고하세요.

expand에서 사용할 식별자 목록 가져오기

expand 함수는 워크플로에서 여러 번 사용하는 식별자 목록이 있을 때 효과적인 솔루션을 제공합니다. 이는 생물정보학에서 매우 흔한 패턴입니다! 하지만 위의 예시들처럼 Snakefile에 이 목록들을 일일이 적는 것이 항상 실용적인 것은 아닙니다. 수십에서 수백 개의 식별자가 있을 수 있기 때문입니다!

식별자 목록은 다양한 방식으로 다른 파일에서 로드할 수 있으며, glob_wildcards를 사용하여 디렉토리의 실제 파일 세트로부터 생성할 수도 있습니다.

파일 또는 디렉토리에서 접근 번호 목록을 로드하는 예시

텍스트 파일에서 접근 번호 목록 불러오기

accessions.txt라는 텍스트 파일에 다음과 같이 간단한 접근 번호 목록이 있다면:

accessions.txt 파일 내용:

GCF_000017325.1
GCF_000020225.1
GCF_000021665.1
GCF_008423265.1

다음 코드는 텍스트 파일의 각 줄을 별개의 ID로 로드합니다:

with open('accessions.txt', 'rt') as fp:
    ACCESSIONS = fp.readlines()
    ACCESSIONS = [ line.strip() for line in ACCESSIONS ]
    
print(f'ACCESSIONS는 길이가 {len(ACCESSIONS)}인 Python 리스트입니다.')
print(ACCESSIONS)

rule all:
    input:
        expand("{acc}.sig", acc=ACCESSIONS)

rule sketch_genome:
    input:
        "genomes/{accession}.fna.gz",
    output:
        "{accession}.sig",
    shell: """
        sourmash sketch dna -p k=31 {input} --name-from-first -o {output}
    """

그리고 각 접근 번호에 대해 sourmash 서명을 빌드합니다.

CSV 파일의 특정 열에서 불러오기

텍스트 파일 대신 여러 열이 있는 CSV 파일이 있고 로드할 ID들이 모두 한 열에 있다면, Python의 pandas 라이브러리를 사용하여 CSV를 읽어올 수 있습니다. 아래 코드에서 pandas.read_csv는 CSV를 pandas DataFrame 객체로 로드하고, 우리는 accession 열을 선택하여 이를 반복 가능한 객체로 사용합니다.

accessions.csv 파일 내용:

accession,information
GCF_000017325.1,genome 1
GCF_000020225.1,genome 2
GCF_000021665.1,genome 3
GCF_008423265.1,genome 4

accessions.csv를 로드하기 위한 Snakefile:

import pandas

CSV_DATAFRAME = pandas.read_csv('accessions.csv')
ACCESSIONS = CSV_DATAFRAME['accession']

print(f'ACCESSIONS는 길이가 {len(ACCESSIONS)}인 pandas Series입니다.')
print(ACCESSIONS)

rule all:
    input:
        expand("{acc}.sig", acc=ACCESSIONS)

rule sketch_genome:
    input:
        "genomes/{accession}.fna.gz",
    output:
        "{accession}.sig",
    shell: """
        sourmash sketch dna -p k=31 {input} --name-from-first -o {output}
    """

설정 파일에서 불러오기

Snakemake는 설정(configuration) 파일의 사용도 지원합니다. Snakefile은 기본 설정 파일의 이름을 제공하며, 이는 다시 명령줄에서 재정의될 수 있습니다.

설정 파일은 접근 번호들을 넣어두기에 좋은 장소이기도 합니다. 다음을 고려해 보세요:

accessions:
- GCF_000017325.1
- GCF_000020225.1
- GCF_000021665.1
- GCF_008423265.1

이것은 다음 Snakefile에서 사용됩니다:

configfile: "config.yml"

ACCESSIONS = config['accessions']
    
print(f'ACCESSIONS는 길이가 {len(ACCESSIONS)}인 Python 리스트입니다.')
print(ACCESSIONS)

rule all:
    input:
        expand("{acc}.sig", acc=ACCESSIONS)

rule sketch_genome:
    input:
        "genomes/{accession}.fna.gz",
    output:
        "{accession}.sig",
    shell: """
        sourmash sketch dna -p k=31 {input} --name-from-first -o {output}
    """

여기서 config.yml은 사람이 읽을 수 있으면서 컴퓨터도 읽을 수 있는 형식인 YAML 파일입니다. 설정 파일에 대해서는 나중에 다시 다루겠습니다!

glob_wildcards를 사용하여 파일 세트에서 ID 또는 접근 번호 불러오기

우리는 와일드카드 장에서 glob_wildcards 명령어를 잠시 소개했습니다: glob_wildcards는 디렉토리에 실제로 존재하는 파일들에 대해 패턴 매칭을 수행합니다.

다음은 실제 파일 이름으로부터 네 개의 접근 번호를 가져오기 위해 glob_wildcards를 사용하는 Snakefile입니다:

GLOB_RESULTS = glob_wildcards("genomes/{acc}.fna.gz")
ACCESSIONS = GLOB_RESULTS.acc

print(f'ACCESSIONS는 길이가 {len(ACCESSIONS)}인 Python 리스트입니다.')
print(ACCESSIONS)

rule all:
    input:
        expand("{acc}.sig", acc=ACCESSIONS)

rule sketch_genome:
    input:
        "genomes/{accession}.fna.gz",
    output:
        "{accession}.sig",
    shell: """
        sourmash sketch dna -p k=31 {input} --name-from-first -o {output}
    """

이는 접근 번호 목록을 가져오는 매우 편리한 방법이지만, 위험할 수도 있습니다. 특히 실수로 파일을 삭제하고도 샘플 하나가 빠졌다는 사실을 알아차리지 못하기 쉽습니다! 그런 이유로 많은 상황에서 로드할 독립적인 파일 목록을 제공하는 것을 권장합니다.

와일드카드와 expand - 맺음말

와일드카드와 결합된 expand는 매우 강력하고 유용합니다. 하지만 와일드카드와 마찬가지로 이러한 강력함에는 약간의 복잡성이 따릅니다. 다음은 이러한 기능들이 어떻게 결합되는지에 대한 간단한 요약입니다.

• expand 함수는 패턴과 채워 넣을 값 리스트를 사용하여 생성할 파일 목록을 만듭니다.
• 규칙의 와일드카드는 이름이 패턴과 일치하는 파일을 생성하기 위한 레시피를 제공합니다.

일반적으로 Snakefile에서 우리는 특정 패턴과 일치하는 파일 목록을 생성하기 위해 expand를 사용하고, 그런 다음 해당 파일들을 실제로 생성하기 위해 와일드카드를 사용하는 규칙을 작성합니다.

expand와 함께 사용할 값 목록은 텍스트 파일, CSV 파일, 설정 파일을 포함한 다양한 소스에서 가져올 수 있습니다. 또한 패턴을 사용하여 실제로 존재하는 파일에서 값 목록을 추출하는 glob_wildcards에서 가져올 수도 있습니다.

링크 및 참고 자료

• expand 함수에 대한 Snakemake 참고 문서
• Python itertools 문서