SegDupAnnotation/RepeatMaskGenome.snakefile at master · ChaissonLab/SegDupAnnotation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
import os
import tempfile
import subprocess
import os.path


# Config
configfile: "sd_analysis.json"

if "asm" not in config and "assembly" in config:
    config["asm"] = config["assembly"]

assembly="assembly.orig.fasta"
asmFai="assembly.orig.fasta.fai"
#asmFai=config["asm"] + ".fai"
asmFaiFile=open(asmFai)
contigs=[l.split()[0] for l in asmFaiFile]
strIdx=[str(i) for i in range(0,len(contigs))]

tempDir=config["temp"]
# Snakemake and working directories
SD = os.path.dirname(workflow.snakefile)

rule all:
    input:
        mask="assembly.repeat_masked.fasta",
        maskedGenomeOut="assembly.repeat_masked.fasta.out"


rule SplitContig:
    input:
        asm=assembly
    output:
        contig="split/to_mask.{index}.fasta",
    params:
        contigName=lambda wildcards: contigs[int(wildcards.index)],
        grid_opts=config["grid_small"]
    shell:"""
mkdir -p split
samtools faidx {input.asm} \"{params.contigName}\" > {output.contig}
"""

rule MaskContig:
    input:
        contig="split/to_mask.{index}.fasta"
    output:
        mask="masked/to_mask.{index}.fasta.masked",
        maskOut="masked/to_mask.{index}.fasta.out"
    params:
        grid_opts=config["grid_repeatmasker"],
        repeatLibrary=config["repeat_library"],
        sd=SD
    shell:"""
mkdir -p masked
TEMP="$TMPDIR/$$_$RANDOM/"
mkdir -p $TEMP
cp \"{input.contig}\" \"$TEMP/to_mask.{wildcards.index}.fasta\" && \
pushd $TEMP &&  \
RepeatMasker {params.repeatLibrary} -pa 8  -s -xsmall \"to_mask.{wildcards.index}.fasta\" && \
popd && \
if [ ! -e $TEMP/to_mask.\"{wildcards.index}\".fasta.masked ]; then
  cp split/to_mask.\"{wildcards.index}\".fasta masked/to_mask.\"{wildcards.index}\".fasta.masked
  cp {params.sd}/repeat_masker.out.header masked/to_mask.\"{wildcards.index}\".fasta.out
else
  cp $TEMP/to_mask.\"{wildcards.index}\".fasta.* masked/ || true
fi
#rm -rf $TEMP
"""

rule SpecialMaskContig:
    input:
        mask="masked/to_mask.{index}.fasta.masked",
    output:
        tt="t2t/to_mask.{index}.fasta.masked",
        ttOut="t2t/to_mask.{index}.fasta.out",
    params:
        grid_opts=config["grid_repeatmasker"],
        repeatLibrary=config["t2t_repeat_library"],
        sd=SD
    shell:"""
if [ {params.repeatLibrary} != "na" ]
then
  mkdir -p t2t
  TEMP="$TMPDIR/$$_$RANDOM/"
  mkdir -p $TEMP
  #
  # Copy input file to temp dir that should have fast IO
  #
  {params.sd}/hardmask {input.mask} $TEMP/to_mask.{wildcards.index}.fasta && \
  pushd $TEMP &&  \
  RepeatMasker -nolow -libdir /home1/mchaisso/miniconda3/share/RepeatMasker/Libraries/Appended -species human -pa 8 -s -xsmall \"to_mask.{wildcards.index}.fasta\" && \
  popd && \

  if [ ! -e $TEMP/to_mask.\"{wildcards.index}\".fasta.masked ]; then
    ls -l masked/to_mask.\"{wildcards.index}\".fasta.masked
    cp masked/to_mask.\"{wildcards.index}\".fasta.masked t2t/to_mask.\"{wildcards.index}\".fasta.masked
    head -3 masked/to_mask.\"{wildcards.index}\".fasta.out > t2t/to_mask.\"{wildcards.index}\".fasta.out
  else
    cp $TEMP/to_mask.\"{wildcards.index}\".fasta.* t2t/ || true
  fi
  rm -rf $TEMP
else
  echo "*** rule SpecialMaskContig made no change (skipped) ***" > t2t/to_mask.\"{wildcards.index}\".fasta.out
  cp {input.mask} t2t/to_mask.\"{wildcards.index}\".fasta.masked >> t2t/to_mask.\"{wildcards.index}\".fasta.out || true
fi
"""

rule MergeMaskerRuns:
    input:
        humLib="masked/to_mask.{index}.fasta.masked",
        humLibOut="masked/to_mask.{index}.fasta.out",
        t2tLib="t2t/to_mask.{index}.fasta.masked",
        t2tLibOut="t2t/to_mask.{index}.fasta.out"
    output:
        comb="comb/to_mask.{index}.fasta.masked",
        combOut="comb/to_mask.{index}.fasta.out",
    params:
        grid_opts=config["grid_small"],
        sd=SD
    shell:"""
{params.sd}/comask {output.comb} {input.humLib} {input.t2tLib}
echo {input.humLibOut} > comb/to_mask.{wildcards.index}.names
echo {input.t2tLibOut} >> comb/to_mask.{wildcards.index}.names
{params.sd}/RepeatMasking/AppendOutFile.py {output.combOut} comb/to_mask.{wildcards.index}.names
"""

rule CombineMask:
    input:
        maskedContigs=expand("comb/to_mask.{index}.fasta.masked", index=strIdx),
        maskedContigsOut=expand("comb/to_mask.{index}.fasta.out", index=strIdx)
    output:
        maskedGenome="assembly.repeat_masked.fasta",
        maskedGenomeOut="assembly.repeat_masked.fasta.out"
    params:
        grid_opts=config["grid_small"],
        outFileNames=" ".join(["comb/to_mask.{}.fasta.out".format(i) for i in strIdx]),
        inFastaNames=" ".join(["comb/to_mask.{}.fasta.masked".format(i) for i in strIdx]),
        sd=SD
    shell:"""
cat {params.inFastaNames} > {output.maskedGenome}
{params.sd}/RepeatMasking/AppendOutFile.py {output.maskedGenomeOut} {params.outFileNames}
"""