shell:
    """
    set -e
    set -x

    export OMPI_MCA_btl=^openib
    export OMP_NUM_THREADS={threads}
    export OMP_PROC_BIND=FALSE

    {params.categ_db} --cluseq_db {params.cluseqdb} \
                      --step "{params.step}" \
                      --mpi_runner "{params.mpi_runner}" \
                      --mmseqs {params.mmseqs_bin} \
                      --aln {params.famsa_bin} \
                      --hhsuite {params.hhsuite} \
                      --reformat {params.reformat} \
                      --consensus {params.consensus} \
                      --hhmake {params.hhmake} \
                      --outdir {params.outdir} \
                      --idir {params.idir} \
                      --clu_hhm {output.clu_hhm} \
                      --threads {threads} 2>{log.err} 1>{log.out}

    """

run:
    shell("echo 'CATEGORY DB DONE'")

run:
    shell("echo 'CATEGORY REFINEMENT DONE'")

run:
    shell("echo 'COMMUNITIES INFERENCE DONE'")

run:
    shell("echo 'CLUSTER CLASSIFICATION DONE'")

shell:
  """
  set -e
  set -x

  export OMPI_MCA_btl=^openib
  export OMP_NUM_THREADS={threads}
  export OMP_PROC_BIND=FALSE

  # HHblits all-clusters vs all-clusters for each category
  CATEG=$(echo -e "eu\ngu\nkwp\nk")

  IN=$(dirname {input.k})
  OUT=$(dirname {output.comm})

  for categ in $CATEG; do

    RES=${{OUT}}/${{categ}}_hhblits.tsv

    if [[ ! -s ${{RES}} ]]; then
        if [[ ! -s {params.hhb_tmp_db}.index  ]]; then
            {params.vmtouch} -f ${{IN}}/${{categ}}*
            {params.mpi_runner} {params.hhblits} -i ${{IN}}/${{categ}}_hhm_db \
                                             -o {params.hhb_tmp_db} \
                                             -n 2 -cpu 1 -v 0 \
                                             -d ${{IN}}/${{categ}}
            mv {params.hhb_tmp_db}.ffdata {params.hhb_tmp_db}
            mv {params.hhb_tmp_db}.ffindex {params.hhb_tmp_db}.index
        fi
        {params.mpi_runner} {params.mmseqs_bin} apply {params.hhb_tmp_db} {params.hhb_tmp_db}.parsed \
            --threads 1 \
            -- {params.hhparse} 2>{log.err}

            sed -e 's/\\x0//g' {params.hhb_tmp_db}.parsed > ${{RES}} 2>{log.err}

            rm -rf {params.hhb_tmp_db} {params.hhb_tmp_db}.index {params.hhb_tmp_db}.dbtype {params.hhb_tmp_db}.parsed* {params.hhb_tmp_db}.ff*
    fi

  done
  """

shell:
  """
  OUT=$(dirname {output.comm})

  if [ ! -s {params.pfam_clan} ]; then
    echo "Dowloading Pfam-A clan information"
    wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam34.0/Pfam-A.clans.tsv.gz -O {params.pfam_clan}
  fi

  # Start cluster community inference
  ./{params.get_comm} -c {params.comm_config} 1>{log.out} 2>{log.err}

  rm -rf ${{OUT}}/tmp

  if [[ {params.db_mode} == "memory" ]]; then
     rm {params.pfam_clan}
  fi

  """

run:
  shell("echo 'COMMUNITIES INFERENCE DONE'")

shell:
    """
    set -e
    set -x

    export OMPI_MCA_btl=^openib
    export OMP_NUM_THREADS={threads}
    export OMP_PROC_BIND=FALSE

    {params.igraph_lib}
    {params.parasail_lib}

    # Run compositional validation in mpi-mode
    if [[ {params.module} == "creation" ]]; then
        DB={params.clseqdb}
    else
        DB={params.new_clseqdb}
    fi

    {params.mpi_runner} {params.mmseqs_bin} apply ${{DB}} {params.outdb} \
        -- {params.cvals} --derep {params.mmseqs_bin} \
                   --msa {params.famsa_bin} \
                   --msaeval {params.odseq_bin} \
                   --ssn {params.parasail_bin} \
                   --gfilter {params.filterg} \
                   --gconnect {params.isconnect} \
                   --seqp {params.seqtk_bin} \
                   --datap {params.datamash} \
                   --stats {params.stats} \
                   --outdir {params.outdir} \
                   --slots {threads} \
                   --threads {threads} 2>{log.err} 1>{log.out}

    # Collect results:
    # collect cluster main compositional validation stats and cluster rejected (bad-aligned) ORFs
    {params.collect} {params.stat_dir} {output.cl_cval} {output.cval_rej} {params.parallel_bin} {params.threads_collect} 2>>{log.err} 1>>{log.out}

    rm -rf {params.outdb} {params.outdb}.index {params.outdb}.dbtype
    rm -rf {params.outdir}/stats {params.outdir}/log_vals {params.outdir}/rejected

    """

run:
    shell("echo 'COMPOSITIONAL VALIDATION DONE'")

shell:
    """

    set -x
    set -e

    # Summary table with cluster db origin (original/shared/new)
    awk -v N={params.data_name} '{{print $1,N,$3}}' {params.or_clu_origin} > {params.clu_origin}
    sed -i 's/ /\\t/g' {params.clu_origin}

    # All gene headers and partiality information
    gzip -c {params.or_partial}  > {params.partial}

    # Spurious and shadow genes information:
    gzip -c {params.or_sp_sh} > {params.sp_sh}

    # All gene Pfam annotations:
    if [ ! -f {params.multi_annot} ]; then
        mv {params.or_multi_annot} {params.multi_annot}
    fi

    # All cluster category annotation files
    ODIR=$(dirname {input.iclu_cat})

    # Combine with new ones
    gzip -c ${{ODIR}}/K_annotations.tsv > {params.rdir}/K_annotations.tsv.gz
    gzip -c ${{ODIR}}/KWP_annotations.tsv > {params.rdir}/KWP_annotations.tsv.gz
    gzip -c ${{ODIR}}/GU_annotations.tsv > {params.rdir}/GU_annotations.tsv.gz

    # Integrated set of cluster categories
    gzip -c {input.iclu_cat} > {output.clu_cat}

    # and the cluster genes
    if [ {params.singl} == true ]; then
        awk -vOFS="\\t" '{{print $2,$3,$1}}' {params.s_categ} | gzip -c > {params.singl_cl_gene_categ}
    fi
    cp {params.or_clu_gene} {params.clu_gene}

    # Integrated set of cluster communities
    gzip -c {input.iclu_com} > {output.clu_com}

    # New integarted cluster HHMs DB and mmseqs profiles
    mkdir -p {params.clu_hhm}
    cp -r {input.iclu_hhm}* {params.clu_hhm}/
    {params.mmseqs_bin} convertprofiledb {input.iclu_hhm} {params.clu_hhm}/clu_hmm_db \
                                         --threads {threads} 2>{log.err}

    # Integrated set of high quality (HQ) clusters
    gzip -c {input.ihq_clu} > {output.hq_clu}

    {params.parser} --clu_or {params.clu_origin} \
                      --cat {output.clu_cat} \
                      --clu_info {params.clu_info} \
                      --comm {output.clu_com} \
                      --hq_clu {output.hq_clu} \
                      --k_annot {params.multi_annot} \
                      --is_singl {params.singl} \
                      --s_cat {params.s_categ} \
                      --anvio {params.data_stage} \
                      --res {output.clu_out_tbl} \
                      --threads {threads} 2>{log.err} 1>{log.out}

    gzip {params.clu_origin}

    # Integrated cluster summary information
    gzip -c {input.iclu_stats}  > {output.clu_stats}

    """

