new: add snp and index bed format file

fei0810 · Apr 8, 2021 · 04034cf · 04034cf
1 parent 1a1181e
commit 04034cf
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 6 deletions.
diff --git 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
@@ -25,4 +25,4 @@ rule QTLseqr_plot:
  log:
  join(dir_path+"/logs", "_".join(samples.bulk.drop_duplicates())+ "_qtlseqr.log")
  script:
- "../scripts/qtlseqr_plot.R"
+ "../scripts/qtlseqr_plot.R"
diff --git 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}
  """