参考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的键前面加上了对应的目录。
如果有大佬看到了怎么嵌套或有更简洁的方法求指路!
2 流程
.sam—samtools
—>.bam
.bam—samtools
—>.sort.bam
.sort.bam—sambamba
—>.dum.bam
.dum.bam—MethylDackel
—>bedGraph
3 小结
- 在想使用多行shell命令可以使用
"""
进行标注 - get_input_sam(wildcards)函数的实际意义是返回键值对中的值(路径)
- all或者最后一步中利用expand定义输出文件的通配符内容
4 文件输入
看了一篇文章里面提到了利用文件建立数据索引的方法。
首先创建一个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,b
samples_table.loc[int(wildcards.sample), "path"]
返回地址
sample=samples_table.index
返回索引号
在b中通过id确定a的输出文件然后自动获得对应地址