run:
    shell("echo 'INTEGRATED DB DONE'")

shell:
    """
    set -e
    set -x

    # Pfam list common domain terms
    if [ ! -s {params.pfam_shared_terms} ]; then
        echo "Dowloading Pfam list of shared domain names"
        wget https://figshare.com/ndownloader/files/31127782 -O {params.pfam_shared_terms}
    fi

    ./{params.funct_valr} --input {input.cl_annot} \
                          --pfam_terms {params.pfam_shared_terms} \
                          --output {output.fval_res} \
                          --functions {params.funct_val_fun} \
                          --threads {threads} 2>{log.err} 1>{log.out}

    if [[ {params.db_mode} == "memory" ]]; then
        rm {params.pfam_shared_terms}
    fi

    """

run:
    shell("echo 'FUNCTIONAL VALIDATION DONE'")

shell:
    """
    set -e
    set -x

    ## 1. Pfam-clan info and multi-domain format for the ORF Pfam annotations

    if [ ! -s {params.pfam_clan} ]; then
        echo "Dowloading Pfam-A clan information"
        wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam34.0/Pfam-A.clans.tsv.gz -O {params.pfam_clan}
    fi

    # Add Pfam clan info to the Pfam annotation results
    join -11 -21 <( awk '{{print $1,$3,$8,$9}}' {input.annot} | sort -k1,1) \
    <( zcat {params.pfam_clan} | awk -vFS='\\t' -vOFS='\\t' '{{print $4,$1,$2}}' \
    | awk -vFS='\\t' -vOFS='\\t' '{{for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i = "no_clan"}};1' \
    | sort -k1,1) > {params.multi_annot}

    # Reorder columns
    # the previous output is: pfam_name - orf - start - end - pfam_acc - pfam_clan
    # we want to get: orf - pfam_name - pfam_acc - pfam_clan - start - end
    awk -vOFS='\\t' '{{print $2,$1,$5,$6,$3,$4}}' {params.multi_annot} > {params.tmp} && mv {params.tmp} {params.multi_annot}

    # Multiple annotations on the same line, separated by “|” (the annotations were ordered first by alignment position)
    sort -k1,1 -k5,6g {params.multi_annot} | \
        awk -vFS='\\t' -vOFS='\\t' '{{print $1,$2,$3,$4}}' | sort -k1,1 | \
        awk -f {params.concat} | sed 's/ /\t/g' > {params.tmp}

    mv {params.tmp} {params.multi_annot}

    gzip -f {params.multi_annot}

    ## 2. Cluster annotations

    # The r script "clu_annot.r" distribute the Pfam annotation in the clusters,
    # creating two sets: "annotated_clusters" and "not_annotated_clusters"
    ./{params.cl_annotr} --pfam_annot {params.multi_annot}.gz \
                         --clusters {input.clu} \
                         --partial {input.partial} \
                         --output_annot {output.cl_annot} \
                         --output_noannot {output.cl_noannot}

    ## 3. Singleton annotations

    join -12 -21 <(sort -k2,2 {params.singl} ) \
    <(sort -k1,1 {params.multi_annot}) > {params.s_annot}

    """

run:
    shell("echo 'CLUSTER PFAM ANNOTATION DONE'")

shell:
  """
  set -x
  set -e

  if [ {params.module} == "creation" ]; then
      DB={params.clseqdb}
      PARTIAL={params.partial}
  else
      DB={params.new_clseqdb}
      PARTIAL={params.new_partial}
  fi


  # Step performed after the cluster validation, to remove:
  # 1. bad clusters (≥ 10% bad-aligned ORFs)
  # 2. shadow clusters (≥ 30% shadow ORFs)
  # 3. single rejected ORFs (shadow, spurious and bad-aligned)

  if [[ ! -s {params.toremove} ]]; then
      # 1st step: from the validation results we already have a table with the good clusters:
      # good_clusters.tsv (in the directory with the validation results)= {input.good}

      # 2nd step: we remove from the good clusters those with ≥ 30% shadows
      join -11 -21 -v2 <(awk '$5>=0.3 && !seen[$2]++{{print $2}}' {input.sp_sh} | sort -k1,1) \
      <( grep -v 'cl_name' {input.good} | awk '{{print $1}}' | sort -k1,1) > {params.good_noshadow}

      ## Retrieve the subdb with the ORFs of these clusters
      {params.mmseqs_bin} createsubdb <(awk '{{print $1}}' {params.good_noshadow}) ${{DB}} {params.tmpdb}

      # 3rd step: retrieve the single rejected/spurious/shadow ORFs and remove them from the remaining clusters
      ## retrieve the bad-aligned sequences in the set of good & no-shadow clusters
      join -11 -21 <(sort -k1,1 {input.cval_rej} ) <(sort -k1,1 {params.good_noshadow}) > {params.rej}

      ## Add the bad-aligned sequences to the spurious and shadows
      awk '$5<0.3 && $6!="FALSE"{{print $1}}' {input.sp_sh} > {params.toremove}
      awk '$5<0.3 && $7!="FALSE"{{print $1}}' {input.sp_sh} >> {params.toremove}
      awk '{{print $2}}' {params.rej} >> {params.toremove}
  fi

  if [[ -s {params.toremove} ]]; then
          # add cluster membership
          join -11 -21 <(sort {params.toremove}) \
          <(awk '{{print $3,$1}}' {params.clu_info} | sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
          > {params.toremove_cl}

          # remove the single orfs from the clusters with less than 10% bad-aligned ORFs and less than 30% shadows.
          if [[ ! -s {params.refdb} ]]; then
              {params.mpi_runner} {params.mmseqs_bin} apply {params.tmpdb} {params.refdb} --threads {threads} \
              -- {params.toremove_sh} {params.toremove_cl} 2>>{log.err} 1>>{log.out}
          fi
          # Create tables with new seqs and new clusters for some stats and checking the numbers
          join -11 -21 <(awk '{{print $1}}' {params.refdb}.index | sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
              <(awk '{{print $1,$3}}' {params.clu_info} \
              | sort -k1,1 --parallel={threads} -T {params.local_tmp}) > {params.tmp1}
          # Remove "bad" ORFs from the refined table
          join -12 -21 -v1 <(sort -k2,2 --parallel={threads} -T {params.local_tmp} {params.tmp1}) \
              <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.toremove}) > {params.tmp} && mv {params.tmp} {params.tmp1}
  else
    mv {params.good_noshadow} {params.tmp1}
    mv {params.tmpdb} {params.refdb}
    mv {params.tmpdb}.index {params.refdb}.index
    mv {params.tmpdb}.dbtype {params.refdb}.dbtype
  fi
  # From the refined clusters select the annotated and the not annotated for the following classification steps

  # annotated (check those left with no-annotated sequences): join with file with all annotated clusters (using the orfs)
  join -11 -21 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.tmp1}) \
    <(awk '{{print $3,$4}}' {params.annot} \
    |  sort -k1,1 --parallel={threads} -T {params.local_tmp}) > {output.ref_annot}

  # Reorder columns to have cl_name - orf - pfam_name
  awk -vOFS='\\t' '{{print $2,$1,$3}}' {output.ref_annot} > {params.tmp} && mv {params.tmp} {output.ref_annot}

  # Find in the refined annotated clusters, those left with no annotated ORFs
  sort -k1,1 --parallel={threads} -T {params.local_tmp} {output.ref_annot} \
    | awk '!seen[$1,$3]++{{print $1,$3}}' \
    | awk 'BEGIN{{getline;id=$1;l1=$1;l2=$2;}}{{if($1 != id){{print l1,l2;l1=$1;l2=$2;}}else{{l2=l2"|"$2;}}id=$1;}}END{{print l1,l2;}}' \
    | grep -v "|" | awk '$2=="NA"{{print $1}}' > {params.new_noannot}

  if [[ -s {params.new_noannot} ]]; then

    # move the clusters left with no annotated member to the not annotated
    join -11 -21 -v1 <(awk '!seen[$1,$2,$3]++' {output.ref_annot} | sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
      <(sort {params.new_noannot}) > {params.tmp}

    join -11 -21 <(awk '!seen[$1,$2,$3]++{{print $1,$2}}' {output.ref_annot} |\
      sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
      <(sort {params.new_noannot}) > {output.ref_noannot}

    # not annotated
    join -12 -21 <(sort -k2,2 --parallel={threads} -T {params.local_tmp} {params.tmp1}) \
      <(awk -vFS='\\t' '$4=="noannot"{{print $1,$2}}' {input.good} |\
       sort -k1,1 --parallel={threads} -T {params.local_tmp}) >> {output.ref_noannot}

    mv {params.tmp} {output.ref_annot}
    sed -i 's/ /\t/g' {output.ref_annot}

  else

    # not annotated
    join -12 -21 <(sort -k2,2 --parallel={threads} -T {params.local_tmp} {params.tmp1}) \
      <(awk '$4=="noannot"{{print $1,$2}}' {input.good} |\
       sort -k1,1 --parallel={threads} -T {params.local_tmp}) > {output.ref_noannot}
  fi

  awk -vOFS='\\t' '{{print $1,$2}}' {output.ref_noannot} > {params.tmp} && mv {params.tmp} {output.ref_noannot}

  # Using the cluster ids retrieve the sub database for not annotated clusters
  {params.mmseqs_bin} createsubdb <(awk '!seen[$1]++{{print $1}}' {output.ref_noannot}) \
                                  {params.refdb} \
                                  {params.refdb_noannot} 2>>{log.err} 1>>{log.out}

  # Add ORF partial information to the refined cluster table

  # Colums partial: 1:orf|2:partial_info
  # Columns tmp3: 1:orf|2:partial_info|3:cl_name
  join -11 -21 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} ${{PARTIAL}}) \
      <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.tmp1}) > {params.tmp2}

  # Add cluster size and ORF length
  # Columns clu_info: 1:cl_name|2:old_repr|3:orf|4:orf_len|5:cl_size
  # Columns intermediate: 1:orf|2:partial_info|3:cl_name|4:old_repr|5:orf_len|6:cl_size
  # Columns tmp4: 1:cl_name|2:orf|3:partial|4:cl_size|5:orf_len
  join -11 -21 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.tmp2}) \
      <(awk '{{print $3,$2,$4,$5}}' {params.clu_info} \
      | sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
      | sort -k3,3 --parallel={threads} -T {params.local_tmp} \
      | awk -vOFS='\\t' '{{print $3,$1,$2,$6,$5}}' > {params.tmp3}

  # Add new representatives from validation
  # Columns tmp5: 1:cl_name|2:orf|3:partial|4:cl_size|5:orf_len|6:new_repr
  join -11 -21 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.tmp3}) \
      <(awk '{{print $1,$2}}' {input.val_res} | sort -k1,1 --parallel={threads} -T {params.local_tmp}) > {params.tmp4}

  # Reorder columns: 1:cl_name|2:new_repres|3:orf|4:cl_size|5:orf_length|6:partial
  awk -vOFS='\\t' '{{print $1,$6,$2,$4,$5,$3}}' {params.tmp4} > {output.ref_clu}

  # cleaning the results...
  rm -rf {params.tmp1} {params.tmp2} {params.tmp3} {params.tmp4}
  rm -rf {params.toremove} {params.rej}
  rm -rf {params.new_noannot} {params.good_noshadow}
  rm -rf {params.tmpdb} {params.tmpdb}.index {params.tmpdb}.dbtype

  if [ {params.module} == "update" ]; then
      rm {params.new_partial}
  fi

  """

run:
    shell("echo 'CLUSTER REFINEMENT DONE'")

shell:
    """
    set -e
    set -x

    ## 1. Pfam-clan info and multi-domain format for the ORF Pfam annotations
    # Download original dataset Pfam annotations
    if [[ ! -s {params.or_multi_annot} ]]; then
        wget https://ndownloader.figshare.com/files/23067146 -O {params.or_multi_annot}
    fi

    if [ ! -s {params.pfam_clan} ]; then
        echo "Dowloading Pfam-A clan information"
        wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam34.0/Pfam-A.clans.tsv.gz -O {params.pfam_clan}
    fi

    # Add Pfam clan info to the Pfam annotation results
    join -11 -21 <( awk '{{print $1,$3,$8,$9}}' {input.annot} | sort -k1,1) \
    <( zcat {params.pfam_clan} | awk -vFS='\\t' -vOFS='\\t' '{{print $4,$1,$2}}' \
    | awk -vFS='\\t' -vOFS='\\t' '{{for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i = "no_clan"}};1' \
    | sort -k1,1) > {params.multi_annot}

    if [[ {params.db_mode} == "memory" ]]; then
        rm {params.pfam_clan}
    fi

    # Reorder columns
    # the previous output is: pfam_name - orf - start - end - pfam_acc - pfam_clan
    # we want to get: orf - pfam_name - pfam_acc - pfam_clan - start - end
    awk -vOFS='\\t' '{{print $2,$1,$5,$6,$3,$4}}' {params.multi_annot} > {params.tmp} && mv {params.tmp} {params.multi_annot}

    # Multiple annotations on the same line, separated by “|” (the annotations were ordered first by alignment position)
    sort -k1,1 -k5,6g {params.multi_annot} | \
        awk -vFS='\\t' -vOFS='\\t' '{{print $1,$2,$3,$4}}' | sort -k1,1 | \
        awk -f {params.concat} | sed 's/ /\t/g' > {params.tmp}

    mv {params.tmp} {params.multi_annot}

    zcat {params.or_multi_annot} | sed 's/ /\\t/g' >> {params.multi_annot}

    # Gene completeness information combined
    # Download original dataset gene completion information
    if [[ ! -s {params.or_partial} ]]; then
        wget https://ndownloader.figshare.com/files/23067005 -O {params.or_partial}
    fi

    # Combine with new completion info
    join -11 -21 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {input.partial}) \
     <( awk '{{print $3}}' {input.clu} | sort -k1,1 --parallel={threads}) > {params.partial}
    join -11 -21 <(zcat {params.or_partial} | \
     sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
     <( awk '{{print $3}}' {input.clu} | sort -k1,1 --parallel={threads}) >> {params.partial} 2>>{log.err}

    sed -i 's/ /\\t/g' {params.partial}

    ## 2. Cluster annotations

    # The r script "clu_annot.r" distribute the Pfam annotation in the clusters,
    # creating two sets: "annotated_clusters" and "not_annotated_clusters"
    ./{params.cl_annotr} --pfam_annot {params.multi_annot} \
                         --clusters {input.clu} \
                         --partial {params.partial} \
                         --output_annot {output.cl_annot} \
                         --output_noannot {output.cl_noannot} 2>>{log.err}

    ## 3. Singleton annotations

    join -12 -21 <(sort -k2,2 {params.singl} ) \
    <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.multi_annot}) > {params.s_annot}

    gzip -f {params.multi_annot}

    """

run:
    shell("echo 'CLUSTER PFAM ANNOTATION DONE'")

shell:
  """
  set -e
  set -x

  #Retrieve old cluster representatives information
  # Not annotated clusters
  #
  join -11 -21 <(awk '{{print $1,$2}}' {input.cval} | sort -k1,1 --parallel={threads}) \
    <(awk '!seen[$1]++{{print $1,$2,"noannot",$4}}' {params.cl_noannot} | sort -k1,1 --parallel={threads}) > {params.val_annot}
  # Annotated clusters
  join -11 -21 <(awk '{{print $1,$2}}' {input.cval} | sort -k1,1 --parallel={threads}) \
    <(awk '!seen[$1]++{{print $1,$2,"annot",$4}}' {params.cl_annot} | sort -k1,1 --parallel={threads} ) >> {params.val_annot}

  awk -vOFS='\\t' '{{print $1,$2,$3,$5,$4}}'  {params.val_annot} \
      > {params.tmp} && mv {params.tmp} {params.val_annot}

  # Combine with functional validation results
   ./{params.val_res} --fval_res {input.fval} \
                      --cval_res {input.cval} \
                      --val_annot {params.val_annot} \
                      --val_res {output.val_res} \
                      --val_stats {params.val_stats} \
                      --good {output.good} \
                      --plots {params.val_plots} 2>{log.err} 1>{log.out}

  rm -rf {params.val_annot} {params.tmp}
  """

run:
    shell("echo 'CLUSTER VALIDATION DONE'")

shell:
    """
    Rscript --vanilla {params.report_maker} --basedir {params.basedir} \
                                            --outdir  {params.outdir} \
                                            --stage {params.stage} \
                                            --name {params.name_data} \
                                            --input {params.input_data} \
                                            --wf_report {params.wf_report} \
                                            --output {output.report} 2>{log.err} 1>{log.out}

    """

shell:
    """
    set -x
    set -e

    if [[ {params.stage} = "contigs" ]]; then

        {params.prodigal_bin} -i {input.contigs} -a {output.fa} -m -p {params.prodigal_mode} -f gff  -o {params.gff_output} -q 2>{log.err} 1>{log.out}

        awk -f {params.rename_orfs} {output.fa} > {params.tmp} && mv {params.tmp} {output.fa}

        awk -f {params.partial_info} {params.gff_output} > {output.partial}

    elif [[ {params.stage} = "genes" ]]; then

        ln -sf {input.contigs} {output.fa}

        ln -sf {params.data_partial} {output.partial}

    elif [[ {params.stage} = "anvio_genes" ]]; then

        ln -sf {input.contigs} {output.fa}

        awk -vOFS="\\t" 'NR>1{{if($6==0) print $1,"00"; else print $1,"11";}}' {params.data_partial} > {output.partial}

    fi

    """

run:
    shell("echo 'GP DONE'")

shell:
    """

    set -x
    set -e

    DIR=$(dirname {output.iclu_hmm})

    mkdir -p ${{DIR}}

    if [[ ! -s {params.or_clu_orig} ]]; then
        wget https://ndownloader.figshare.com/files/23066966 -O {params.or_clu_orig}
    fi
    # Summary table with cluster db origin (original/shared/new)
    join -11 -21 <(zcat {params.or_clu_orig} | awk '{{print $1,$2}}' | sort -k1,1 --parallel={threads} ) \
        <(awk '{{print $1,$3}}' {params.original} | sort -k1,1) > {params.clu_origin}
    join -11 -21 <(zcat {params.or_clu_orig} | awk '{{print $1,$2}}' |  sort -k1,1 --parallel={threads} ) \
        <(awk '{{print $1,$3}}' {params.shared} | sort -k1,1) > {params.clu_origin}.temp
    awk -vN={params.data_name} '{{print $1,$2"_"N,$3}}' {params.clu_origin}.temp >> {params.clu_origin}
    awk -vN={params.data_name} '{{print $1,N,$3}}' {params.new} >> {params.clu_origin}
    rm {params.clu_origin}.temp
    sed -i 's/ /\t/g' {params.clu_origin}
    gzip {params.clu_origin}
    # Clean the mmseqs_clustering/ folder
    mv {params.all} {params.final}

    rm {params.original} {params.shared} {params.new} {params.new_cluseqdb}*

    # All gene headers and partiality information
    cat {params.partial} <(zcat {params.or_partial}) | gzip > {params.ipartial}

    # Spurious and shadow genes information:
    gzip -c {params.sp_sh} > {params.isp_sh}

    # All gene Pfam annotations:
    cp {params.multi_annot} {params.imulti_annot}

    # All cluster category annotation files
    ODIR=$(dirname {params.or_clu_cat})
    NDIR=$(dirname {input.clu_cat})

    # Download original dataset category annotations
    if [[ ! -s ${{ODIR}}/K_annotations.tsv.gz ]]; then
        wget https://ndownloader.figshare.com/files/23063648 -O ${{ODIR}}/K_annotations.tsv.gz
        wget https://ndownloader.figshare.com/files/23067074 -O ${{ODIR}}/KWP_annotations.tsv.gz
        wget https://ndownloader.figshare.com/files/23067080 -O ${{ODIR}}/GU_annotations.tsv.gz
    fi
    # Combine with new ones
    cat <(zcat ${{ODIR}}/K_annotations.tsv.gz) ${{NDIR}}/K_annotations.tsv | gzip > {params.idir}/K_annotations.tsv.gz
    cat <(zcat ${{ODIR}}/KWP_annotations.tsv.gz) ${{NDIR}}/KWP_annotations.tsv | gzip > {params.idir}/KWP_annotations.tsv.gz
    cat <(zcat ${{ODIR}}/GU_annotations.tsv.gz) ${{NDIR}}/GU_annotations.tsv | gzip > {params.idir}/GU_annotations.tsv.gz
    # rm ${{ODIR}}/*_annotations.tsv.gz

    # Integrated set of cluster categories
    # Download original gene cluster catgeory info
    if [[ ! -s {params.or_clu_cat} ]]; then
        wget https://ndownloader.figshare.com/files/23067140 -O {params.or_clu_cat}
    fi

    # and the cluster genes
    # Download original gene cluster catgeory info
    if [[ ! -s {params.or_clu_gene} ]]; then
        wget https://ndownloader.figshare.com/files/24121865 -O {params.or_clu_gene}
    fi

    if [[ {params.eval_shared} == "true" ]]; then
        # GC-categories
        join -11 -21 -v1 <(zcat {params.or_clu_cat} |\
         sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
         <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {input.clu_cat}) > {params.tmpl2}
        sed -i 's/ /\\t/g' {params.tmpl2}
        cat {params.tmpl2} {input.clu_cat} | gzip > {output.iclu_cat}

         # GC-genes-categories: get the new genes in the good clusters (shared)
         # remove shared GCs from the old table
         join -11 -21 -v1 <(zcat {params.or_clu_gene} |\
          sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
          <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {input.clu_cat}) > {params.tmpl1}
         sed -i 's/ /\\t/g' {params.tmpl1}
         cat {params.tmpl1} <(zcat {params.clu_gene}) | gzip > {params.iclu_gene}

         rm {params.tmpl1}
    else
        # GC-categories
        cat <(zcat {params.or_clu_cat}) {input.clu_cat} | gzip > {output.iclu_cat}

        # GC-genes-categories: get the new genes in the good clusters (shared)
        ### double join, first cluDB_info, then orf_seqs.txt to get the new ones only
        join -11 -21 <(zcat {params.or_clu_cat} | sort -k1,1) \
            <(awk '{{print $1,$3}}' {params.clu_info} \
            | sort -k1,1 --parallel={threads} -T {params.local_tmp} ) > {params.tmpl1}
        join -13 -21 <(sort -k3,3 --parallel={threads} -T {params.local_tmp} {params.tmpl1}) \
            <(sort -k1,1 <(awk '{{print $1}}' {params.partial})) > {params.tmpl2}
        # add to original clu genes and new clu genes
        cat <(awk -vOFS='\\t' '{{print $2,$3,$1}}' {params.tmpl2}) \
            <(zcat {params.or_clu_gene}) > {params.tmpl1}
        cat {params.tmpl1} <(zcat {params.clu_gene}) | gzip > {params.iclu_gene}
        rm {params.tmpl1} {params.tmpl2}
    fi

    if [ {params.singl} == "true" ]; then
        if [[ ! -s {params.or_singl} ]]; then
            wget https://figshare.com/ndownloader/files/31012435 -O {params.or_singl}
        fi
        join -11 -21 <(zcat {params.clu_origin}.gz | awk '$3==1{{print $1}}' |\
            sort -k1,1 --parallel={threads} -T {params.local_tmp} ) \
            <(zcat {params.or_singl} | sort -k1,1 --parallel={threads} -T {params.local_tmp}) > {params.tmpl1}

        cat <(sed 's/ /\t/g' {params.tmpl1}) \
         <( awk -vOFS="\\t" '{{print $2,$3,$1}}' {params.s_categ}) > {params.singl_cl_gene_categ}

        gzip {params.singl_cl_gene_categ}
     fi

    # Integrated set of cluster communities
    # to avoid having overlapping communities names, the dataset origin is appended to the name
    if [[ ! -s {params.or_clu_com} ]]; then
        wget https://ndownloader.figshare.com/files/23067134 -O {params.or_clu_com}
    fi

    ./{params.integr_comm} --name {params.data_name} \
                           --comm {input.clu_com} \
                           --ocomm {params.or_clu_com} \
                           --icomm {output.iclu_com} \
                           --shared {params.eval_shared} 2>>{log.err} 1>>{log.out}

    gzip {output.iclu_com}

    # Integrated cluster summary information
    if [[ ! -s {params.or_clu_stats} ]]; then
        wget https://figshare.com/ndownloader/files/31003162 -O {params.or_clu_stats}
    fi
    if [[ {params.eval_shared} == "true" ]]; then
        join -11 -21 -v1 <(zcat {params.or_clu_stats} | sort -k1,1 --parallel={threads} -T {params.local_tmp}) \
            <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {input.clu_stats}) > {params.tmpl1}
        sed -i 's/ /\\t/g' {params.tmpl1}
        cat {input.clu_stats} {params.tmpl1} | gzip > {output.iclu_stats}
    else
        cat {input.clu_stats} <(zcat {params.or_clu_stats} | awk -vOFS='\\t' 'NR>1{{print $0}}') | gzip > {output.iclu_stats}
    fi

    # Integrated set of high quality (HQ) clusters
    if [[ ! -s {params.or_hq_clu} ]]; then
        wget https://ndownloader.figshare.com/files/23067137 -O {params.or_hq_clu}
    fi

    if [[ {params.eval_shared} == "true" ]]; then
        join -11 -21 -v1 <(zcat {params.or_hq_clu} | sort -k1,1 ) \
            <(sort -k1,1 {input.hq_clu}) > {params.tmpl1}
        sed -i 's/ /\\t/g' {params.tmpl1}
        cat {input.hq_clu} {params.tmpl1} | gzip > {output.ihq_clu}
    else
        cat {input.hq_clu} <(zcat {params.or_hq_clu} | awk -vOFS='\\t' 'NR>1{{print $0}}' ) | gzip > {output.ihq_clu}
    fi

    # New integarted cluster HMMs DB (for MMseqs profile searches)
    # Download and uncompress the existing GC profiles
    if [[ ! -s {params.or_clu_hhm} ]]; then
        wget https://figshare.com/ndownloader/files/30998305 -O {params.or_profiles}.tar.gz
        tar -C {params.or_dir} -xzvf {params.or_profiles}.tar.gz
        rm {params.or_profiles}.tar.gz
    fi

    if [[ {params.eval_shared} == "true" ]]; then
        {params.mmseqs_bin} createsubdb <(awk '{{print $1}}' {params.tmpl2}) \
                                        {params.or_clu_hhm} {params.or_clu_hhm}.left
        {params.mmseqs_bin} concatdbs {input.clu_hhm} {params.or_clu_hhm}.left {params.iclu_hhm} \
                                       --threads 1 --preserve-keys 2>{log.err}
        # Create a comprehensive profile cluster DB in MMseqs format to perform profile searches
        {params.mmseqs_bin} convertprofiledb {params.iclu_hhm} {output.iclu_hmm} \
                                             --threads {threads} 2>{log.err}
        rm {params.or_clu_hhm}.left*
    else
        {params.mmseqs_bin} concatdbs {input.clu_hhm} {params.or_clu_hhm} {params.iclu_hhm} \
                                      --threads 1 --preserve-keys 2>{log.err}
        # Create a comprehensive profile cluster DB in MMseqs format to perform profile searches
        {params.mmseqs_bin} convertprofiledb {params.iclu_hhm} {output.iclu_hmm} \
                                             --threads {threads} 2>{log.err}
    fi
    """

run:
    shell("echo 'INTEGRATED DB DONE'")

shell:
    """
    set -e
    set -x

    # Create sequences database for the clustering
    {params.mmseqs_bin} createseqfiledb {params.seqdb} {params.cludb} {params.cluseqdb} --threads {threads} 2>{log.err}

    # To convert this cluster results tab separated file in wide format (repres member member member ..)
    awk -f {params.awk_wide} {input.clu} > {params.wide}

    # Cluster name rep and number of ORFs (size)
    awk -vFS='\\t' -vOFS='\\t' '{{print NR,$1,NF-1}}' {params.wide} > {params.info1}

    if [[ -s {params.cludb}.0 ]]; then
        {params.concat_clu} {params.cludb} {threads}
    fi

    if [[ -s {params.cluseqdb}.0 ]]; then
        {params.concat_clu} {params.cluseqdb} {threads}
    fi

    # Cluster naming
    # The official cluster names are going to be based on the line number of the wide formatted file
    # We are also going to produce a correspondence file to access the clusters in the MMseqs2 indices
    # Retrieving MMseqs-entry-name and sequence header from the clu_seqDB

    {params.mpi_runner} {params.mmseqs_bin} apply {params.cluseqdb} {params.namedb} -- {params.naming} --threads {threads}

    # Join the sequences with the representatives: this way we combine mmseqs-name with cluster-name (row number)
    join -12 -22 <(sed -e 's/\\x0//g' {params.namedb} | sort -k2,2 --parallel={threads} -T {params.local_tmp}) \
    <(sort -k2,2 --parallel={threads} -T {params.local_tmp} {params.info1} ) > {output.index}.tmp

    awk -vOFS='\\t' '{{print $2,$3}}' {output.index}.tmp > {output.index}

    # Rename the cluster sequence DB with the cluster names that are goiung to be used in the next steps
    {params.mmseqs_bin} renamedbkeys {output.index} {params.cluseqdb} {params.cluseqdb}.new

    mv {params.cluseqdb}.new {params.cluseqdb}
    mv {params.cluseqdb}.new.index {params.cluseqdb}.index
    mv {params.cluseqdb}.new.dbtype {params.cluseqdb}.dbtype

    # Join with the long format cluster file
    join -11 -22 <(sort --parallel={threads} -T {params.local_tmp} -k1,1 {input.clu} ) \
        <(sort --parallel={threads} -T {params.local_tmp} -k2,2 {params.info1} ) > {params.tmp}

    # Add length info
    sed -e 's/\\x0//g' {params.cluseqdb} | {params.seqtk_bin} comp | cut -f1,2 > {params.length}

    join -11 -22 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.length}) \
    <(sort -k2,2 --parallel={threads} -T {params.local_tmp} {params.tmp}) > {output.clu_info}

    # Reorder fields (cl_name rep orf length size)
    sort -k4,4n --parallel={threads} -T {params.local_tmp} {output.clu_info} | \
     awk -vOFS='\\t' '{{print $4,$3,$1,$2,$5}}' > {params.tmp}

    cp {params.tmp} {output.clu_info}

    # Subset the clusters based on their size

    # Clusters with more than 1 member
    awk -vOFS='\\t' '$3>=2{{print $1,$2}}' {params.info1} > {params.tmp1}
    join -12 -21 <(sort -k2,2 --parallel={threads} -T {params.local_tmp} {params.tmp1}) \
     <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {input.clu} ) > {params.tmp}

    awk -vOFS='\\t' '{{print $2,$1,$3}}' {params.tmp} > {output.clusters}

    # Singletons
    awk -vOFS='\\t' '$3=="1"{{print $1,$2}}' {params.info1} > {output.singl}

    rm {params.tmp} {params.tmp1} {params.namedb} {params.namedb}.index {params.namedb}.dbtype {output.index}.tmp

    rm {params.wide} {params.length}

    """

run:
    shell("echo 'CLUSTERING RESULTS DONE'")

shell:
    """
    set -x
    set -e

    export OMPI_MCA_btl=^openib
    export OMP_NUM_THREADS={threads}
    export OMP_PROC_BIND=FALSE

    {params.mmseqs_bin} createdb {input.orfs} {params.seqdb} 2>{log.err} 1>{log.out}

      {params.mmseqs_bin} cluster \
      {params.seqdb} \
      {params.cludb} \
      {params.mmseqs_tmp} \
      --local-tmp {params.mmseqs_local_tmp} \
      --threads {threads} \
      -c {params.mmseqs_cov} \
      --cov-mode {params.mmseqs_cov_mode} \
      --min-seq-id {params.mmseqs_id} \
      -s {params.mmseqs_ens} \
      --mpi-runner "{params.mmseqs_mpi_runner}" 2>>{log.err} 1>>{log.out}

    {params.mmseqs_bin} createtsv {params.seqdb} {params.seqdb} {params.cludb} {output.clu}
    """

run:
    shell("echo 'MMSEQS2 CLUSTERING DONE'")

shell:
    """
    set -x
    set -e

    export OMPI_MCA_btl=^openib
    export OMP_NUM_THREADS={threads}
    export OMP_PROC_BIND=FALSE

    # Create directory for original data to be integrated
    mkdir -p {params.or_dir}

    # Download original mmseqs clustering that you want to use for the integration
    if [[ ! -s {params.or_cludb}.index ]]; then
        wget https://ndownloader.figshare.com/files/23066651 -O {params.or_dir}/mmseqs_clustering.tar.gz
        tar -C {params.or_dir} -zxvf {params.or_dir}/mmseqs_clustering.tar.gz
        rm {params.or_dir}/mmseqs_clustering.tar.gz
    fi
    {params.mmseqs_bin} createdb {input.new_seqs} {params.new_seqdb} 2>{log.err} 1>{log.out}

    # Create symbolic link between original cluDB and seqDB (required by mmeseqs to run the update)

    ln -sf {params.or_seqdb}_h {params.or_cludb}_h
    ln -sf {params.or_seqdb}_h.index {params.or_cludb}_h.index
    ln -sf {params.or_seqdb}_h.dbtype {params.or_cludb}_h.dbtype
    ln -sf {params.or_seqdb}.lookup {params.or_cludb}.lookup

    # Concat new and old seq DBs
    {params.mmseqs_bin} concatdbs {params.or_seqdb} {params.new_seqdb} \
        {params.conc_seqdb} \
        --threads 1  2>>{log.err} 1>>{log.out}
    {params.mmseqs_bin} concatdbs {params.or_seqdb}_h {params.new_seqdb}_h \
        {params.conc_seqdb}_h \
        --threads 1  2>>{log.err} 1>>{log.out}

    # Update the clusterDB:
    {params.mmseqs_bin} clusterupdate \
      {params.or_seqdb} \
      {params.conc_seqdb} \
      {params.or_cludb} \
      {params.updt_seqdb} \
      {params.updt_cludb} \
      {params.mmseqs_tmp} \
      --local-tmp {params.mmseqs_local_tmp} \
      --mpi-runner "{params.mmseqs_mpi_runner}" \
      --threads {threads} \
      -c {params.mmseqs_cov} \
      --cov-mode {params.mmseqs_cov_mode} \
      --min-seq-id {params.mmseqs_id} \
      -s {params.mmseqs_ens} 2>>{log.err} 1>>{log.out}

    {params.mmseqs_bin} createtsv \
      {params.updt_seqdb} \
      {params.updt_seqdb} \
      {params.updt_cludb} \
      {output.updt_clu} \
      --threads {threads} 2>>{log.err} 1>>{log.out}

     #mv {params.conc_seqdb} {params.updt_seqdb}
     #mv {params.conc_seqdb}.dbtype {params.updt_seqdb}.dbtype
     #mv {params.conc_seqdb}_h {params.updt_seqdb}_h
     #mv {params.conc_seqdb}_h.dbtype {params.updt_seqdb}_h.dbtype

     #rm {params.conc_seqdb}*

    """

run:
    shell("echo 'MMSEQS2 CLUSTERING DONE'")

