【发布时间】：2018-12-19 14:03:27
【问题描述】：
我有一个使用多个床文件的分区遗传性的有效 Snakefile。这会使用snakemake -np 生成一个完美的作业列表，因此该文件只需要稍作调整（我希望！）。
我的问题出现在下面的merge_peaks 规则中。
现阶段，我有25个bed文件，需要运行5次merge_peaks的调用，每个扩展名ext=[100,200,300,400,500]调用一次，所以每次只需要调用包含相关扩展名的5个bed文件即可。
例如对于下面的merge_peaks输出文件peak_files/Fullard2018_peaks.mrgd.blrm.100.bed，我只需要以下5个带有ext=100的床文件作为输入：
peak_files/fullard2018_NpfcATAC_1.blrm.100.bed
peak_files/fullard2018_NpfcATAC_2.blrm.100.bed
peak_files/fullard2018_NpfcATAC_3.blrm.100.bed
peak_files/fullard2018_NpfcATAC_4.blrm.100.bed
peak_files/fullard2018_NpfcATAC_5.blrm.100.bed
这是我的配置文件：
samples:
fullard2018_NpfcATAC_1:
fullard2018_NpfcATAC_2:
fullard2018_NpfcATAC_3:
fullard2018_NpfcATAC_4:
fullard2018_NpfcATAC_5:
index: /home/genomes_and_index_files/hg19.chrom.sizes
这里是蛇文件：
# read config info into this namespace
configfile: "config.yaml"
print (config['samples'])
rule all:
input:
expand("peak_files/{sample}.blrm.{ext}.bed", sample=config['samples'], ext=[100,200,300,400,500]),
expand("LD_annotation_files/Fullard2018.{ext}.{chr}.l2.ldscore.gz", sample=config['samples'], ext=[100,200,300,400,500], chr=range(1,23))
rule annot2bed:
input:
folder = "Reference/baseline"
params:
file = "Reference/baseline/baseline.{chr}.annot.gz"
output:
"LD_annotation_files/baseline.{chr}_no_head.bed"
shell:
"zcat {params.file} | tail -n +2 |awk -v OFS=\"\t\" '{{print \"chr\"$1, $2-1, $2, $3, $4}}' "
"| sort -k1,1 -k2,2n > {output}"
rule extend_bed:
input:
"peak_files/{sample}_peaks.blrm.narrowPeak"
output:
"peak_files/{sample}.blrm.{ext}.bed"
params:
ext = "{ext}",
index = config["index"]
shell:
"bedtools slop -i {input} -g {params.index} -b {params.ext} > {output}"
rule merge_peaks:
input:
expand("peak_files/{sample}.blrm.{ext}.bed", sample = config['samples'], ext=[100,200,300,400,500])
output:
"peak_files/Fullard2018_peaks.mrgd.blrm.{ext}.bed"
shell:
"cat {input} | bedtools sort -i stdin | bedtools merge -i stdin > {output}"
rule intersect_mybed:
input:
annot = rules.annot2bed.output,
mybed = rules.merge_peaks.output
output:
"LD_annotation_files/Fullard2018.{ext}.{chr}.annot.gz"
params:
uncompressed = "LD_annotation_files/Fullard2018.{ext}.{chr}.annot"
shell:
"echo -e \"CHR\tBP\tSNP\tCM\tANN\" > {params.uncompressed}; "
"/share/apps/bedtools intersect -a {input.annot} -b {input.mybed} -c | "
"sed 's/^chr//g' | awk -v OFS=\"\t\" '{{print $1, $2, $4, $5, $6}}' >> {params.uncompressed}; "
"gzip {params.uncompressed}"
rule ldsr:
input:
annot = "LD_annotation_files/Fullard2018.{ext}.{chr}.annot.gz",
bfile_folder = "Reference/1000G_plinkfiles",
snps_folder = "Reference/hapmap3_snps"
output:
"LD_annotation_files/Fullard2018.{ext}.{chr}.l2.ldscore.gz"
conda:
"envs/p2-ldscore.yaml"
params:
bfile = "Reference/1000G_plinkfiles/1000G.mac5eur.{chr}",
ldscores = "LD_annotation_files/Fullard2018.{ext}.{chr}",
snps = "Reference/hapmap3_snps/hm.{chr}.snp"
log:
"logs/LDSC/Fullard2018.{ext}.{chr}_ldsc.txt"
shell:
"export MKL_NUM_THREADS=2;" # Export arguments are workaround as ldsr uses all available cores
"export NUMEXPR_NUM_THREADS=2;" # Numbers must match cores parameter in cluster config
"Reference/ldsc/ldsc.py --l2 --bfile {params.bfile} --ld-wind-cm 1 "
"--annot {input.annot} --out {params.ldscores} --print-snps {params.snps} 2> {log}"
当前发生的情况是，每次调用的所有 25 个床文件都被输入到合并峰值规则中，而我每次只需要输入带有相关扩展名的 5 个。我正在努力研究如何正确使用扩展功能或替代方法，以便在每次调用规则时仅包含和合并相关的床文件。
我认为这个问题会问类似的问题，但不是我想要的，因为它不使用配置文件 - Snakemake: rule for using many inputs for one output with multiple sub-groups
我喜欢 Snakemake，但我的蟒蛇有点冒险。
非常感谢。
【问题讨论】：
标签： expand many-to-one snakemake