shell:
  '''
  mkdir -p {params.wd}
  mkdir -p {params.out}
  awk '{{print $1}}' {input.filtered_gff} | sort | uniq > {params.wd}_contigs
  split -l {params.split} {params.wd}_contigs {params.wd}
  for file in {params.wd}??; do awk -v file=$file '{{print $1 "\t" file}}' $file | awk -F'/|\t' '{{print $1 "\t" $NF}}'; done > {params.wd}_bedtools_list
  {input.bamdeal} modify bamAssign -i {params.bam} -a {params.wd}_bedtools_list -q 0 -o {params.wd}
  for file in {params.wd}??; do ./scripts/bedtools.sh $file {params.gff} {params.wd}; done > {output.coverage}
  {params.evenness_yes} for file in {params.wd}??; do ./scripts/evenness.sh $file {params.wd}; done
  {params.evenness_yes} cat {params.wd}??_evenness > {params.evenness}
  {params.evenness_yes} rm {params.wd}??_evenness
  rm {params.bam}
  rm -r {params.wd}
  rm {params.wd}_contigs
  rm {params.wd}_bedtools_list
  rm {params.wd}*coverage
  rm {params.wd}*gff
  rm {params.wd}??
  '''

shell:
  '''
  touch {output.touch}
  wget -P {params.location} {params.profiles}
  wget -P {params.location} {params.ko_list}
  tar -C {params.location} -xf {params.location}/profiles.tar.gz
  gunzip {params.location}/ko_list.gz
  rm {params.location}/profiles.tar.gz
  '''

shell:
  '''
  grep '^\*' {input.ko} | awk '{{print $2 "\t" $3}}' > {output.out}
  {params.splitting_multiples}{params.evenness_no}./scripts/splitting_multiples.sh {input.coverage} {output.out} {output.out} #include if splitting multiples and not using evenness
  {params.splitting_multiples}{params.evenness_yes}./scripts/splitting_multiples_evenness.sh {input.coverage} {output.out} {output.out} {params.evenness} {params.pid} #include if splitting multiples and using evenness
  {params.grouping_multiples}awk '{{print $1}}' {output.out} | sort | uniq -d > {output.out}_duplicates #include if grouping multiples
  {params.grouping_multiples}cat {output.out} | parallel --colsep '\t' ./scripts/multiples.sh {{1}} {output.out} {output.out}_duplicates | sort | uniq > {output.out}_multiples #include if grouping multiples
  {params.grouping_multiples}awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {output.out}_multiples {input.coverage} | awk 'BEGIN{{OFS=FS}} $3 == "" {{$3 = "NoHit"}} 1' > {output.out}_coverage #include if grouping multiples
  {params.grouping_multiples}rm {output.out}_duplicates #include if grouping multiples
  {params.grouping_multiples}rm {output.out}_multiples #include if grouping multiples
  {params.grouping_multiples}{params.evenness_no}awk '{{print $2 "\t" $3}}' {output.out}_coverage | awk '{{a[$2]+=$1}}END{{for(i in a) print i,a[i]}}' > {output.KOunt} #include if grouping multiples and not using evenness
  {params.grouping_multiples}{params.evenness_yes}awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {params.evenness} {output.out}_coverage | gawk -v pid="{params.pid}" '$4>pid' > {output.out}_coverage && awk '{{print $2 "\t" $3}}' {output.out}_coverage | awk '{{a[$2]+=$1}}END{{for(i in a) print i,a[i]}}' > {output.KOunt} #include if grouping multiples and using evenness
  mkdir -p {params.split_keggs}
  awk '{{print $3}}' {output.out}_coverage > {output.out}_coverage_tmp && while read l; do mkdir -p {params.split_keggs}/"$l"; done < {output.out}_coverage_tmp && rm {output.out}_coverage_tmp
  awk '{{print $1>f"{params.split_keggs}" $3 "/{wildcards.id}_" $3}}' {output.out}_coverage
  for file in {params.split_keggs}*/{wildcards.id}_*; do ./scripts/extracting_genes.sh $file {input.filtered_aa}; done
  grep NoHit {output.out}_coverage | awk '{{print $1}}' > {output.nohit}
  rm {output.out}_coverage
  '''

