forked from EisenRa/holoflow
-
Notifications
You must be signed in to change notification settings - Fork 0
/
tmp_metagenomics_IA.py
133 lines (97 loc) · 4.77 KB
/
tmp_metagenomics_IA.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
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('-c', help="config file", dest="config_file", required=False)
parser.add_argument('-l', help="pipeline log file", dest="log", required=False)
parser.add_argument('-t', help="threads", dest="threads", required=True)
args = parser.parse_args()
in_f=args.input_txt
path=args.work_dir
cores=args.threads
# retrieve current directory
file = os.path.dirname(sys.argv[0])
curr_dir = os.path.abspath(file)
if not (args.config_file):
config = os.path.join(os.path.abspath(curr_dir),"workflows/preparegenomes/config.yaml")
else:
config=args.config_file
if not (args.log):
log = os.path.join(path,"Holoflow_metagenomics.log")
else:
log=args.log
#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)
data['logpath'] = str(log)
dump = yaml.dump(data, config_file)
###########################
## Functions
###########################
###########################
###### 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_06-BinScaffolding"
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
# 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=str(file[2]) # current input file path and name
desired_filename=os.path.join(str(in_dir),''+str(file[0])+'_'+str(read)+'.fastq') # desired input file path and name specified in input.txt
if not (os.path.exists(str(desired_filename))):
print(filename == desired_filename)
print(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 an output file based on input.txt info to a list for Snakemake command
output_files+=(path+"/"+final_temp_dir+"/"+file[0]+"/Scaffolded_bins ")
# Add stats output file only once per sample
#output_files+=(path+"/MIA_01-Assembly/"+file[0]+".stats ")
# change for
#####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 = 'module unload gcc/5.1.0 && module load tools anaconda3/4.4.0 && snakemake -s '+path_snkf+' -k '+out_files+' --configfile '+config+' --cores '+cores+''
subprocess.check_call(mtg_snk_Cmd, shell=True)
print("Have a nice run!\n\t\tHOLOFOW Metagenomics starting")
###########################
#### Workflows running
###########################
# 2 # Metagenomics workflow
run_metagenomics(in_f, path, config, cores)