misc-scripts/align-clusters.nf at master · maxemil/misc-scripts · 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
params.cluster = ""
params.output_alignments = "alignments"
params.output_faa = "cluster_faa"

cluster = Channel.fromPath(params.cluster)

process sortClustersBySize {
  input:
  file faa from cluster

  output:
  file "small_cluster/${faa.simpleName}.faa" into small_faa optional true
  file "${faa.simpleName}.faa" into normal_faa optional true

  publishDir "${params.output_faa}", mode: 'copy', pattern: 'small_cluster/*'
  tag {"${faa.simpleName}"}

  script:
  """
  if [ "\$(grep -c '^>' $faa)" -lt 4 ]
  then
    mkdir small_cluster
    mv $faa small_cluster/${faa.simpleName}.faa
  else
    cp -L $faa ${faa.simpleName}.faa
  fi
  """
}

normal_faa.into{normal_faa_alignment; normal_faa_filering}


process alignCluster {
  input:
  file faa from normal_faa_alignment

  output:
  file "${faa.simpleName}.divvy.aln" into divvyied_alignments
  file "${faa.simpleName}.aln" into untreated_alignments

  // publishDir "${params.output_alignments}", mode: 'copy'
  tag {"${faa.simpleName}"}
  cpus 30

  script:
  """
  sed -i 's/*//g' $faa
  prequal $faa
  mafft-einsi --anysymbol --thread ${task.cpus} ${faa}.filtered > "${faa.simpleName}.aln"
  divvier -divvy -divvygap "${faa.simpleName}.aln"
  mv ${faa.simpleName}.aln.divvy.fas ${faa.simpleName}.divvy.aln
  """
}

filter_alignments = normal_faa_filering.combine(divvyied_alignments)

process filterAlignmentSize{
  input:
  set file(faa), file(divvy) from filter_alignments

  output:
  file "$divvy" into filtered_divvyied_alignments optional true
  file "$faa" into filtered_faa optional true
  file "small_cluster/*.faa" into filtered_small_cluster optional true

  publishDir "${params.output_alignments}", mode: 'copy', pattern: '*.aln'
  publishDir "${params.output_faa}", mode: 'copy', pattern: 'small_cluster/*'
  publishDir "${params.output_faa}", mode: 'copy', pattern: '*.faa'
  stageInMode 'copy'
  tag {"${faa.simpleName}"}

  when:
  "${faa.simpleName}" == "${divvy.simpleName}"

  script:
  """
  #! /usr/bin/env python3
  from Bio import SeqIO
  import os

  os.makedirs("small_cluster")
  missing_data = ['-', 'X']

  divvy = {rec.id:rec for rec in SeqIO.parse("$divvy", 'fasta')}
  faa = {rec.id:rec for rec in SeqIO.parse("$faa", 'fasta')}
  all_missing_data = []
  for recid, rec in divvy.items():
        if all([c in missing_data for c in rec.seq]):
            all_missing_data.append(recid)
  [divvy.pop(recid) for recid in all_missing_data]
  small = "" if len(divvy.keys()) > 3 else 'small_cluster/'
  with open("{}$faa".format(small), 'w') as faa_file, open("{}$divvy".format(small), 'w') as divvy_file:
      for recid, rec in faa.items():
          if recid in divvy.keys():
              SeqIO.write(divvy[recid], divvy_file, 'fasta')
              SeqIO.write(rec, faa_file, 'fasta')
          else:
              with open("small_cluster/{}.faa".format(recid.replace('.', '_')), 'w') as single:
                  SeqIO.write(rec, single, 'fasta')
  if small:
    os.remove("$divvy")
    os.remove("$faa")
  """

}