shell:
  '''
  grep '^\*' {input.ko} | awk '{{print $2 "\t" $3}}' > {output.out}
  {params.splitting_multiples}{params.evenness_no}./scripts/splitting_multiples.sh {input.coverage} {output.out} {output.out} #include if splitting multiples and not using evenness
  {params.splitting_multiples}{params.evenness_yes}./scripts/splitting_multiples_evenness.sh {input.coverage} {output.out} {output.out} {params.evenness} {params.pid} #include if splitting multiples and using evenness
  {params.grouping_multiples}awk '{{print $1}}' {output.out} | sort | uniq -d > {output.out}_duplicates #include if grouping multiples
  {params.grouping_multiples}cat {output.out} | parallel --colsep '\t' ./scripts/multiples.sh {{1}} {output.out} {output.out}_duplicates | sort | uniq > {output.out}_multiples #include if grouping multiples
  {params.grouping_multiples}awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {output.out}_multiples {input.coverage} | awk 'BEGIN{{OFS=FS}} $3 == "" {{$3 = "NoHit"}} 1' > {output.out}_coverage #include if grouping multiples
  {params.grouping_multiples}rm {output.out}_duplicates #include if grouping multiples
  {params.grouping_multiples}rm {output.out}_multiples #include if grouping multiples
  {params.grouping_multiples}{params.evenness_no}awk '{{print $2 "\t" $3}}' {output.out}_coverage | awk '{{a[$2]+=$1}}END{{for(i in a) print i,a[i]}}' > {output.KOunt} #include if grouping multiples and not using evenness
  {params.grouping_multiples}{params.evenness_yes}awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {params.evenness} {output.out}_coverage | gawk -v pid="{params.pid}" '$4>pid' > {output.out}_coverage && awk '{{print $2 "\t" $3}}' {output.out}_coverage | awk '{{a[$2]+=$1}}END{{for(i in a) print i,a[i]}}' > {output.KOunt} #include if grouping multiples and using evenness
  grep NoHit {output.out}_coverage | awk '{{print $1}}' > {output.nohit}
  rm {output.out}_coverage
  '''

shell:
  '''
  grep '^\*' {input.ko} | awk '{{print $2 "\t" $3}}' > {output.out}
  {params.splitting_multiples}{params.evenness_no}./scripts/splitting_multiples.sh {input.coverage} {output.out} {output.out} #include if splitting multiples and not using evenness
  {params.splitting_multiples}{params.evenness_yes}./scripts/splitting_multiples_evenness.sh {input.coverage} {output.out} {output.out} {params.evenness} {params.pid} #include if splitting multiples and using evenness
  {params.grouping_multiples}awk '{{print $1}}' {output.out} | sort | uniq -d > {output.out}_duplicates #include if grouping multiples
  {params.grouping_multiples}cat {output.out} | parallel --colsep '\t' ./scripts/multiples.sh {{1}} {output.out} {output.out}_duplicates | sort | uniq > {output.out}_multiples #include if grouping multiples
  {params.grouping_multiples}awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {output.out}_multiples {input.coverage} | awk 'BEGIN{{OFS=FS}} $3 == "" {{$3 = "NoHit"}} 1' > {output.out}_coverage #include if grouping multiples
  {params.grouping_multiples}rm {output.out}_duplicates #include if grouping multiples
  {params.grouping_multiples}rm {output.out}_multiples #include if grouping multiples
  {params.grouping_multiples}{params.evenness_no}awk '{{print $2 "\t" $3}}' {output.out}_coverage | awk '{{a[$2]+=$1}}END{{for(i in a) print i,a[i]}}' > {output.KOunt} #include if grouping multiples and not using evenness
  {params.grouping_multiples}{params.evenness_yes}awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {params.evenness} {output.out}_coverage | gawk -v pid="{params.pid}" '$4>pid' > {output.out}_coverage && awk '{{print $2 "\t" $3}}' {output.out}_coverage | awk '{{a[$2]+=$1}}END{{for(i in a) print i,a[i]}}' > {output.KOunt} #include if grouping multiples and using evenness
  rm {output.out}_coverage
  '''

shell:
  '''
  mkdir -p {params.tmp}
  if [ ! -e {params.kegglist} ]; then cat {input} | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  for i in {input}; do cp $i {params.tmp}; done
  ./scripts/abundance_matrix.sh {params.tmp} {params.list} {params.kegglist} > {output.results}
  rm -r {params.tmp}
  if [ -e {outdir}Touch/kegg_fasta ]; then rm {params.kegglist}; fi
  '''

shell:
  '''
  mkdir -p {params.tmp}
  if [ ! -e {params.kegglist} ]; then cat {input} | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  for i in {input}; do cp $i {params.tmp}; done
  ./scripts/abundance_matrix_noc.sh {params.tmp} {params.list} {params.kegglist} > {output.results}
  rm -r {params.tmp}
  if [ -e {outdir}Touch/kegg_fasta ]; then rm {params.kegglist}; fi
  '''

shell:
  '''
  mkdir -p {params.tmp}
  if [ ! -e {params.kegglist} ]; then cat {input} | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  for i in {input}; do cp $i {params.tmp}; done
  ./scripts/abundance_matrix_noref.sh {params.tmp} {params.list} {params.kegglist} > {output.results}
  rm -r {params.tmp}
  if [ -e {outdir}Touch/kegg_fasta ]; then rm {params.kegglist}; fi
  '''

shell:
  '''
  mkdir -p {params.cd}
  mkdir -p {params.cdhit}
  if [ ! -e {params.kegglist} ]; then cat {input} | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  cat {params.kegglist} | parallel "cat {params.split_dir}{{1}}/*faa > {params.cdhit}{{1}}.faa"
  rm -r {params.split_dir}
  if [ -e {params.mega} ]; then rm -r {params.mega}; fi
  touch {output}
  if [ -e {params.results} ]; then rm {params.kegglist}; fi
  '''

shell: 
  '''
  mkdir -p {params.lists}
  (cd {params.sample} && ls K*_1.0) > keggslist
  {params.grouping_multiples}(cd {params.sample} && ls Multiple*1.0) >> keggslist #include if grouping multiples
  split -d --number=l/{params.num} keggslist {params.lists}/keggslist
  rm keggslist
  touch {output}
  '''

shell:
  '''
  mkdir -p {params.mm}
  mkdir -p {params.mmseq}
  while read file; do ./scripts/mmseq_keggs.sh $file {params.cdhit} {params.mmseq}; done < {params.list}
  while read file; do ./scripts/count_mmseq_clusters.sh {params.mmseq}"$file"_0.9_all; done < {params.list}
  while read file; do ./scripts/count_mmseq_clusters.sh {params.mmseq}"$file"_0.5_all; done < {params.list}
  rm {params.list}
  touch {output}
  '''

shell:
  '''
  mkdir -p {params.mm}
  mkdir -p {params.mmseq}
  ./scripts/mmseq_keggs.sh NoHit.faa_1.0 {params.cdhit} {params.mmseq}
  ./scripts/count_mmseq_clusters.sh {params.mmseq}NoHit.faa_1.0_0.9_all
  ./scripts/count_mmseq_clusters.sh {params.mmseq}NoHit.faa_1.0_0.5_all
  touch {output}
  '''

shell:
  '''
  echo -e "KEGG,Number of clusters,Number of singletons,Number with multiple genes,Number of proteins\n$(cat {params.mmseq}*txt)" > {output}
  rm {params.mmseq}*txt
  rm -r {params.list}
  '''

