forked from anttonalberdi/holoflow
-
Notifications
You must be signed in to change notification settings - Fork 1
/
holoflow.py
233 lines (149 loc) · 7.81 KB
/
holoflow.py
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
import argparse
import subprocess
import os
import sys
import ruamel.yaml
###########################
#Argument parsing
###########################
parser = argparse.ArgumentParser(description='Runs holoflow pipeline.')
parser.add_argument('-f', help="input.txt file", dest="input_txt", required=True)
parser.add_argument('-d', help="temp files directory path", dest="work_dir", required=True)
parser.add_argument('-w', help="chosen workflow", dest="workflow", required=True)
parser.add_argument('-c', help="config file", dest="config_file", required=True)
parser.add_argument('-t', help="threads", dest="threads", required=True)
args = parser.parse_args()
in_f=args.input_txt
path=args.work_dir
workflow=args.workflow
config=args.config_file
cores=args.threads
# retrieve current directory
file = os.path.dirname(sys.argv[0])
curr_dir = os.path.abspath(file)
#Append current directory to .yaml config for standalone calling
yaml = ruamel.yaml.YAML()
yaml.explicit_start = True
with open(str(config), 'r') as config_file:
data = yaml.load(config_file)
with open(str(config), 'w') as config_file:
data['holopath'] = str(curr_dir)
dump = yaml.dump(data, config_file)
###########################
## Functions
###########################
###########################
###### PREPARE GENOMES FUNCTIONS
###########################
###### PREPROCESSING FUNCTIONS
def in_out_preprocessing(path,in_f):
"""Generate output names files from input.txt. Rename and move
input files where snakemake expects to find them if necessary."""
# Define input directory and create it if not exists "00-InputData"
in_dir = os.path.join(path,"PPR_00-InputData")
if not os.path.exists(in_dir):
os.makedirs(in_dir)
with open(in_f,'r') as in_file:
# Generate desired output file names from input.txt
read = 0
output_files=''
final_temp_dir="PPR_03-MappedToReference"
lines = in_file.readlines() # Read input.txt lines
for file in lines:
if not (file.startswith('#')):
file = file.strip('\n').split(' ') # Create a list of each line
read+=1 # every sample will have two reads, keep the name of the file but change the read
# Add an output file based on input.txt info to a list for Snakemake command
output_files+=(path+"/"+final_temp_dir+"/"+file[0]+"_"+str(read)+".fastq ")
# Move files to new dir "00-InputData" and change file names for 1st column in input.txt
# if the current input file names do not match the designed ones in input.txt
filename=file[2] # current input file path and name
desired_filename='"'+in_dir+'/'+file[0]+'_'+str(read)+'.fastq"' # desired input file path and name specified in input.txt
if not ((filename == desired_filename) and (os.path.exists(str(desired_filename)))):
if filename.endswith('.gz'): # uncompress input file if necessary
uncompressCmd='gunzip -c '+filename+' > '+desired_filename+''
subprocess.check_call(uncompressCmd, shell=True)
else: # else just move the input file to "00-InputData" with the new name
copyfilesCmd='cp '+filename+' '+desired_filename+''
subprocess.check_call(copyfilesCmd, shell=True)
if read == 2:
read=0 # two read files for one sample finished, new sample
# Add stats output file only once per sample
output_files+=(path+"/"+final_temp_dir+"/"+file[0]+".stats ")
return output_files
def run_preprocessing(in_f, path, config, cores):
"""Run snakemake on shell"""
# Define output names
out_files = in_out_preprocessing(path,in_f)
curr_dir = os.path.dirname(sys.argv[0])
holopath = os.path.abspath(curr_dir)
path_snkf = os.path.join(holopath,'workflows/preprocessing/Snakefile')
# Run snakemake
prep_snk_Cmd = 'snakemake -s '+path_snkf+' '+out_files+' --configfile '+config+' --cores '+cores+''
subprocess.check_call(prep_snk_Cmd, shell=True)
print("Have a nice run!\n\t\tHOLOFOW Preprocessing starting")
###########################
###### METAGENOMICS FUNCTIONS
def in_out_metagenomics(path,in_f):
"""Generate output names files from input.txt. Rename and move
input files where snakemake expects to find them if necessary."""
in_dir = os.path.join(path,"PPR_03-MappedToReference")
if not os.path.exists(in_dir):
os.makedirs(in_dir)
with open(in_f,'r') as in_file:
# Paste desired output file names from input.txt
read = 0
output_files=''
final_temp_dir="MIA_03-Binning"
lines = in_file.readlines() # Read input.txt lines
for file in lines:
if not (file.startswith('#')):
file = file.strip('\n').split(' ') # Create a list of each line
read+=1 # every sample will have two reads, keep the name of the file but change the read
# Add an output file based on input.txt info to a list for Snakemake command
output_files+=(path+"/"+final_temp_dir+"/"+file[0]+"_dastool/"+file[0])
# Move files to new dir "PPR_03-MappedToReference/" and change file names for 1st column in input.txt
# if the current input file names do not match the designed ones in input.txt
filename=file[2] # current input file path and name
desired_filename='"'+in_dir+'/'+file[0]+'_'+str(read)+'.fastq"' # desired input file path and name specified in input.txt
if not ((filename == desired_filename) and (os.path.exists(str(desired_filename)))):
if filename.endswith('.gz'): # uncompress input file if necessary
uncompressCmd='gunzip -c '+filename+' > '+desired_filename+''
subprocess.check_call(uncompressCmd, shell=True)
else: # else just move the input file to "00-InputData" with the new name
copyfilesCmd='cp '+filename+' '+desired_filename+''
subprocess.check_call(copyfilesCmd, shell=True)
if read == 2: # two read files for one sample finished, new sample
read=0
# Add stats output file only once per sample
output_files+=(path+"/"+final_temp_dir+"/"+file[0]+".stats ")
return output_files
def run_metagenomics(in_f, path, config, cores):
"""Run snakemake on shell"""
# Define output names
out_files = in_out_metagenomics(path,in_f)
curr_dir = os.path.dirname(sys.argv[0])
holopath = os.path.abspath(curr_dir)
path_snkf = os.path.join(holopath,'workflows/metagenomics/individual_assembly/Snakefile')
# Run snakemake
mtg_snk_Cmd = 'snakemake -s '+path_snkf+' '+out_files+' --configfile '+config+' --cores '+cores+''
subprocess.check_call(mtg_snk_Cmd, shell=True)
print("Have a nice run!\n\t\tHOLOFOW Metagenomics starting")
###########################
###### GENOMICS FUNCTIONS
###########################
#### Snakemake pipeline run - load required modules
###########################
load_modulesCmd='module unload gcc/5.1.0 && module load tools anaconda3/4.4.0'
subprocess.check_call(load_modulesCmd, shell=True)
###########################
#### Workflows running
###########################
# 0 # Prepare genomes workflow
# 1 # Preprocessing workflow
if workflow == "preprocessing":
run_preprocessing(in_f, path, config, cores)
# 2 # Metagenomics workflow
if workflow == "metagenomics": # DATA HAS TO BE PREPROCESSED!
run_metagenomics(in_f, path, config, cores)
# 3 # Genomics workflow