diff --git a/PRO-seq/PRO-seq.snakefile b/PRO-seq/PRO-seq.snakefile index 96be2f378569298f184d6f6f1c8e7f559228f2c5..8155de603a54802996defe710fa64a2a56db07b3 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 4fb828b0cebb432aebbb78bff43750e9cf434139..2321af30200539013f2c005ac7d0e0261b850f91 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} """