shell:
  '''
  mkdir -p {params.wd}
  ./scripts/extracting_fastas.sh {input.nohit} {input.faa} > {params.faa}
  ./scripts/extracting_fastas.sh {input.nohit} {input.fa} > {params.fa}
  touch {output}
  '''

shell:
  '''
  mkdir -p {params.wd}
  ./scripts/extracting_fastas.sh {input.nohit} {input.faa} > {params.faa}
  ./scripts/extracting_fastas.sh {input.nohit} {input.fa} > {params.fa}
  touch {output}
  '''

shell:
  '''
  mkdir -p {params.wd}
  ./scripts/extracting_fastas.sh {input.nohit} {input.faa} > {params.faa}
  touch {output}
  '''

shell:
  '''
  diamond blastp --query {params.faa} --threads {params.threads} --outfmt {params.of} -d {params.db} > {output.raw}
  mmseqs easy-search {params.fa} {params.rna} {output.mmseq} {params.tmp} --threads {params.threads} --search-type 3 --format-output {params.ofm}
  cat {output.raw} {output.mmseq} > {params.cat}
  rm -r {params.tmp}
  awk -F '\t' '{{ print $0 "\t" ($4/$5) * 100 }}' {params.cat} | gawk -F '\t' '$NF>{params.min_qc} {{print $0}}' | gawk -F '\t' '$NF<{params.max_qc} {{print $0}}' | gawk -F '\t' '$3>{params.min_pid} {{print $0}}' > {output.fil}
  awk '! a[$1]++' {output.fil} > {output.fil}_tmp
  awk '{{print $1 "," $3 "," $10 "\t" $2}}' {output.fil}_tmp > {output.fil}_tmp.1
  awk '{{print $1 "," $3 "," $10 "\t" $2}}' {output.fil} > {output.fil}.1
  awk 'NR==FNR{{a[$1];next}}$1 in a' {output.fil}_tmp.1 {output.fil}.1 | tr ',' '\t' > {output.best}
  rm {output.fil}_tmp.1
  rm {output.fil}.1
  awk -F '_' '{{print $0 "\t" $NF}}' {output.best} | awk '{{print $1 "\t" $NF}}' | sort | uniq > {output.ko}
  awk '{{print $1}}' {output.ko} | uniq -c | awk '{{print $2 "\t" $1}}' > {output.ko}_count
  awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {input.cov} {output.ko} | awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {output.ko}_count - | awk '{{print $1 "\t" $2 "\t" ($3/$4)}}' > {output.ko}_coverage
  awk '{{a[$2]+=$3}}END{{for(i in a) print i,a[i]}}' {output.ko}_coverage | awk -F'-' '{{print $0 "\t" NF-1}}' | awk -v OFMT='%.5f' '{{print $1"\t",($2/$3)}}' | cut -d- -f2- | awk '{{gsub(/-/,"\t"$2",")}}1' | tr ',' '\n' | awk -v OFMT='%.5f' '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {output.KOunt}
  awk '{{print $1}}' {output.ko}_coverage > {output.ko}_coverage_tmp
  cat {input.nohit} {output.ko}_coverage_tmp | sort | uniq -u > {output.no}
  rm {output.ko}_coverage_tmp
  rm {output.ko}_coverage
  rm {output.ko}_count
  '''

