Commit 55168169 authored by Blaise Li's avatar Blaise Li
Browse files

Minor tweaks.

parent 58285aac
......@@ -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"
......
......@@ -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}
"""
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment