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
"""
1 路径我的数据有很多的样本,每个样本里也有细分的数据。
刚开始想通过嵌套,先获得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文件
id | name |
---|---|
1 | /home/…/…sam |
2 | /home/…/…sam |
… | … |
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,bsamples_table.loc[int(wildcards.sample), "path"]
返回地址sample=samples_table.index
返回索引号
在b中通过id确定a的输出文件然后自动获得对应地址
你是否还在寻找稳定的海外服务器提供商?创新互联www.cdcxhl.cn海外机房具备T级流量清洗系统配攻击溯源,准确流量调度确保服务器高可用性,企业级服务器适合批量采购,新人活动首月15元起,快前往官网查看详情吧
文章标题:snakemake[3:一个简单的WGBS流程]-创新互联
网页网址:http://myzitong.com/article/dihodd.html