Skip to content
Snippets Groups Projects
Commit 7569bf56 authored by Blaise Li's avatar Blaise Li
Browse files

Tried to use bedtools intersect for counting.

It currently fails when libraries have too many rRNA reads:
https://github.com/arq5x/bedtools2/issues/565

Also added lfc distribution and improved MA plot.
parent 7110f2c7
No related branches found
No related tags found
No related merge requests found
from snakemake.shell import shell
from libworkflows import SHELL_FUNCTIONS
# Define functions to be used in shell portions
shell.prefix(SHELL_FUNCTIONS)
cmd = """
# make bedtools crash if it uses too much memory
# http://cr.yp.to/daemontools/softlimit.html
# sudo apt install daemontools
converter="/pasteur/entites/Mhe/Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/genes_id2name.pickle"
# Put fifos in a temporary directory
workdir=$(mktemp -d)
annot="${{workdir}}/annot.bed"
# Doesn't work if fifo:
# The annot file may have only 3 columns instead of 6
# Maybe it is read before the line is complete?
#mkfifo ${{annot}}
genome_file="${{workdir}}/chrom_sizes.txt"
mkfifo ${{genome_file}}
cleanup()
{{
rm -rf ${{workdir}}
}}
> {snakemake.log.err}
# Write the fifo(s) in the background
samtools view -H {snakemake.input.sorted_bam} | awk '$1 == "@SQ" {{print $2"\t"$3}}' | sed 's/SN://g' | sed 's/LN://' > ${{genome_file}} 2>> {snakemake.log.err} || {{ cleanup && error_exit "making ${{genome_file}} failed"; }} &
#gff_merge_transcripts.sh {snakemake.params.annot} > ${{annot}} 2>> {snakemake.log.err} || {{ cleanup && error_exit "making merged transcript failed"; }}
#cmd="bedtools intersect -a ${{annot}} -b {snakemake.input.sorted_bam} -sorted -g ${{genome_file}} -c | tee {snakemake.output.counts} | id2name.py ${{converter}} > {snakemake.output.counts_converted}"
cmd="softlimit -a 4294967296 bedtools intersect -a {snakemake.params.annot} -b {snakemake.input.sorted_bam} -sorted -g ${{genome_file}} -c | tee {snakemake.output.counts} | id2name.py ${{converter}} > {snakemake.output.counts_converted}"
echo ${{cmd}} > {snakemake.log.log}
eval ${{cmd}} 1>> {snakemake.log.log} 2>> {snakemake.log.err} || {{ cleanup && error_exit "bedtools intersect failed"; }}
cleanup
"""
shell(cmd)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment