Merge branch 'master' of https://github.com/chaissonlab/SegDupAnnotation

Author: mchaisso

RepeatMaskGenome.snakefile

@@ -10,7 +10,9 @@ configfile: "sd_analysis.json"
 if "asm" not in config and "assembly" in config:
     config["asm"] = config["assembly"]
 
-asmFai=config["asm"] + ".fai"
+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))]
@@ -27,7 +29,7 @@ rule all:
 
 rule SplitContig:
     input:
-        asm=config["asm"]
+        asm=assembly
     output:
         contig="split/to_mask.{index}.fasta",
     params:
@@ -78,26 +80,31 @@ rule SpecialMaskContig:
         tmpdir=tempDir,
         sd=SD
     shell:"""
-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 [ {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
+  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
-  cp $TEMP/to_mask.\"{wildcards.index}\".fasta.* t2t/ || true
+  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
-rm -rf $TEMP
-
 """
 
 rule MergeMaskerRuns:

SegDupAnalysis.snakefile

@@ -34,7 +34,7 @@ pos=[]
 subs=["all", "high_ident"]
 
 
-localrules: all, AnnotateResolvedTandemDups, GetUniqueGencodeUnresolvedDupGenes,  IntersectGenesWithFullSDList, FullDupToBed12, FullDupToLinks, MakeWMBed, MaskFile, ConvertHMMCopyNumberToCollapsedDuplications, SortSedef, FilterSedef, CountMaskedSedef, RemoveSedefTooMasked, MakeSedefGraph, MakeSedefGraphTable, FilterByGraphClusters, FullDupToBed12, FiltDupToBed12, GetUniqueGencodeUnresolvedDupGenesCN, GetUniqueGencodeUnresolvedDupGenes, GetGencodeMulticopy, GetGencodeMappedInDup, GetSupportedMulticopy,FindResolvedDuplicatedGenes, Bed12ToBed6, CombineGenesWithCollapsedDups, CombineDuplicatedGenes, MinimapGeneModelBed, FilterGencodeBed12, FindGenesInResolvedDups, SelectOneIsoform, SplitSplicedAndSingleExon, AnnotateLowCoverageFlanks, UnionMasked,GetNamedFasta, SelectDups, SortDups, GetDepthOverDups, FilterLowDepthDups, GetFullGeneCountTable, AddCollapsedGenes, GetCombinedTable, SelectDupsOneIsoform, GetFinalMerged, DupsPerContig, GetAllMultiGenes, AnnotateHighIdentity, GetTotalMasked, AnnotateResolvedTandemDups, GeneCountFact, GetFullGeneCountTable, FilterMultiExonBed, MappedSamIdentityDups, RemoveOriginal, RemoveBams, MakeSedefIntv, HighestIdentPairs, SelectHighIdent, GetCollapseByRange, GetCollapsedMask, GetCN
+localrules: all, AnnotateResolvedTandemDups, GetUniqueGencodeUnresolvedDupGenes,  IntersectGenesWithFullSDList, FullDupToBed12, FullDupToLinks, MakeWMBed, MaskFile, ConvertHMMCopyNumberToCollapsedDuplications, SortSedef, FilterSedef, CountMaskedSedef, RemoveSedefTooMasked, MakeSedefGraph, MakeSedefGraphTable, FilterByGraphClusters, FullDupToBed12, FiltDupToBed12, GetUniqueGencodeUnresolvedDupGenesCN, GetUniqueGencodeUnresolvedDupGenes, GetGencodeMulticopy, GetGencodeMappedInDup, GetSupportedMulticopy,FindResolvedDuplicatedGenes, Bed12ToBed6, CombineGenesWithCollapsedDups, CombineDuplicatedGenes, MinimapGeneModelBed, FilterGencodeBed12, FindGenesInResolvedDups, SelectOneIsoform, SplitSplicedAndSingleExon, AnnotateLowCoverageFlanks, UnionMasked,GetNamedFasta, SelectDups, SortDups, GetDepthOverDups, FilterLowDepthDups, GetFullGeneCountTable, AddCollapsedGenes, GetCombinedTable, SelectDupsOneIsoform, GetFinalMerged, DupsPerContig, GetAllMultiGenes, AnnotateHighIdentity, GetTotalMasked, AnnotateResolvedTandemDups, GeneCountFact, GetGeneCountTableAbbreviatedNames, FilterMultiExonBed, MappedSamIdentityDups, RemoveOriginal, RemoveBams, MakeSedefIntv, HighestIdentPairs, SelectHighIdent, GetCollapseByRange, GetCollapsedMask, GetCN
 
 
 
@@ -64,7 +64,8 @@ rule all:
         comb_with_depth="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth",
         fact="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.fact",
         gene_count="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.gene_count",
-        gene_count_2column="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.gene_count_multi_single",            
+        gene_count_2column="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.gene_count_multi_single",
+        gene_count_abbrvNames="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.gene_count.abbrv_names",
         asmMask=expand("{asm}.count_masked", asm=["assembly.orig.fasta", "assembly.masked.fasta", "assembly.repeat_masked.fasta", "assembly.union_masked.fasta"]),
         uniqueDupGenes="gencode.mapped.bam.bed12.dups.unique",
         sedef_high_uniq="sedef_out/all/final.sorted.bed.uniq.high",        
@@ -1428,7 +1429,7 @@ rule MappedSamIdentity:
         grid_opts=config["grid_small"],
         sd=SD,
     shell:"""
- {params.sd}/hmcnc/HMM/samToBed {input.mappedsam} --reportAccuracy > {output.mappedsambed}
+ {params.sd}/hmcnc/src/samToBed {input.mappedsam} --reportAccuracy > {output.mappedsambed}
 """
 
 rule AddDepthCopyNumber:
@@ -1466,7 +1467,8 @@ rule AnnotateOriginal:
     output:
         annot="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.sam.bed.dups.annot_orig",
     params:
-        sd=SD
+        sd=SD,
+        grid_opts=config["grid_medium"] # KEON
     shell:"""
 {params.sd}/AnnotateOriginal.py {input.mappedsambeddups} > {output.annot}
 """
@@ -1688,6 +1690,19 @@ cat {input.depth_filt} | awk '{{ if (NR == 1) {{ print "gene\\tresolved\\tcollap
 cat {output.gene_count_2column} | awk '{{ if (NR == 1) {{ print "gene\\tcopies";}} else {{  print $1"\\t"$2+$3;}} }}' > {output.gene_count}
 """
 
+rule GetGeneCountTableAbbreviatedNames:
+    input:
+        gene_count="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.gene_count",
+    output:
+        gene_count_abbrvNames="gencode.mapped.bam.bed12.multi_exon.fasta.named.mm2.dups.one_isoform.txt.combined.depth.filt.gene_count.abbrv_names",
+    params:
+        sd=SD
+    shell:"""
+awk 'BEGIN {{ FS = "[ |\t]"}} \
+NR == 1 {{print}} \
+NR > 1 {{print $1,"\t" $3}}' {input.gene_count} > {output.gene_count_abbrvNames}
+"""
+
 rule cramBam:
     input:
         bam=config['bam'],