在规则中的函数中使用通配符

如何解决在规则中的函数中使用通配符

我想对管道中的每个样本进行特定的修剪。

我在下面的配置中有一个示例列表：

#######################
#      Config data    # 
#######################

samples:
  3373-1_CCGCGGTT-CTAGCGCT-AHV5HLDSXY_L004:
  3373-2_TTATAACC-TCGATATC-AHV5HLDSXY_L004:
  3373-3_GGACTTGG-CGTCTGCG-AHV5HLDSXY_L004:
  3373-4_AAGTCCAA-TACTCATA-AHV5HLDSXY_L004: 
  ....

我的配置中还有另一个列表，其中包含样本和引物作为 trim5 和 trim3 的字典。

trim5:
  3373-1_CCGCGGTT-CTAGCGCT-AHV5HLDSXY_L004: GATCGGAAGAGCACACGTCTGAACTCCAGTCACCCGCGGTTATCTCGTATGCCGTCTTCTGCTTG 
  3373-2_TTATAACC-TCGATATC-AHV5HLDSXY_L004: GATCGGAAGAGCACACGTCTGAACTCCAGTCACTTATAACCATCTCGTATGCCGTCTTCTGCTTG
  3373-3_GGACTTGG-CGTCTGCG-AHV5HLDSXY_L004: GATCGGAAGAGCACACGTCTGAACTCCAGTCACGGACTTGGATCTCGTATGCCGTCTTCTGCTTG
  3373-4_AAGTCCAA-TACTCATA-AHV5HLDSXY_L004: GATCGGAAGAGCACACGTCTGAACTCCAGTCACAAGTCCAAATCTCGTATGCCGTCTTCTGCTTG
 .....

所以我做了一个函数，允许选择链接到这个样本的引物。

# get_index_trim5/3 allow to verify that we use the correct adaptater depending of the correct sequence
def get_index_trim5():
      for trim3 in config['trim5']:
             if {wildcards.sample} == trim5:
                print (config['trim5'])
                print (config['trim5'][trim5])
                    return ( config['trim5'][trim5] )
             else: 
                    continue 
         
            

# get_index_trim5/3 allow to verify that we use the correct adaptater depending of the correct sequence
def get_index_trim3():
      for trim3 in config['trim3']:      
            if {wildcards.sample} == trim3:
                print (config['trim3'])
                print (config['trim3'][trim3])
                    return ( config['trim3'][trim3] )
             else: 
                    continue
            

  

rule cutadapt_remove_adaptater_trimm:
    priority: 0
    input:
        reads=["../resources/sequences/{sample}_R1.fastq.gz","../resources/sequences/{sample}_R2.fastq.gz"]
    output:
        fastq1= "../results/trimmed/{sample}_R1_trimmed.fastq.gz",fastq2= "../results/trimmed/{sample}_R2_trimmed.fastq.gz",qc= "../results/trimmed/{sample}_qc.txt"
    params:
        adapters="-g %s -a %s -a %s -a %s  -G %s -A %s -A %s -A %s -A %s  "%(get_index_trim5(),config['adaptaters']['adap_R1'],config['adaptaters']['PolyAG'],get_index_trim3(),get_index_trim5(),config['adaptaters']['adap_R2'],config['adaptaters']['PolyG']),extra="--minimum-length 100 -q 20"
    log:
        "../results/logs/trimmed/{sample}_trimmed.log"
    benchmark : 
        "../results/benchmarks/trimmed/{sample}_trimmed_benchmark.txt"
    message:
        """
        --- Trimming on {wildcards.sample} {params.adapters} in process ---
        """
    threads: 4
    resources:
        mem_mb=25000
    shell:
      "cutadapt {params.adapters} {params.extra} -o {output.fastq1} -p {output.fastq2} -j {threads} {input.reads} > {output.qc} 2> {log}"

问题是我需要在我的函数中检查通配符是否对应于引物，但我无法将通配符.sample 作为参数传递给我的函数。 我在使用通配符时遇到了一些困难

你能帮我吗？

提前致谢，祝您有美好的一天， 又名

######################################### ##############

            EDIT Solution

######################################### ###############

你好

我终于找到了解决方案。 第一个问题是我使用了函数 get_index_trim5() 但我们不需要 ()。

然后第二个问题是我试图用 %s 调用一个函数，它返回了一个函数的对象，而不是我想要的键。

所以，它并不完美，但它解决了我的问题，我决定在我的 params 部分将我的适配器变量拆分为不同的变量：R1 和 R2 用于我的经典适配器，没有函数 get_index_trim5 和 { {1}} 和另外两个变量直接使用我的函数。

就像我打电话给我的参数一样。在我的 cutadapt bash 命令中。

get_index_trim3

阿卡

解决方法

我无法将通配符.sample 作为参数传递给我的函数。

您可以通过将自定义函数与 lambda 表达式结合来解决此问题：

def get_index_trim3(sample,config):
    for trim3 in config['trim3']:      
        if sample == trim3:
             print (config['trim3'])
             print (config['trim3'][trim3])
                 return ( config['trim3'][trim3] )
         else: 
             continue

rule cutadapt:
    input:
        lambda wildcards: get_index_trim3(wildcards.sample,config)
    ...

顺便说一下，我不知道您的具体情况，但通常您不需要指定完整的适配器序列。使用第一个公共基础作为所有这些基础的统称，例如：

cutadapt -a GATCGGAAGAGCAC ...

也许试试看它是否有效...