shell:
  '''
  diamond blastp --query {params.faa} --threads {params.threads} --outfmt {params.of} -d {params.db} > {output.raw}
  awk -F '\t' '{{ print $0 "\t" ($4/$5) * 100 }}' {output.raw} | gawk -F '\t' '$NF>{params.min_qc} {{print $0}}' | gawk -F '\t' '$NF<{params.max_qc} {{print $0}}' | gawk -F '\t' '$3>{params.min_pid} {{print $0}}' > {output.fil}
  awk '! a[$1]++' {output.fil} > {output.fil}_tmp
  awk '{{print $1 "," $3 "," $10 "\t" $2}}' {output.fil}_tmp > {output.fil}_tmp.1
  awk '{{print $1 "," $3 "," $10 "\t" $2}}' {output.fil} > {output.fil}.1
  awk 'NR==FNR{{a[$1];next}}$1 in a' {output.fil}_tmp.1 {output.fil}.1 | tr ',' '\t' > {output.best}
  rm {output.fil}_tmp.1
  rm {output.fil}.1
  awk -F '_' '{{print $0 "\t" $NF}}' {output.best} | awk '{{print $1 "\t" $NF}}' | sort | uniq > {output.ko}
  awk '{{print $1}}' {output.ko} | uniq -c | awk '{{print $2 "\t" $1}}' > {output.ko}_count
  awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {input.cov} {output.ko} | awk 'NR==FNR{{a[$1]=$2;next}}{{print $0,a[$1]}}' {output.ko}_count - | awk '{{print $1 "\t" $2 "\t" ($3/$4)}}' > {output.ko}_coverage
  awk '{{a[$2]+=$3}}END{{for(i in a) print i,a[i]}}' {output.ko}_coverage | awk -F'-' '{{print $0 "\t" NF-1}}' | awk -v OFMT='%.5f' '{{print $1"\t",($2/$3)}}' | cut -d- -f2- | awk '{{gsub(/-/,"\t"$2",")}}1' | tr ',' '\n' | awk -v OFMT='%.5f' '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {output.KOunt}
  awk '{{print $1}}' {output.ko}_coverage > {output.ko}_coverage_tmp
  cat {input.nohit} {output.ko}_coverage_tmp | sort | uniq -u > {output.no}
  rm {output.ko}_coverage_tmp
  rm {output.ko}_coverage
  rm {output.ko}_count
  '''

shell:
  '''
  mkdir -p {params.nohit}
  awk -F '=|\t|;' '{{print $1 "_" $10 "\t" ($4-1) "\t" $5}}' {input.gff} > {params.bed}
  awk 'NR==FNR{{a[$1]=$2;b[$1]=$3;next}}{{print $0 "\t" a[$1] "\t" b[$1]}}' {params.bed} {input.nohit} > {params.nobed}_tmp
  rm {params.bed}
  awk '{{print $1}}' {params.nobed}_tmp | rev | cut -d "_" -f 2- | rev > {params.nobed}_tmp.1
  awk '{{print $2 "\t" $3}}' {params.nobed}_tmp > {params.nobed}_tmp.2
  paste {params.nobed}_tmp.1 {params.nobed}_tmp.2 > {params.nobed}
  rm {params.nobed}_tmp
  rm {params.nobed}_tmp.1
  rm {params.nobed}_tmp.2
  bedtools intersect -a {params.nobed} -b {params.bam} -bed -wb | awk '{{print $7}}' | sort | uniq > {params.reads}
  rm {params.nobed}
  rm {params.bam}
  gawk -F '/' '$NF==1' {params.reads} > {params.reads}_R1
  gawk -F '/' '$NF==2' {params.reads} > {params.reads}_R2
  zcat {params.r1} | seqtk subseq - {params.reads}_R1 > {params.fastq}_R1
  zcat {params.r2} | seqtk subseq - {params.reads}_R2 > {params.fastq}_R2
  gzip {params.fastq}_R1
  gzip {params.fastq}_R2
  rm {params.reads}_R1
  rm {params.reads}_R2
  cat {params.fastq}_R1.gz {params.fastq}_R2.gz > {params.fastq}
  rm {params.fastq}_R1.gz
  rm {params.fastq}_R2.gz
  touch {output}
  '''

