У меня есть работающий файл Snakefile для разделенного наследования, использующий несколько файлов кроватей.Это создает идеальный список заданий с использованием snakemake -np, поэтому этот файл требует лишь незначительной настройки (надеюсь!).
Моя проблема возникает в правиле merge_peaks ниже.
На этом этапе у меня есть 25 файлов кроватей, и мне нужно выполнить 5 вызовов merge_peaks, по одному вызову для каждого расширения ext=[100,200,300,400,500], поэтому мне нужно только 5 файлов кроватей, содержащих соответствующее расширение, для вызовакаждый раз.
Например, для следующего 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
Вот Snakefile:
# 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: правило для использования множества входов для одного выхода с несколькими подгруппами
Я люблюЗмея, но мой питон немного рискованный.
Большое спасибо.