Snakemake: Как мне использовать функцию, которая принимает подстановочный знак и возвращает значение?

Question

У меня есть файлы cram(bam), которые я хочу разделить по группе чтения. Это требует чтения заголовка и извлечения идентификаторов группы чтения.

У меня есть эта функция, которая делает это в моем файле Snakemake:

def identify_read_groups(cram_file):
    import subprocess
    command = 'samtools view -H ' + cram_file + ' | grep ^@RG | cut -f2 | cut -f2 -d":" '
    read_groups = subprocess.check_output(command, shell=True)
    read_groups = read_groups.split('\n')[:-1]
    return(read_groups)

У меня есть все это правило:

rule all:
input:
    expand('cram/RG_bams/{sample}.RG{read_groups}.bam', read_groups=identify_read_groups('cram/{sample}.bam.cram'))

И это правило на самом деле сделать раскол:

rule split_cram_by_rg:
input:
    cram_file='cram/{sample}.bam.cram',
    read_groups=identify_read_groups('cram/{sample}.bam.cram')
output:
    'cram/RG_bams/{sample}.RG{read_groups}.bam'
run:
    import subprocess
    read_groups = open(input.readGroupIDs).readlines()
    read_groups = [str(rg.replace('\n','')) for rg in read_groups]
    for rg in read_groups:
        command = 'samtools view -b -r ' + str(rg) + ' ' + str(input.cram_file) + ' > ' + str(output)
        subprocess.check_output(command, shell=True)

Я получаю эту ошибку при выполнении пробного запуска:

[E::hts_open_format] fail to open file 'cram/{sample}.bam.cram'
samtools view: failed to open "cram/{sample}.bam.cram" for reading: No such file or directory
TypeError in line 19 of /gpfs/gsfs5/users/mcgaugheyd/projects/nei/mcgaughey/EGA_EGAD00001002656/Snakefile:
a bytes-like object is required, not 'str'
File "/gpfs/gsfs5/users/mcgaugheyd/projects/nei/mcgaughey/EGA_EGAD00001002656/Snakefile", line 37, in <module>
  File "/gpfs/gsfs5/users/mcgaugheyd/projects/nei/mcgaughey/EGA_EGAD00001002656/Snakefile", line 19, in identify_read_groups

{sample} не передается в функцию.

Как мне решить эту проблему? Я открыт для других подходов, если я не делаю это "змеиным путем".

==============

РЕДАКТИРОВАТЬ 1

Хорошо, у первого набора примеров, которые я дал, было много проблем.

Вот лучший (?) Набор кода, который, я надеюсь, демонстрирует мою проблему.

import sys
from os.path import join

shell.prefix("set -eo pipefail; ")

def identify_read_groups(wildcards):
    import subprocess
    cram_file = 'cram/' + wildcards + '.bam.cram'
    command = 'samtools view -H ' + cram_file + ' | grep ^@RG | cut -f2 | cut -f2 -d":" '
    read_groups = subprocess.check_output(command, shell=True)
    read_groups = read_groups.decode().split('\n')[:-1]
    return(read_groups)

SAMPLES, = glob_wildcards(join('cram/', '{sample}.bam.cram'))
RG_dict = {}
for i in SAMPLES:
    RG_dict[i] = identify_read_groups(i)

rule all:
    input:
        expand('{sample}.boo.txt', sample=list(RG_dict.keys()))

rule split_cram_by_rg:
    input:
        file='cram/{sample}.bam.cram',
        RG = lambda wildcards: RG_dict[wildcards.sample]
    output:
        expand('cram/RG_bams/{{sample}}.RG{input_RG}.bam') # I have a problem HERE. How can I get my read groups values applied here? I need to go from one cram to multiple bam files split by RG (see -r in samtools view below). It can't pull the RG from the input.
    shell:
        'samtools view -b -r {input.RG} {input.file} > {output}'