shell:
  '''
  mkdir -p {params.nohit}
  awk -F '=|\t|;' '{{print $1 "_" $10 "\t" ($4-1) "\t" $5}}' {input.gff} > {params.bed}
  awk 'NR==FNR{{a[$1]=$2;b[$1]=$3;next}}{{print $0 "\t" a[$1] "\t" b[$1]}}' {params.bed} {input.nohit} > {params.nobed}_tmp
  rm {params.bed}
  awk '{{print $1}}' {params.nobed}_tmp | rev | cut -d "_" -f 2- | rev > {params.nobed}_tmp.1
  awk '{{print $2 "\t" $3}}' {params.nobed}_tmp > {params.nobed}_tmp.2
  paste {params.nobed}_tmp.1 {params.nobed}_tmp.2 > {params.nobed}
  rm {params.nobed}_tmp
  rm {params.nobed}_tmp.1
  rm {params.nobed}_tmp.2
  bedtools intersect -a {params.nobed} -b {params.bam} -bed -wb | awk '{{print $7}}' | sort | uniq > {params.reads}
  rm {params.nobed}
  rm {params.bam}
  gawk -F '/' '$NF==1' {params.reads} > {params.reads}_R1
  gawk -F '/' '$NF==2' {params.reads} > {params.reads}_R2
  zcat {params.r1} | seqtk subseq - {params.reads}_R1 > {params.fastq}_R1
  zcat {params.r2} | seqtk subseq - {params.reads}_R2 > {params.fastq}_R2
  gzip {params.fastq}_R1
  gzip {params.fastq}_R2
  rm {params.reads}_R1
  rm {params.reads}_R2
  cat {params.fastq}_R1.gz {params.fastq}_R2.gz > {params.fastq}
  rm {params.fastq}_R1.gz
  rm {params.fastq}_R2.gz
  touch {output}
  '''

shell:
  '''
  mkdir -p {params.wd}
  diamond blastx --query {params.fq} --threads {params.threads} --outfmt {params.of} -d {input.db} > {params.nohit}
  mmseqs easy-search {params.fq} {input.rna} {params.mmseq} {params.tmp} --threads {params.threads} --search-type 3 --format-output {params.ofm}
  cat {params.nohit} {params.mmseq} > {params.cat}
  rm {params.nohit}
  rm {params.mmseq}
  rm -r {params.tmp}
  ./scripts/unmapped_reads_ko_100.sh {params.cat} {params.wd} {output.besthits} {input.bg} {output.ko}
  awk '{{print $1}}' {output.besthits} | sort | uniq | cat - {params.reads} | sort | uniq -u > {params.missing}
  rm {params.cat}
  rm {params.reads}
  '''

shell:
  '''
  mkdir -p {params.wd}
  diamond blastx --query {params.fq} --threads {params.threads} --outfmt {params.of} -d {input.db} > {params.nohit}
  mmseqs easy-search {params.fq} {input.rna} {params.mmseq} {params.tmp} --threads {params.threads} --search-type 3 --format-output {params.ofm}
  cat {params.nohit} {params.mmseq} > {params.cat}
  rm {params.nohit}
  rm {params.mmseq}
  rm -r {params.tmp}
  ./scripts/unmapped_reads_ko_100.sh {params.cat} {params.wd} {output.besthits} {input.bg} {output.ko}
  awk '{{print $1}}' {output.besthits} | sort | uniq | cat - {params.reads} | sort | uniq -u > {params.missing}
  rm {params.cat}
  rm {params.reads}
  '''

shell:
  '''
  mkdir -p {params.wd}
  diamond blastx --query {params.fq} --threads {params.threads} --outfmt {params.of} -d {input.db} > {params.nohit}
  ./scripts/unmapped_reads_ko_100.sh {params.nohit} {params.wd} {output.besthits} {input.bg} {output.ko}
  rm {params.nohit}
  '''