shell:
    """
    set -e
    set -x

    mkdir -p {params.mmseqs_tmp}

    # Create sequences database for the clustering
    # This DB can be accessed with ffindex, to extractt separated fasta files for each cluster and perform operations on them
    # Ex: ffindex_apply_mpi cluseqdb cluseqdb.index -- your_program/script
    if [[ ! -s {params.cluseqdb}.index ]]; then

        {params.mmseqs_bin} createseqfiledb {params.seqdb} {params.cludb} {params.cluseqdb} --threads {threads} 2>{log.err}

        #Check if sequences contain "*", and if yes, clean them using the following command
        #{params.mpi_runner} {params.mmseqs_bin} apply \
        #    {params.cluseqdb} \
        #    {params.cluseqdb}_clean  \
        #    --threads {threads} \
        #    -- {params.clean_seqs} 2>>{log.err}
    fi

    if [[ -s {params.cludb}.0 ]]; then
        {params.concat_clu} {params.cludb} {threads}
    fi

    if [[ -s {params.cluseqdb}.0 ]]; then
        {params.concat_clu} {params.cluseqdb} {threads}
    fi

    # To convert this cluster results tab separated file in wide format (repres member member member ..)
    awk -f {params.awk_wide} <(awk '!seen[$2]++' {input.clu}) > {params.wide} 2>>{log.err}

    # Cluster name rep and number of ORFs (size)
    awk -vFS='\\t' -vOFS='\\t' '{{print NR,$1,NF-1}}' {params.wide} > {params.info1}


    # Join with the original clustering using the representative sequences
    join -12 -22 <(sort -k2,2 {params.orig_info1}) \
      <(sort -k2,2 {params.info1}) > {params.rename_or} 2>>{log.err}

    join -12 -22 -v2 <(sort -k2,2 {params.orig_info1}) \
      <(sort -k2,2 {params.info1}) > {params.rename_new} 2>>{log.err}

    # Store the information about changes in the original Clusters
    sort -k2,2n {params.rename_or} | awk -vOFS='\\t' '{{print $2,$1,$3,$5}}' > {params.or_new_info}
    # Order and re-join old and new cluster ids/NAMES (the new ids will start from the last of the original clustering)
    OR=$(wc -l {params.or_new_info} | cut -d ' ' -f1)

    awk -vOR=$OR -vOFS='\\t' '!seen[$0]++{{print NR+OR,$1,$3}}' {params.rename_new} > {params.new_info1} 2>>{log.err}

    cat <(awk -vOFS='\\t' '!seen[$0]++{{print $1,$2,$4}}' {params.or_new_info} ) \
      {params.new_info1} > {params.info1}

    # Cluster naming
    # The official cluster names are going to be based on the line number of the wide formatted file
    # We are also going to produce a correspondence file to access the clusters in the MMseqs2 indices
    # Retrieving MMseqs-entry-name and sequence header from the clu_seqDB
    if [[ ! -s {params.index} ]]; then
        {params.mpi_runner} {params.mmseqs_bin} apply {params.cluseqdb} {params.namedb} --threads {threads} -- {params.naming}

        # Join the sequences with the representatives: this way we combine mmseqs-name with cluster-name (row number)
        join -12 -22 <(sed -e 's/\\x0//g' {params.namedb} | sort -k2,2 --parallel={threads} -T {params.local_tmp}) \
            <(sort -k2,2 --parallel={threads} -T {params.mmseqs_tmp} {params.info1} ) > {params.index}
        awk -vOFS='\\t' '!seen[$0]++{{print $2,$3}}' {params.index} > {params.tmp} && mv {params.tmp} {params.index}
    fi

    # Rename the cluster sequence DB with the cluster names that are goiung to be used in the next steps
    {params.mmseqs_bin} renamedbkeys {params.index} {params.cluseqdb} {params.cluseqdb}.new

    mv {params.cluseqdb}.new {params.cluseqdb}
    mv {params.cluseqdb}.new.index {params.cluseqdb}.index
    mv {params.cluseqdb}.new.dbtype {params.cluseqdb}.dbtype

    # Join with the long format cluster file
    join -11 -22 <(sort --parallel={threads} -T {params.local_tmp} -k1,1 {input.clu}) \
        <(sort --parallel={threads} -T {params.local_tmp} -k2,2 {params.info1} ) > {params.tmp1}

    # Add length info
    sed -e 's/\\x0//g' {params.cluseqdb} | {params.seqtk_bin} comp | cut -f1,2 > {params.length}

    join -11 -22 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.length}) \
     <(sort -k2,2 --parallel={threads} -T {params.mmseqs_tmp} {params.tmp1}) > {params.tmp}

    # Reorder fields (cl_name rep orf length size)
    sort -k4,4n --parallel={threads} -T {params.local_tmp} {params.tmp} | \
     awk -vOFS='\\t' '{{print $4,$3,$1,$2,$5}}' > {params.clu_info}

    ## Retrieve the different cluster sets:
    # 1. Only-original clusters: cluDB_original_name_rep_size.tsv

    awk -vOFS='\\t' '$3==$4{{print $1,$2,$3}}' {params.or_new_info} | awk '!seen[$0]++' > {params.original}

    # 2. Shared original-new clusters: cluDB_shared_name_rep_size.tsv

    awk -vOFS='\\t' '$3!=$4{{print $1,$2,$4}}' {params.or_new_info} | awk '!seen[$0]++' > {params.shared}

    # 3. Only-new clusters: cluDB_new_name_rep_size.tsv
    # we already have this set: {params.new_info1}

    ## We are going to re-validate and classified the clusters in new_cluDB and

    if [[ {params.eval_shared} == "true" ]]; then

        if [[ {params.shared_type} == "all" ]]; then
            # Evaluate all the shared GCs (GCs that got integrated with new genes)
            cat {params.new_info1} \
                <(awk -vFS='\\t' -vOFS='\\t' '!seen[$0]++{{print $1,$2,$3}}' {params.shared} ) > {params.info2}
        fi

        if [[ {params.shared_type} == "discarded" ]]; then
            # Add the previously discarded shared GCs
            # download if not there the original cluster - category table
            if [[ ! -s {params.or_clu_cat} ]]; then
                wget https://ndownloader.figshare.com/files/23067140 -O {params.or_clu_cat}
            fi
            # anti join shared GCs with the original cluster-category to get the discarded GCs
            join -11 -21 -v1 <(sort -k1,1 --parallel={threads} -T {params.local_tmp} {params.shared}) \
                <(zcat {params.or_clu_cat} | awk '{{print $1}}' |\
                 sort -k1,1 --parallel={threads} -T {params.local_tmp}) > {params.shared}.temp

            cat {params.new_info1} \
                <(awk -vFS='\\t' -vOFS='\\t' '!seen[$0]++{{print $1,$2,$3}}' {params.shared}.temp ) > {params.info2}
            rm {params.shared}.temp
        fi

        if [[ {params.shared_type} == "increase" ]]; then
            # Or the shared GCs with more than % new genes
            # Ex: increase of 30% in number of sequences
            awk '{{print $1,$2,$4,($4-$3)/$3}}' {params.or_new_info} > {params.or_new_info}.temp
            cat {params.new_info1} \
                <(awk -vFS='\\t' -vOFS='\\t' '!seen[$0]++ && $4>=0.3{{print $1,$2,$3}}' {params.or_new_info}.temp ) > {params.info2}
            rm {params.or_new_info}.temp
        fi
    else
        # Evaluate also the previously singletons now integrated in GCs with new genes
        cat {params.new_info1} \
            <(awk -vFS='\\t' -vOFS='\\t' '$3==1 && $4>1 && !seen[$0]++{{print $1,$2,$4}}' {params.or_new_info} ) > {params.info2}
    fi

    # Subset the cluseqDB:
    if [[ ! -s {params.new_cluseqdb} ]]; then
        {params.mmseqs_bin} createsubdb <(awk '{{print $1}}' {params.info2}) \
                                        {params.cluseqdb} \
                                        {params.new_cluseqdb} 2>>{log.err}
    fi

    # Subset the clusters based on their size
    # Clusters with more than 1 member
    awk -vOFS='\\t' '$3>=2{{print $1,$2}}' {params.info2} > {params.tmp1}
    join -12 -21 <(sort -k2,2 {params.tmp1}) \
    <(sort -k1,1 --parallel={threads} -T {params.local_tmp} <(awk '!seen[$2]++' {input.clu} )) > {params.tmp}

    awk -vOFS='\\t' '!seen[$0]++{{print $2,$1,$3}}' {params.tmp} > {output.clusters}

    # Singletons
    awk -vOFS='\\t' '$3=="1"{{print $1,$2}}' {params.info2} | awk '!seen[$0]++' > {output.singl}

    rm {params.tmp} {params.tmp1} {params.namedb} {params.namedb}.index {params.namedb}.dbtype
    rm {params.wide} {params.length} {params.rename_or} {params.rename_new}

    """

run:
    shell("echo 'CLUSTERING RESULTS DONE'")

shell:
    """

    set -x
    set -e

    DB=$(basename {params.orig_db})

    if [[ ${{DB}} == "agnostosDB" ]]; then
        # Download agnostosDB ecological analysis data
        wget https://ndownloader.figshare.com/files/23066879 -O ecological.tar.gz
        tar -xzvf ecological.tar.gz

        # Download agnostosDB phylogentic analysis data
        wget https://ndownloader.figshare.com/files/23066864 -O phylogenetic.tar.gz
        tar -xzvf phylogenetic.tar.gz

        # Download agnostosDB experimental analysis data
        wget https://ndownloader.figshare.com/files/23066864 -O phylogenetic.tar.gz
        tar -xzvf experimental.tar.gz
    fi

    # Run R script to retrieve a general table containing a summary of the integration
    # add contextual data if the original DB is the agnostosDB
    awk -vOFS='\\t' '{{split($1,a,"_\\\+|_-"); print a[1],$1}}' {input.genes} > {params.contig}

    ./{params.parser} --clu_or {params.clu_origin} \
                      --contig {params.contig} \
                      --cat {input.cat} \
                      --clu_info {params.clu_info} \
                      --name {params.new_data_name} \
                      --comm {params.comm} \
                      --hq_clu {params.hq_clu} \
                      --k_annot {params.k_annot} \
                      --kwp_annot {params.kwp_annot} \
                      --gu_annot {params.gu_annot} \
                      --orig_db {params.orig_db} \
                      --is_singl {params.singl} \
                      --s_categ {params.singl_cat} \
                      --anvio {params.data_stage} \
                      --threads {threads} 2>{log.err} 1>{log.out}

    """

run:
    shell("echo 'Parsing of results DONE'")

run:
    shell("echo 'ANNOTATION DONE'")

shell:
    """
    set -x
    set -e

    export OMPI_MCA_btl=^openib
    export OMP_NUM_THREADS={threads}
    export OMP_PROC_BIND=FALSE

    # 1. Detection of spurious ORFs

    if [ ! -s {params.antifamdb} ]; then
        echo "Dowloading AntiFam database"
        DB=$(dirname {params.antifamdb})
        wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/AntiFam/current/Antifam.tar.gz -P ${{DB}}
        tar xvfz ${{DB}}/Antifam.tar.gz --directory ${{DB}}
        rm ${{DB}}/Antifam.tar.gz
        {params.hmmpress} {params.antifamdb}
    fi

    NANTI=$(grep -c '^NAME' {params.antifamdb})
    NSEQS=$(grep -c '^>' {input.fasta})
    N=$(($NSEQS * $NANTI))

    # Run hmmsearch (MPI mode)
    {params.mpi_runner} {params.hmmer_bin} --mpi --cut_ga -Z "${{N}}" --domtblout {params.hmmout} -o {params.hmmlog} {params.antifamdb} {input.fasta} 2>{log.err} 1>{log.out}

    if [[ {params.db_mode} == "memory" ]]; then
        DB=$(dirname {params.antifamdb})
        rm -rf ${{DB}}/AntiFam*
    fi

    # Parse the results
    grep -v '^#' {params.hmmout} > {params.spur}.tmp || true > {params.spur}.tmp 2>>{log.err}

    rm {params.hmmout} {params.hmmlog}

    if [[ -s {params.spur}.tmp ]]; then
        awk '$13<=1e-05 && !seen[$1]++{{print $1}}' {params.spur}.tmp > {params.spur} 2>>{log.err}
        rm {params.spur}.tmp
    else
        mv {params.spur}.tmp {params.spur}
    fi

    # 2. Detection of shadow ORFs
    ./{params.shadowr} --orfs {params.partial} \
                       --shadows {params.all_shad} \
                       --threads {threads} 2>>{log.err}

    ## 2.1 Parsing of results adding cluster information
    ## Add cluster info on the first column of ORFs
    join -12 -23 <(sort --parallel={threads} -k2,2 {params.all_shad}) \
            <(sort --parallel={threads} -T {params.local_tmp} -k3,3 {params.clu_info} ) > {params.tmp1}

    ## reorder columns...
    awk -vOFS='\\t' '{{print $6,$7,$8,$9,$1,$11,$12,$13,$3,$4,$5,$2,$10}}' {params.tmp1} > {params.tmp2} && mv {params.tmp2} {params.tmp1}

    ## Add cluster info on the second column of ORFs
    join -11 -23 <(sort --parallel={threads} -k1,1 {params.tmp1}) \
            <(sort --parallel={threads} -T {params.local_tmp} -k3,3 {params.clu_info} ) > {params.tmp2}

    ## reorder columns...
    awk -vOFS='\\t' '{{print $1,$14,$15,$16,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13}}' {params.tmp2} > {params.tmp1}

    ### 2.2 Filter to get only the real shadow ORFs
    ### print ORFs and "TRUE" or "BOTH" if they fall in the definition of shadows
    awk -vFS='\\t' -vOFS='\\t' '{{if($4>$11) print $8,"TRUE"; else if($11>$4) print $1,"TRUE"; else if($4==$11 && $3<$10) print1; else if($4==$11 && $3>$10) print $1,"TRUE";}}' {params.tmp1} > {params.tmp2}
    awk -vFS='\\t' -vOFS='\\t' '{{if($4==$11 && $3==$10) print $1,"BOTH";}}' {params.tmp1} >> {params.tmp2}
    awk -vFS='\\t' -vOFS='\\t' '{{if($4==$11 && $3==$10) print $8,"BOTH";}}' {params.tmp1} >> {params.tmp2}

    ### re-join with the cluster information table
    join -13 -21 -a1 -a2 <(sort --parallel={threads} -T {params.local_tmp} -k3,3 {params.clu_info} ) \
      <(sort --parallel={threads} -k1,1 <(awk '!seen[$0]++' {params.tmp2}) ) > {params.only_shad}

    ### Add a "FALSE" to the non-shadow ORFs
    awk -vOFS='\\t' '!seen[$0]++{{if($6=="TRUE" || $6=="BOTH") print $2,$1,$3,$4,$5,$6; else print $2,$1,$3,$4,$5,"FALSE";}}' {params.only_shad} > {params.tmp1}
    mv {params.tmp1} {params.only_shad}

    # 3. Combine spurious and shadow information
    join -12 -21 -a1 <(sort --parallel={threads} -k2,2 {params.only_shad}) \
        <(awk '{{print $1,"TRUE"}}' {params.spur} | sort --parallel={threads} -k1,1) > {params.tmp1}

    awk '{{if($7=="TRUE") print $1,$4,$2,$3,$5,$6,$7; else print $1,$4,$2,$3,$5,$6,"FALSE";}}' {params.tmp1} | \
        awk '{{print $1,$2,$3,$5,$6,$7}}' > {output.sp_sh}.tmp

    # calculate proportion of shadows per cluster
    # grep 'TRUE\|BOTH' {params.only_shad} | \
    awk '($6 == "TRUE" || $6 == ""BOTH) && (!seen[$0]++){{a[$1"\\t"$5]+=1}}END{{for(i in a) print i,a[i]}}' | awk '{{print $1,$2,$3/$2}}' > {params.tmp1}

    # join with the whole table on cluster ID
    join -13 -21 -a1 <(sort --parallel={threads} -T {params.local_tmp} -k3,3 {output.sp_sh}.tmp) \
        <(sort --parallel={threads} -k1,1 {params.tmp1}) > {params.tmp2}

    # Add 0 for the clusters with no shadows and
    # reorder column output table: orf - length - cl_name - cl_size - prop_shadow - is.shadow - is.spurious
    awk -vOFS='\\t' '{{ for(i=1; i<=7; i++) if($i ~ /^ *$/) $i = 0 }};1' {params.tmp2} | \
        awk -vOFS='\\t' '{{print $2,$3,$1,$4,$7,$5,$6}}' > {output.sp_sh}.tmp

    if [[ {params.module} == "update" ]]; then
        # Download original dataset spurious and shadow info and combine with the new dataset
        if [[ ! -s {params.or_sp_sh} ]]; then
            wget https://ndownloader.figshare.com/files/23067131 -O {params.or_sp_sh}
        fi
        if [[ {params.eval_shared} == "true" ]]; then
            join -11 -21 -v1 <(zcat {params.or_sp_sh} |\
             sort -k1,1 --parallel {threads} -T {params.local_tmp}) \
            <(awk '{{print $1}}' {output.sp_sh}.tmp |\
             sort -k1,1 --parallel {threads} -T {params.local_tmp} ) > {params.tmp1}
             sed -i 's/ /\\t/g' {params.tmp1}
             cat {output.sp_sh}.tmp {params.tmp1} > {output.sp_sh}
        else
            cat {output.sp_sh}.tmp <(zcat {params.or_sp_sh}) > {output.sp_sh}
        fi
        rm {output.sp_sh}.tmp
    else
        mv {output.sp_sh}.tmp {output.sp_sh}
    fi

    rm {params.tmp1} {params.tmp2} {params.all_shad} {params.spur} {params.only_shad}
    """

run:
    shell("echo 'Spurious and Shadows detection DONE'")

shell:
    """
    Rscript --vanilla {params.report_maker} --basedir {params.basedir} \
                                            --outdir  {params.outdir} \
                                            --name {params.name_data} \
                                            --input {params.input_data} \
                                            --stage {params.stage} \
                                            --wf_report {params.wf_report} \
                                            --output {output.report} 2>{log.err} 1>{log.out}

    """