From 551681691735931e01a2459cedd9fe4a0aeb8cc5 Mon Sep 17 00:00:00 2001 From: Blaise Li <blaise.li__git@nsup.org> Date: Thu, 11 Jan 2018 09:40:20 +0100 Subject: [PATCH] Minor tweaks. --- PRO-seq/PRO-seq.snakefile | 23 ++++++----------------- small_RNA-seq/small_RNA-seq.snakefile | 8 ++------ 2 files changed, 8 insertions(+), 23 deletions(-) diff --git a/PRO-seq/PRO-seq.snakefile b/PRO-seq/PRO-seq.snakefile index 96be2f3..8155de6 100644 --- a/PRO-seq/PRO-seq.snakefile +++ b/PRO-seq/PRO-seq.snakefile @@ -433,19 +433,14 @@ rule fuse_bams: rule compute_coverage: input: - bam = rules.fuse_bams.output.sorted_bam, + sorted_bam = rules.fuse_bams.output.sorted_bam, output: coverage = OPJ(output_dir, "{trimmer}", aligner, "mapped_C_elegans", "{lib}_{rep}_on_C_elegans_coverage.txt"), params: genomelen = genomelen, - log: - coverage_log = OPJ(log_dir, "{trimmer}", "compute_coverage", "{lib}_{rep}.log"), - coverage_err = OPJ(log_dir, "{trimmer}", "compute_coverage", "{lib}_{rep}.err"), shell: """ - # genomelen=$(samtools view -H {input} | grep "^@SQ" | awk -F ":" '{{sum+=$3}} END {{print sum}}') - # genomelen=$(samtools view -H {input} | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{{sum+=$1}} END {{print sum}}') - bases=$(samtools depth {input} | awk '{{sum += $3}} END {{print sum}}') || error_exit "samtools depth failed" + bases=$(samtools depth {input.sorted_bam} | awk '{{sum += $3}} END {{print sum}}') || error_exit "samtools depth failed" python3 -c "print(${{bases}} / {params.genomelen})" > {output.coverage} """ @@ -540,17 +535,11 @@ rule htseq_count_reads: annot = biotype2annot, message: "Counting {wildcards.orientation} {wildcards.biotype} reads for {wildcards.lib}_{wildcards.rep} with htseq-count." + benchmark: + OPJ(log_dir, "{trimmer}", "htseq_count_reads", "{lib}_{rep}_{biotype}_{orientation}_benchmark.txt") log: - log = OPJ(log_dir, "{trimmer}", "htseq_count_reads", "{lib}_{rep}.log"), - err = OPJ(log_dir, "{trimmer}", "htseq_count_reads", "{lib}_{rep}.err") - #shell: - # """ - # annot="/pasteur/entites/Mhe/Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/{wildcards.biotype}.gtf" - # converter="/pasteur/entites/Mhe/Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/genes_id2name.pickle" - # cmd="htseq-count -f bam -s {params.stranded} -a 0 -t transcript -i gene_id {input.sorted_bam} ${{annot}} | tee {output.counts} | id2name.py ${{converter}} > {output.counts_converted}" - # echo ${{cmd}} - # eval ${{cmd}} > {log} - # """ + log = OPJ(log_dir, "{trimmer}", "htseq_count_reads", "{lib}_{rep}_{biotype}_{orientation}.log"), + err = OPJ(log_dir, "{trimmer}", "htseq_count_reads", "{lib}_{rep}_{biotype}_{orientation}.err") wrapper: "file:///pasteur/homes/bli/src/bioinfo_utils/snakemake_wrappers/htseq_count_reads" diff --git a/small_RNA-seq/small_RNA-seq.snakefile b/small_RNA-seq/small_RNA-seq.snakefile index 4fb828b..2321af3 100644 --- a/small_RNA-seq/small_RNA-seq.snakefile +++ b/small_RNA-seq/small_RNA-seq.snakefile @@ -803,18 +803,14 @@ rule sam2indexedbam: rule compute_coverage: input: - rules.sam2indexedbam.output.sorted_bam, + sorted_bam = rules.sam2indexedbam.output.sorted_bam, output: coverage = OPJ(output_dir, aligner, "mapped_C_elegans", "{lib}_{rep}", "{read_type}_on_C_elegans_coverage.txt"), params: genomelen = genomelen, - # Currently not used - #log: - # log = OPJ(log_dir, "compute_coverage", "{lib}_{rep}_{read_type}.log"), - # err = OPJ(log_dir, "compute_coverage", "{lib}_{rep}_{read_type}.err"), shell: """ - bases=$(samtools depth {input} | awk '{{sum += $3}} END {{print sum}}') || error_exit "samtools depth failed" + bases=$(samtools depth {input.sorted_bam} | awk '{{sum += $3}} END {{print sum}}') || error_exit "samtools depth failed" python3 -c "print(${{bases}} / {params.genomelen})" > {output.coverage} """ -- GitLab