shell:
  '''
  gawk -F '/' '$NF==1' {params.missing} > {params.missing}_R1
  gawk -F '/' '$NF==2' {params.missing} > {params.missing}_R2
  zcat {params.fq} | seqtk subseq - {params.missing}_R1 | gzip > {params.r1}
  zcat {params.fq} | seqtk subseq - {params.missing}_R2 | gzip > {params.r2}
  rm {params.fq}
  rm {params.missing}
  rm {params.missing}_R1
  rm {params.missing}_R2
  kallisto quant -i {input.ref} -o {params.kal} -t {params.threads} {params.r1} {params.r2}
  awk -v len=$(cat {params.r1} {params.r2} | awk 'NR%4==2{{sum+=length($0)}}END{{print sum/(NR/4)}}') 'NR>1{{print $1"\t"(len*$4)/$2}}' {params.kal}/abundance.tsv | gawk '$2!=0' > {params.kal}/tmp
  awk '{{print $1}}' {params.kal}/tmp | awk -F'-' '{{print NF-1}}' | paste {params.kal}/tmp - | awk -v OFMT='%.10f' -v OFS='\t' '{{print $1,($2/$3)}}' | awk '{{gsub(/-/,"\t"$2",")}}1' | cut -d, -f2- | sed 's/,$//g' | tr ',' '\n' | awk -v OFMT='%.10f' '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {output.ko}
  rm {params.kal}/tmp
  rm {params.r1}
  rm {params.r2}
  rm -r {params.kal}
  gzip {input.besthits}
  '''

shell:
  '''
  mkdir -p {params.wd}
  mkdir -p {params.mm}
  mkdir -p {params.mm2}
  mkdir -p {params.tmp}
  zcat {params.r1} {params.r2} | awk 'NR%4==1 {{print substr($1,2)}}' | cat - {params.mapped} | sort --parallel=8 | uniq -u > {params.unmapped}
  rm {params.mapped}
  zcat {params.r1} {params.r2} | seqtk subseq - {params.unmapped} | gzip > {params.fastq}
  diamond blastx --query {params.fastq} --threads {params.threads} --outfmt 6 qseqid sseqid pident qlen slen qstart qend sstart send evalue -d {input.db} > {params.dia}
  mmseqs easy-search {params.fastq} {input.rna} {params.mmseq} {params.tmp} --threads {params.threads} --search-type 3 --format-output {params.ofm}
  rm -r {params.tmp}
  cat {params.dia} {params.mmseq} > {params.cat}
  rm {params.dia}
  rm {params.mmseq}
  ./scripts/unmapped_reads_ko_100.sh {params.cat} {params.wd} {output.besthits} {input.bg} {output.ko}
  rm {params.cat}
  '''

shell:
  '''
  mkdir -p {params.wd}
  zcat {params.r1} {params.r2} | awk 'NR%4==1 {{print substr($1,2)}}' | cat - {params.mapped} | sort | uniq -u > {params.unmapped}
  rm {params.mapped}
  zcat {params.r1} {params.r2} | seqtk subseq - {params.unmapped} | gzip > {params.fastq}
  rm {params.unmapped}
  diamond blastx --query {params.fastq} --threads {params.threads} --outfmt 6 qseqid sseqid pident qlen slen qstart qend sstart send evalue -d {input.db} > {params.dia}
  ./scripts/unmapped_reads_ko_100.sh {params.dia} {params.wd} {output.besthits} {input.bg} {output.ko}
  rm {params.dia}
  '''

shell:
  '''
  mkdir -p {params.nohit}
  awk '{{print $1}}' {input.besthits} | sort | uniq | cat - {params.unmapped} | sort | uniq -u > {params.still_missing}
  rm {params.unmapped}
  gawk -F '/' '$NF==1' {params.still_missing} > {params.still_missing}_R1
  gawk -F '/' '$NF==2' {params.still_missing} > {params.still_missing}_R2
  zcat {params.fq} | seqtk subseq - {params.still_missing}_R1 | gzip > {params.r1}
  zcat {params.fq} | seqtk subseq - {params.still_missing}_R2 | gzip > {params.r2}
  rm {params.fq}
  rm {params.still_missing}_R1
  rm {params.still_missing}_R2
  kallisto quant -i {input.ref} -o {params.kal} -t {params.threads} {params.r1} {params.r2}
  awk -v len=$(zcat {params.r1} {params.r2} | awk 'NR%4==2{{sum+=length($0)}}END{{print sum/(NR/4)}}') 'NR>1{{print $1"\t"(len*$4)/$2}}' {params.kal}/abundance.tsv | gawk '$2!=0' > {params.kal}/tmp
  awk '{{print $1}}' {params.kal}/tmp | awk -F'-' '{{print NF-1}}' | paste {params.kal}/tmp - | awk -v OFMT='%.10f' -v OFS='\t' '{{print $1,($2/$3)}}' | awk '{{gsub(/-/,"\t"$2",")}}1' | cut -d, -f2- | sed 's/,$//g' | tr ',' '\n' | awk -v OFMT='%.10f' '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {output.ko}
  rm {params.kal}/tmp
  rm {params.r1}
  rm {params.r2}
  gzip {input.besthits}
  '''