rule merge_RG_bams_into_one_bam:
    input:
        rules.split_cram_by_rg.output
    output:
        '{sample}.boo.txt'
    message:
        'echo {input}'
    shell:
        'samtools merge {input} > {output}' #not working
        """

==============

РЕДАКТИРОВАТЬ 2

Становится НАМНОГО ближе, но в настоящее время борется с расширением, правильно создавая файлы Lane Bam и сохраняя подстановочные знаки

Я использую этот цикл для создания промежуточных имен файлов:

for sample in SAMPLES:
    for rg_id in list(return_ID(sample)):
        out_rg_bam.append("temp/lane_bam/{}.ID{}.bam".format(sample, rg_id))

return_ID это функция, которая берет шаблон подстановки и возвращает список групп чтения, которые содержит образец

Если я использую out_rg_bam в качестве входных данных для правила слияния, тогда ВСЕ файлы объединяются в объединенный элемент, а не разделяются sample,

Если я использую expand('temp/realigned/{{sample}}.ID{rg_id}.realigned.bam', sample=SAMPLES, rg_id = return_ID(sample)) тогда rg_id применяется к каждому образцу. Поэтому, если у меня есть два образца (a,b) с группами чтения (0,1) и (0,1,2), я получу a0, a1, a0, a1, a2 и b0, b1, b0, b1, Би 2.

Answer 1

Я собираюсь дать более общий ответ, чтобы помочь другим, которые могут найти эту тему. Snakemake применяет подстановочные знаки только к строкам в разделах "input" и "output", когда строки перечислены непосредственно, например:

input:
    '{sample}.bam'

Если вы пытаетесь использовать функции, как вы были здесь:

input:
    read_groups=identify_read_groups('cram/{sample}.bam.cram')

Замена подстановочного знака не будет выполнена. Вы можете использовать лямбда-функцию и сделать замену самостоятельно:

input:
    read_groups=lambda wildcards: identify_read_groups('cram/{sample}.bam.cram'.format(sample=wildcards.sample))

Answer 2

Попробуйте это: я использую id = 0, 1, 2, 3 для именования выходного файла BAM в зависимости от того, сколько readgroup для файла BAM.

## this is a regular function which takes the cram file, and get the read-group to 
## construct your rule all
## you actually just need the number of @RG, below can be simplified  
def get_read_groups(sample):
    import subprocess
    cram_file = 'cram/' + sample + '.bam.cram'
    command = 'samtools view -H ' + cram_file + ' | grep ^@RG | cut -f2 | cut -f2 -d":" '
    read_groups = subprocess.check_output(command, shell=True)
    read_groups = read_groups.decode().split('\n')[:-1]
    return(read_groups)

SAMPLES, = glob_wildcards(join('cram/', '{sample}.bam.cram'))
RG_dict = {}
for sample in SAMPLES:
    RG_dict[sample] = get_read_groups(sample)

outbam = []
for sample in SAMPLES:
    read_groups = RG_dict[sample]
    for i in range(len(read_groups)):
        outbam.append("{}.RG{}.bam".format(sample, id))


rule all:
    input:
        outbam

## this is the input function, only uses wildcards as argument 
def identify_read_groups(wildcards):
    import subprocess
    cram_file = 'cram/' + wildcards.sample + '.bam.cram'
    command = 'samtools view -H ' + cram_file + ' | grep ^@RG | cut -f2 | cut -f2 -d":" '
    read_groups = subprocess.check_output(command, shell=True)
    read_groups = read_groups.decode().split('\n')[:-1]
    return(read_groups[wildcards.id])

rule split_cram_by_rg:
    input:
        cram_file='cram/{sample}.bam.cram',
        read_groups=identify_read_groups
    output:
        'cram/RG_bams/{sample}.RG{id}.bam'  
    run:
        import subprocess
        read_groups = input.read_groups
        for rg in read_groups:
            command = 'samtools view -b -r ' + str(rg) + ' ' + str(input.cram_file) + ' > ' + str(output)
            subprocess.check_output(command, shell=True)

когда используете змеиный мейк, продумывайте путь снизу вверх. Сначала определите, что вы хотите сгенерировать в правиле all, а затем создайте правило, чтобы создать окончательное all.

Answer 3

Ваше правило не может иметь подстановочные знаки. Это не подстановочная зона.

РЕДАКТИРОВАТЬ 1

Я набрал этот псевдокод в Notepad++, он не предназначен для компиляции, просто пытаясь создать фреймворк. Я думаю, что это больше, чем вы после.

Используйте функцию внутри раскрытия, чтобы сгенерировать список имен файлов, которые затем будут использоваться для управления всем правилом конвейера Snakemake. Переменные baseSuffix и basePrefix просто для того, чтобы дать вам представление о передаче строки, здесь разрешены аргументы. При возврате списка строк вам нужно будет их распаковать, чтобы Snakemake правильно прочитал результат.

def getSampleFileList(String basePrefix, String baseSuffix){
    myFileList = []
    ListOfSamples = *The wildcard glob call*
    for sample in ListOfSamples:
        command = "samtools -h " + sample + "SAME CALL USED TO GENERATE LIST OF HEADERS"
        for rg in command:
            myFileList.append(basePrefix + sample + ".RG" + rg + baseSuffix)
}


basePreix = "cram/RG_bams/"
baseSuffix = ".bam" 

rule all:
    input:
        unpack(expand("{fileName}", fileName=getSampleFileList(basePrefix, baseSuffix)))


rule processing_rg_files:
    input:
        'cram/RG_bams/{sample}.RG{read_groups}.bam'
    output:
        'cram/RG_TXTs/{sample}.RG{read_groups}.txt'
    run:
        "Let's pretend this is useful code"

Конец редактирования

Если бы это не было во всем правиле, вы бы использовали встроенные функции

Так что я не уверен, чего вы пытаетесь достичь. Согласно моим догадкам, прочитайте ниже некоторые заметки о вашем коде.

rule all:
input:
    expand('cram/RG_bams/{sample}.RG{read_groups}.bam', read_groups=identify_read_groups('cram/{sample}.bam.cram'))

Пробный запуск завершается неудачно, когда он вызывает функцию "identif_read_groups" внутри правила all call. Он передается в вызов вашей функции в виде строки, а не подстановочного знака.

Технически, если вызов samtools не завершился неудачно, а вызов функции identifi_read_groups (cram_file) вернул список из 5 строк, он расширился бы примерно до следующего:

rule all:
    input:
        'cram/RG_bams/{sample}.RG<output1FromFunctionCall>.bam',
        'cram/RG_bams/{sample}.RG<output2FromFunctionCall>.bam',
        'cram/RG_bams/{sample}.RG<output3FromFunctionCall>.bam',
        'cram/RG_bams/{sample}.RG<output4FromFunctionCall>.bam',
        'cram/RG_bams/{sample}.RG<output5FromFunctionCall>.bam'

Но термин "{sample}" на этом этапе предварительной обработки Snakemake считается строкой. Как вам нужно было обозначать подстановочные знаки в функции расширения с {{}}.

Посмотрите, как я обращаюсь к каждой переменной Snakemake, которую я объявляю для своего правила для всех входных вызовов и не использую подстановочные знаки:

expand("{outputDIR}/{pathGVCFT}/tables/{samples}.{vcfProgram}.{form[1][varType]}{form[1][annotated]}.txt", outputDIR=config["outputDIR"], pathGVCFT=config["vcfGenUtil_varScanDIR"], samples=config["sample"], vcfProgram=config["vcfProgram"], form=read_table(StringIO(config["sampleFORM"]), " ").iterrows())

В этом случае read_table возвращает 2-мерный массив в форму. Snakemake хорошо поддерживается Python. Мне это нужно для сопряжения разных аннотаций с разными вариантами.

Все ваше правило должно быть строкой или списком строк в качестве входных данных. Вы не можете иметь подстановочные знаки в своем правиле "все". Эти правила используют все входные строки, которые Snakemake использует для генерации совпадений с другими символами подстановки. Создайте все имя файла в вызове функции и верните его, если вам нужно.

Я думаю, что вы должны просто превратить это в нечто вроде этого:

rule all:
input:
    expand("{fileName}", fileName=myFunctionCall(BecauseINeededToPass, ACoupleArgs))

Также рассмотрите возможность обновления, чтобы оно было более общим:

rule split_cram_by_rg:
    input:
        cram_file='cram/{sample}.bam.cram',
        read_groups=identify_read_groups('cram/{sample}.bam.cram')

У него может быть два или более подстановочных знака (почему мы любим Snakemake). Вы можете получить доступ к подстановочным знакам позже в директиве Python "run" через объект подстановочных знаков, поскольку, похоже, вы захотите использовать их для каждого цикла. Я думаю, что подстановочные знаки ввода и вывода должны совпадать, так что, возможно, попробуйте это так же.

rule split_cram_by_rg:
    input:
        'cram/{sample}.bam.cram'
    output:
        expand('cram/RG_bams/{{sample}}.RG{read_groups}.bam', read_groups=myFunctionCall(BecauseINeededToPass, ACoupleArgs))
    ...
    params:
          rg=myFunctionCall(BecauseINeededToPass, ACoupleArgs)
    run:
        command = 'Just an example ' +  + str(params.rg)

Опять же, я не уверен, что вы пытаетесь сделать, я не уверен, что мне нравится идея вызова функции дважды, но, эй, она запустится;P Также обратите внимание на использование подстановочного знака "sample" в директиве input внутри строки {} и в выходной директиве внутри раскрытия {{}}.

Пример доступа к групповым символам в вашей директиве run

Пример вызова функций в тех местах, о которых вы не думаете. Я схватил поля VCF, но это могло быть что угодно. Я использую внешний конфигурационный файл здесь.