snakemake[3:一个简单的WGBS流程]-创新互联

参考DNA甲基化数据分析全流程
搭建了一个简单的管道：
先上代码

samples:
    cd133/a: SAM/cd133/a.sam
    cd133/b: SAM/cd133/b.sam
    cd133/c: SAM/cd133/c.sam
    cd34/a: SAM/cd34/d.sam
    cd34/a: SAM/cd34/e.sam
    cd34/a: SAM/cd34/f.sam
   

configfile:"config.yaml"

rule all:
    input:
        expand("SAM/{sample}/bedGraph",group=config["samples"])

def get_input_sam(wildcards):
    return config["samples"][wildcards.sample]

    
rule samToBam:
    input:
        get_input_sam
    output:
        "SAM/{sample}.bam"
    shell:
        "samtools view -S -b -o {output} {input}"

rule bamToSortBam:
    input:
        "SAM/{sample}.bam"
    output:
        "SAM/{sample}.sort.bam"
    shell:
        "samtools sort -m 2000000000 {input} {output}"

rule duplication:
    input:
        "SAM/{sample}.sort.bam"
    output:
        "SAM/{sample}.dup.bam"
    shell:
        "sambamba markdup -t 5 {input} {output} 2>markdup_report.txt"

rule getBadGraph:
    input:
        "SAM/{sample}.dup.bam"
    output:
        "SAM/{group}/bedGraph"
    shell:
        """
        MethylDackel extract /home/data/GRCh37.primary_assembly.genome.fa {input} -o {output}
        MethylDackel extract /home/data/GRCh37.primary_assembly.genome.fa {input} -o {output} --CHG
        MethylDackel extract /home/data/GRCh37.primary_assembly.genome.fa {input} -o {output} --CHH
        """

我的数据有很多的样本，每个样本里也有细分的数据。
刚开始想通过嵌套，先获得groups里的group，然后利用group的名字建立键值对对下一级目录进行访问的，但搞了一晚上发现好像有点难，于是直接在samples的键前面加上了对应的目录。
如果有大佬看到了怎么嵌套或有更简洁的方法求指路！

.sam—samtools—>.bam
.bam—samtools—>.sort.bam
.sort.bam—sambamba—>.dum.bam
.dum.bam—MethylDackel—>bedGraph

• 在想使用多行shell命令可以使用"""进行标注
• get_input_sam(wildcards)函数的实际意义是返回键值对中的值（路径）
• all或者最后一步中利用expand定义输出文件的通配符内容

看了一篇文章里面提到了利用文件建立数据索引的方法。

首先创建一个csv文件

import pandas as pd
samples_table = pd.read_csv("index.csv").set_index("id", drop=False)
print(samples_table.loc[1,"path"])
print(samples_table.index)
rule all:
    input:
         "results/awesome2"
         
def sam_from_sample(wildcards):
  return samples_table.loc[int(wildcards.sample), "path"]

rule a:
    input: 
        sam_from_sample
    output:
        "results/awesome/{sample}_R1.fq"
    shell:
        "echo "
        
rule b:
    input: 
        expand("results/awesome/{sample}_R1.fq",sample=samples_table.index)
    output:
        "results/awesome2"
    shell:
        "echo "

※利用all才可以串联规则a，b
samples_table.loc[int(wildcards.sample), "path"]返回地址
sample=samples_table.index返回索引号
在b中通过id确定a的输出文件然后自动获得对应地址

你是否还在寻找稳定的海外服务器提供商？创新互联www.cdcxhl.cn海外机房具备T级流量清洗系统配攻击溯源，准确流量调度确保服务器高可用性，企业级服务器适合批量采购，新人活动首月15元起，快前往官网查看详情吧