shell:
  '''
  mkdir {params.tmp}
  for i in {input.kofam}; do cp $i {params.tmp}; done
  for i in {input.nohit_p}; do cp $i {params.tmp}; done
  for i in {input.nohit_r}; do cp $i {params.tmp}; done
  for i in {input.kal}; do cp $i {params.tmp}; done
  for i in {input.unmap}; do cp $i {params.tmp}; done
  for i in {input.kal_unmap}; do cp $i {params.tmp}; done
  sed -i '/NoHit/d' {params.tmp}/*
  if [ ! -e {params.kegglist} ]; then cat {params.tmp}* | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  for i in {params.id}; do echo $i; done > {params.tmp}samples
  while read file; do cat {params.tmp}"$file"_* | awk '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {params.tmp}"$file"_catted; done < {params.tmp}samples
  ./scripts/abundance_matrix_catted.sh {params.tmp} {params.list} {params.kegglist} > {output.ann_results}
  rm {params.kegglist}
  rm -r {params.tmp}
  '''

shell:
  '''
  mkdir {params.tmp}
  for i in {input.kofam}; do cp $i {params.tmp}; done
  for i in {input.nohit_p}; do cp $i {params.tmp}; done
  for i in {input.nohit_r}; do cp $i {params.tmp}; done
  for i in {input.kal}; do cp $i {params.tmp}; done
  for i in {input.unmap}; do cp $i {params.tmp}; done
  for i in {input.kal_unmap}; do cp $i {params.tmp}; done
  sed -i '/NoHit/d' {params.tmp}/*
  if [ ! -e {params.kegglist} ]; then cat {params.tmp}* | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  for i in {params.id}; do echo $i; done > {params.tmp}samples
  while read file; do cat {params.tmp}"$file"_* | awk '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {params.tmp}"$file"_catted; done < {params.tmp}samples
  ./scripts/abundance_matrix_catted.sh {params.tmp} {params.list} {params.kegglist} > {output.ann_results}
  rm {params.kegglist}
  rm -r {params.tmp}
  '''

shell:
  '''
  mkdir {params.tmp}
  for i in {input.kofam}; do cp $i {params.tmp}; done
  for i in {input.nohit_p}; do cp $i {params.tmp}; done
  for i in {input.nohit_reads}; do cp $i {params.tmp}; done
  for i in {input.unmap}; do cp $i {params.tmp}; done
  sed -i '/NoHit/d' {params.tmp}/*
  if [ ! -e {params.kegglist} ]; then cat {params.tmp}* | awk '{{print $1}}' | sort | uniq > {params.kegglist}; fi
  for i in {params.id}; do echo $i; done > {params.tmp}samples
  while read file; do cat {params.tmp}"$file"_* | awk '{{a[$1]+=$2}}END{{for(i in a) print i,a[i]}}' > {params.tmp}"$file"_catted; done < {params.tmp}samples
  ./scripts/abundance_matrix_catted.sh {params.tmp} {params.list} {params.kegglist} > {output.ann_results}
  rm -r {params.tmp}
  rm {params.kegglist}
  '''