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

Less stringent remapping parameters.

parent e8104160
No related branches found
No related tags found
No related merge requests found
......@@ -26,7 +26,7 @@ mpl.rcParams["font.family"] = "sans-serif"
import pandas as pd
import matplotlib.pyplot as plt
from libhts import median_ratio_to_pseudo_ref_size_factors, plot_histo
from libhts import make_empty_bigwig, median_ratio_to_pseudo_ref_size_factors, plot_histo
from libworkflows import get_chrom_sizes, cleanup_and_backup
from libworkflows import last_lines, ensure_relative, SHELL_FUNCTIONS, warn_context
from libworkflows import feature_orientation2stranded, read_feature_counts, sum_feature_counts
......@@ -41,7 +41,8 @@ alignment_settings = {
# Lower stringency settings, to remap the unmapped
realignment_settings = {
# Allow more mismatches in the seed
"bowtie2": "--local -L 6 -i S,1,0.8 -N 1",
# Reduce minimal mismatch and gap open penalties
"bowtie2": "--local -L 6 -i S,1,0.8 -N 1 --mp 6,1 --rdg 4,3",
# TODO: Find how to be less stringent with crac
"crac": "-k 20 --stranded --use-x-in-cigar"}
......@@ -90,7 +91,7 @@ COND_COLUMNS = pd.DataFrame(CONDITIONS).assign(
LIB_TYPE = config["lib_type"]
#TRIMMERS = ["fastx_clipper"]
TRIMMERS = ["cutadapt"]
COUNTERS = ["feature_count"]
#COUNTERS = ["feature_count"]
ORIENTATIONS = ["fwd", "rev", "all"]
WITH_ADAPT = ["adapt_deduped", "adapt_nodedup"]
POST_TRIMMING = ["noadapt_deduped"] + WITH_ADAPT
......@@ -669,33 +670,32 @@ rule make_normalized_bigwig:
size = pd.read_table(input.norm_file).T.loc[wildcards.norm][0] / 1000000
#scale = 1 / pd.read_table(input.summary, index_col=0).loc[
# wildcards.norm_file].loc[f"{wildcards.lib}_{wildcards.rep}"]
assert size > 0
# TODO: make this a function of deeptools version
no_reads = """Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.
"""
zero_bytes = """needLargeMem: trying to allocate 0 bytes (limit: 100000000000)
bam2bigwig.sh: bedGraphToBigWig failed
"""
try:
shell("""
bam2bigwig.sh {input.bam} {params.genome_binned} \\
{wildcards.lib}_{wildcards.rep} {wildcards.orientation} %s \\
%f {output.bigwig_norm} \\
> {log.log} 2> {log.err} \\
|| error_exit "bam2bigwig.sh failed"
""" % (LIB_TYPE[-1], size))
except CalledProcessError as e:
if last_lines(log.err, 2) in {no_reads, zero_bytes}:
bw_out = pyBigWig.open(output.bigwig_norm, "w")
bw_out.addHeader(list(chrom_sizes.items()))
for (chrom, chrom_len) in bw_out.chroms().items():
bw_out.addEntries(chrom, 0, values=np.nan_to_num(np.zeros(chrom_len)[0::10]), span=10, step=10)
bw_out.close()
#with open(output.bigwig_norm, "w") as bwfile:
# bwfile.write("")
else:
raise
assert size >= 0, f"{size} is not positive"
if size == 0:
make_empty_bigwig(output.bigwig_norm, chrom_sizes)
else:
# TODO: make this a function of deeptools version
no_reads = """Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.
"""
zero_bytes = """needLargeMem: trying to allocate 0 bytes (limit: 100000000000)
bam2bigwig.sh: bedGraphToBigWig failed
"""
try:
shell("""
bam2bigwig.sh {input.bam} {params.genome_binned} \\
{wildcards.lib}_{wildcards.rep} {wildcards.orientation} %s \\
%f {output.bigwig_norm} \\
> {log.log} 2> {log.err} \\
|| error_exit "bam2bigwig.sh failed"
""" % (LIB_TYPE[-1], size))
except CalledProcessError as e:
if last_lines(log.err, 2) in {no_reads, zero_bytes}:
make_empty_bigwig(output.bigwig_norm, chrom_sizes)
#with open(output.bigwig_norm, "w") as bwfile:
# bwfile.write("")
else:
raise
onsuccess:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment