diff --git a/tritimap/Snakefile b/tritimap/Snakefile index a6500bf..6de5deb 100644 --- a/tritimap/Snakefile +++ b/tritimap/Snakefile @@ -186,6 +186,8 @@ if config['module'] == "only_mapping": join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates()) + "_snpindex_input.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_output.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indelinfo.txt"), + join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indel.bed"), + join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_snp.bed"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_region.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_raw_region.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_snpinfo.txt"), @@ -211,6 +213,8 @@ elif config['module'] == "all": join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates()) + "_snpindex_input.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_output.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indelinfo.txt"), + join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indel.bed"), + join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_snp.bed"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_region.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_raw_region.txt"), join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_snpinfo.txt"), diff --git a/tritimap/rules/qtlseqr_plot.smk b/tritimap/rules/qtlseqr_plot.smk index 390ed39..5a54c19 100644 --- a/tritimap/rules/qtlseqr_plot.smk +++ b/tritimap/rules/qtlseqr_plot.smk @@ -25,4 +25,4 @@ rule QTLseqr_plot: log: join(dir_path+"/logs", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr.log") script: - "../scripts/qtlseqr_plot.R" \ No newline at end of file + "../scripts/qtlseqr_plot.R" diff --git a/tritimap/rules/vcf_pretreat.smk b/tritimap/rules/vcf_pretreat.smk index 8f97fba..f245efb 100644 --- a/tritimap/rules/vcf_pretreat.smk +++ b/tritimap/rules/vcf_pretreat.smk @@ -26,18 +26,23 @@ rule vcf2tab: fi """ -rule getRegionIndel: +rule getRegionIndelandSNPbed: input: indel = join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates()) + "_indel.txt"), + snp = join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_snpinfo.txt"), region=join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_region.txt") params: dirpath = join(dir_path+"/06_regionout") output: - join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indelinfo.txt") - message: "\nGet indel information in candidate region. Input file: {input}\n" + indelinfo = join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indelinfo.txt"), + indelbed = join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_indel.bed"), + snpbed = join(dir_path+"/06_regionout", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr_filter_snp.bed") + message: "\nGet candidate region indel and snp bed file. Input file: {input.indel}; {input.snp}; {input.region}\n" log: - join(dir_path+"/logs", "_".join(samples.bulk.drop_duplicates()) + "_getregionindel.log") + join(dir_path+"/logs", "_".join(samples.bulk.drop_duplicates()) + "_get_region_indel_snp_bed.log") shell:""" cat {input.region} | grep -v 'start' | cut -f1,3,4 | sort -k1,1 -k2,2n > {params.dirpath}/temp.candidateregion.indel.bed - bedmap --skip-unmapped --delim '\t' --header --echo {input.indel} {params.dirpath}/temp.candidateregion.indel.bed | cat <(head -n1 {input.indel}) - > {output} && rm {params.dirpath}/temp.candidateregion.indel.bed + bedmap --skip-unmapped --delim '\t' --header --echo {input.indel} {params.dirpath}/temp.candidateregion.indel.bed | cat <(head -n1 {input.indel}) - | cut -f1,3- | sed 's/^#CHROM/CHROM/' > {output.indelinfo} && rm {params.dirpath}/temp.candidateregion.indel.bed + cat {output.indelinfo} | grep -v '^CHROM'| awk 'BEGIN{{OFS="\t"}}{{print $1,$2-1,$2,$3,$4}}' |sort -k1,1 -k2,2n > {output.indelbed} + cat {input.snp} | grep -v '^CHROM'| awk 'BEGIN{{OFS="\t"}}{{print $1,$2-1,$2,$3,$4}}' |sort -k1,1 -k2,2n > {output.snpbed} """ \ No newline at end of file