New upstream version 2.0.3+dfsg

#!/bin/bash

#conda-env update -f environment/packages.yaml

source activate qcumber

python3 QCumber.py $*

source deactivate qcumber

#!/usr/bin/env python3

__author__ = 'LieuV'

__version__ = "2.0.3"

import configparser

import argparse

import re

import getpass

import warnings

import os

import json

import sys

from itertools import groupby

from os.path import  splitext, join, isfile, dirname, abspath, exists, isdir, basename

from collections import OrderedDict

import subprocess

from pandas import read_csv

import snakemake

import datetime

import yaml

def get_basename(abs_name):

	return os.path.basename(os.path.splitext(os.path.splitext(os.path.basename(abs_name))[0])[0])

def getFilenameWithoutExtension(string, getBase=False):

	if getBase:

		string = os.path.basename(string)

	string = os.path.splitext(string)[0]

	i = 0

	while os.path.splitext(string)[-1] in [".gz", ".gzip", ".zip", ".bz", ".fasta", ".fastq", ".bam"]:

		string = os.path.splitext(string)[0]

	return string

def get_setname(filename, base=True, grouping=True):

	if not grouping:

		return get_basename(filename)

	try:

		if base:

			filename = get_basename(filename)

		sep = iter([";", "\t"])

		if arguments["rename"]:

			rename_dict = read_csv(arguments["rename"], index_col=0, header=None)

			while rename_dict.columns.__len__() == 0 or sep.__length_hint__() != 0:

				rename_dict = read_csv(arguments["rename"], index_col=0, header=None, sep=next(sep))

			if rename_dict.columns.__len__() == 0:

				warnings.warn(

					"Problems reading rename-file %s. Valid delimiters are ';', ',' or '\\t'." % arguments["rename"],

					UserWarning, stacklevel=2)

			new_name = [(x, rename_dict.ix[x][1]) for x in rename_dict.index if filename.startswith(x)]

			if len(new_name) != 1:

				warnings.warn("Could not find unique renames. Found %s for %s" % (new_name, filename), UserWarning,

				              stacklevel=2)

				filename = getFilenameWithoutExtension(filename, True)

			else:

				filename = filename.replace(new_name[0][0], new_name[0][1])

	except:

		print("Couldnt rename sample files.")

	try:

		paired_reads_pattern = base_pattern

		setname_pattern = re.search(paired_reads_pattern, os.path.basename(filename))

		if setname_pattern:

			return setname_pattern.group("setname")

		else:

			return filename

	except:

		print("Problems getting samplenames: %s" % filename)

		return filename

def get_defaults():

	if arguments["trimBetter"] == "assembly":

		# if arguments["forAssembly"]:

		if not arguments["trimBetter_threshold"]:

			arguments["trimBetter_threshold"] = parameter["forAssembly." + arguments['technology']][

				"trimBetter_threshold"]

		if not arguments["trimOption"]:

			arguments["trimOption"] = parameter["forAssembly." + arguments['technology']]["trimOption"]

	elif arguments["trimBetter"] == "mapping":

		# elif arguments["forMapping"]:

		if not arguments["trimBetter_threshold"]:

			arguments["trimBetter_threshold"] = parameter["forMapping." + arguments['technology']][

				"trimBetter_threshold"]

		if not arguments["trimOption"]:

			arguments["trimOption"] = parameter["forMapping." + arguments['technology']]["trimOption"]

	if arguments["trimBetter_threshold"] =="default":

		arguments["trimBetter_threshold"] = parameter["Trimmomatic"]["trimBetter_threshold"]

	if not arguments["trimOption"]:

		arguments["trimOption"] = parameter["Trimmomatic"]["trimOption"]

	if arguments["trimBetter_threshold"]:

		arguments["trimBetter_threshold"] = parameter["Trimmomatic"]["trimBetter_threshold"]

def get_version(cmd, jar=None):

	try:

		return re.match(r"(?P<version>.*)\n",subprocess.check_output(cmd, shell=True, stderr=subprocess.PIPE).decode("utf-8")).group("version")

	except:

		try:

			try:

				return subprocess.check_output(cmd, shell=True, stderr=subprocess.PIPE).decode("utf-8")

			except:

				from debpackageinfo import get_upstream_version

				return get_upstream_version(jar)

		except:

			return "NaN"

def get_input():

	bam_ext = [x.strip(" ") for x in parameter["Fileextension"]["bam"].split(",")]

	fastq_ext = [x.strip(" ") for x in parameter["Fileextension"]["fastq"].split(",")]

	sample_dict = OrderedDict()

	all_files = []

	name_dict = {}

	join_reads = {}

	type = "PE"

	join_lanes = OrderedDict()

	if arguments["r1"]:

		assert os.path.getsize(arguments["r1"]) != 0, "File %s is empty." % arguments["r1"]

		if any([arguments["r1"].endswith(ext) for ext in bam_ext]):

			arguments["technology"] = "IonTorrent"

		else:

			arguments["technology"] = "Illumina"

		if arguments["r2"]:

			assert os.path.getsize(arguments["r2"]) != 0, "File %s is empty." % arguments["r2"]

			sample_dict[get_setname(arguments["r1"])] = [abspath(arguments["r1"]), abspath(arguments["r2"])]

			name_dict[get_basename(arguments["r1"]) ] = get_setname(arguments["r1"]) + "_R1"

			name_dict[get_basename(arguments["r2"]) ] = get_setname(arguments["r2"]) + "_R2"

		else:

			type = "SE"

			sample_dict[get_setname(arguments["r1"])] = [abspath(arguments["r1"])]

			name_dict[get_basename(arguments["r1"])] = get_setname(arguments["r1"])

	else:

		if os.path.isdir(arguments["input"][0]):

			for root, dirs, files in os.walk(arguments["input"][0]):

				for file in files:

					if any([file.endswith(ext) for ext in fastq_ext + bam_ext]):

						if os.path.getsize(join(root, file)) != 0:

							all_files.append(join(root, file))

						else:

							warnings.warn("Skip empty file %s" % file, stacklevel=2)

		else:

			all_files = arguments["input"]

			assert all_files, "Check input again. No files found for pattern %s " % arguments["input"]

		if len([x for x in all_files if any([ext in x for ext in bam_ext])]) == len(all_files):

			arguments["technology"] = "IonTorrent"

		else:

			arguments["technology"] = "Illumina"

		if (len(all_files) == 0):

			sys.exit(str(arguments["input"]) + " does not contain fastq or bam files.")

		# find read pairs

		all_files = sorted(list(all_files))

		if all([re.search(base_pattern, x) for x in all_files]):

			for setname, files in groupby(all_files, key=lambda x: re.search(base_pattern, x).group("setname")):

				read_pairs = dict()

				setname = get_setname(setname)

				for lane, lane_file in groupby(list(files), key=lambda x: re.search(base_pattern, x).group("lane")):

					read_pairs[lane] = []

					for readgroup, readfiles in groupby(list(lane_file),

					                                    key=lambda x: re.search(base_pattern, x).group("read")):

						readfiles = list(readfiles)

						if len(readfiles) != 0:

							if len(readfiles) > 1:

								concat_reads = "QCResults/tmp/join_reads/" + "_".join(

									re.search(base_pattern, basename(readfiles[0])).groups()[:-2]) + "_000.fastq.gz"

								join_reads[concat_reads] = [abspath(x) for x in readfiles]

								readfiles = concat_reads

							else:

								readfiles = abspath(readfiles[0])

							read_pairs[lane].append(readfiles)

				# Multiple lanes

				if len(read_pairs) > 1:

					join_lanes[setname] = []

				for key in sorted(read_pairs.keys()):

					samplename = setname + "_" + key

					sample_dict[samplename] = read_pairs[key]

					if len(read_pairs[key])==2:#type == "PE":

						name_dict[get_basename(read_pairs[key][0])] = get_setname(samplename) + "_R1"

						name_dict[get_basename(read_pairs[key][1])] = get_setname(samplename) + "_R2"

					else:

						type = "SE"

						name_dict[get_basename(read_pairs[key][0])] = get_setname(samplename, grouping = False)

					if len(read_pairs) > 1:

						join_lanes[setname].append(samplename)

		else:  # treat each file as sample

			print("Treat files as single end")

			sample_dict = OrderedDict([[get_setname(x, grouping=False), [abspath(x)]] for x in all_files])

			name_dict = dict([get_basename(x), get_setname(x, grouping=False)] for x in all_files)

			type = "SE"

	return type, sample_dict, join_lanes, name_dict, join_reads

def check_input_validity():

	if arguments["reference"]:

		if not os.path.exists(arguments["reference"]):

			sys.exit("Reference does not exist.")

		try:

			seq_record = open(arguments["reference"], "r").readline()

			assert seq_record.startswith(">"), "Reference file is not valid."

		except:

			sys.exit("Reference file is not valid.")

	#

	# Check validity of Kraken DB

	if not arguments["nokraken"]:

		if not os.path.exists(arguments["kraken_db"]):

			sys.exit(arguments["kraken_db"] + " does not exist. Enter a valid database for kraken")

		else:

			if "database.kdb" not in os.listdir(arguments["kraken_db"]):

				sys.exit("database " + arguments["kraken_db"] + " does not contain necessary file database.kdb")

	#

	# Check input

	if arguments["input"]:

		if not all([os.path.exists(x) for x in arguments["input"]]):

			sys.exit(str(arguments["input"]) + " does not exist.")

	else:

		if not arguments["r1"]:

			sys.exit("Pleaser enter an input file (--input or -1/-2)")

		if not os.path.isfile(arguments["r1"]):

			sys.exit(arguments["r1"] + " does not exist. Input file required. Use option -input or -1 / -2.")

		if arguments["r2"]:

			if not os.path.isfile(arguments["r2"]):

				sys.exit(arguments["r2"] + " does not exist.")

	if arguments["trimBetter_threshold"] and not arguments["trimBetter"]:

		sys.exit("--trimBetter must be set to use --trimbetter_threshold. Add --trimBetter to your command or remove --trimbetter_threshold.")

if __name__ == "__main__":

	config = configparser.ConfigParser()

	config.read(join(dirname(abspath(__file__)), "config", "config.txt"))

	kraken_db = ""

	if "DEFAULT" in config:

		if "kraken_db" in config["DEFAULT"].keys():

			kraken_db = config["DEFAULT"]["kraken_db"]

	parser = argparse.ArgumentParser(  formatter_class=argparse.RawDescriptionHelpFormatter,description= """\

	--------------------------------------------------------------------

			    	         < QCumber >

		      < Quality control and read trimming of NGS data >

		   https://gitlab.com/RKIBioinformaticsPipelines/QCumber/

	--------------------------------------------------------------------""",

	epilog = "Example usage: QCumber-2 --input fastq_folder --reference reference.fasta")

	# ------------------------------------------------< INPUT >---------------------------------------------------------#

	group_input = parser.add_argument_group("Input")

	group_input.add_argument('--input', '-i', dest='input',

	                         help="input sample folder. Illumina filenames should end with _<lane>_<R1|R2>_number, e.g. Sample_12_345_R1_001.fastq, to find the right paired set.",

	                         required=False, nargs="+")

	group_input.add_argument('--read1', '-1', dest='r1', help="Read 1 file", required=False)

	group_input.add_argument('--read2', '-2', dest='r2', help="Read 2 file", required=False)

	group_input.add_argument('--technology', '-T', dest='technology', choices=['Illumina', "IonTorrent"],

	                         required=False,

	                         help="If not set, automatically determine technology and search for fastq and bam files. Set technology to IonTorrent if all files are bam-files, else set technology to Illumina.")

	group_input.add_argument('--adapter', '-a', dest='adapter',

	                         choices=['TruSeq2-PE', 'TruSeq2-SE', 'TruSeq3-PE', 'TruSeq3-SE', 'TruSeq3-PE-2',

	                                  'NexteraPE-PE'], help = "Adapter name for trimming. Default: all")

	mapping_exclusion = group_input.add_mutually_exclusive_group()

	mapping_exclusion.add_argument('--reference', '-r', dest='reference', required=False, help="Map reads against reference. Reference needs to be in fasta-format.")

	mapping_exclusion.add_argument('--index', '-I', dest='index', required=False, help="Bowtie2 index if available.")

	group_input.add_argument('--kraken_db', '-d', dest='kraken_db', help='Kraken database', required=False,

	                         default=kraken_db)

	group_optional = parser.add_argument_group("Optional steps")

	group_optional.add_argument('--sav', '-w', dest='sav', required=False,

	                            help="Illumina folder for SAV. Requires RunInfo.xml, RunParameter.xml and Interop folder.")

	group_optional.add_argument('--trimBetter', choices=["assembly", "mapping", "default"],

	                            help="Optimize trimming parameter using 'Per sequence base content' from fastqc. Not recommended for amplicons.")

	group_optional.add_argument('--nokraken', '-K', action="store_true")

	#group_optional.add_argument('--nomapping', '-M', action="store_true")

	group_optional.add_argument('--notrimming', '-Q', action="store_true")

	group_params = parser.add_argument_group("Parameter settings")

	group_params.add_argument('--illuminaclip', '-L', dest='illuminaclip', default="2:30:10",

	                          help='Illuminaclip option: <leading quality>:<trailing quality>:<sliding window>. Default: %(default)s')

	group_params.add_argument('--minlen', '-m', default=50, dest='minlen',

	                          help='Minlen parameter for Trimmomatic. Drops read short than minlen. Default: %(default)s', type=int)

	group_params.add_argument('--trimOption', '-O', dest="trimOption",

	                          help='Additional Trimmomatic input. Default (if trimBetter is not set): SLIDINGWINDOW:4:20',

	                          type=str)

	group_params.add_argument('--trimBetter_threshold', '-b', dest='trimBetter_threshold',

	                          help="Set -trimBetter to use this option.Default setting for Illumina: 0.15 and for IonTorrent: 0.25.",

	                          required=False, type=float)

	group_output = parser.add_argument_group("Output")

	group_output.add_argument('--output', '-o', dest='output', default="")

	group_output.add_argument('--rename', '-R', dest="rename", required=False,

	                          help="TSV File with two columns: <old sample name> <new sample name>")

	group_output.add_argument('--save_mapping', '-S', action="store_true", default=False)

	parser.add_argument('--threads', '-t', dest='threads', default=4, type=int,

	                    help="Number of threads. Default:  %(default)s")

	parser.add_argument('--config', '-c', dest='config', help="Configfile to run pipeline. Additional parameters in the commandline will override arguments in configfile.")

	parser.add_argument('--version', '-v', action='version', version='%(prog)s v' + __version__)

	arguments, unknown_args = parser.parse_known_args()

	arguments  = vars(arguments)

	if arguments["config"]:

		config_args = yaml.load(open(arguments["config"],"r"))

		add_args = [x.replace("-","") for x in sys.argv if x.startswith("-")]

		add_args.extend(list(set(arguments.keys()) - set(config_args.keys()) )  )

		new_args = dict((key, arguments[key]) for key in add_args)

		config_args.update(new_args)

		arguments = config_args

	arguments["output"] = os.path.abspath(arguments["output"])

	if arguments["reference"]:

		arguments["reference"] = os.path.abspath(arguments["reference"])

	if arguments["sav"]:

		arguments["sav"] = os.path.abspath(arguments["sav"])

	if arguments["rename"]:

		arguments["rename"] = os.path.abspath(arguments["rename"])

	parameter = configparser.ConfigParser(strict=False, interpolation=None)

	parameter.read(join(dirname(abspath(__file__)), "config", "parameter.txt"))

	check_input_validity()

	base_pattern = r"(?P<setname>.*)_(?P<lane>L\d{3})_(?P<read>(R1|R2))_(?P<nfiles>\d{3}).*"

	type, samples, joined_samples, name_dict, join_reads = get_input()

	get_defaults()

	os.makedirs(arguments["output"], exist_ok=True)

	yaml.dump(arguments, open(os.path.join(arguments["output"], "config.yaml"), "w"), default_flow_style=False)

	# additional infos

	general_information = OrderedDict()

	general_information["User"] = getpass.getuser()

	general_information["QCumber"] = __version__

	general_information["QCumber_path"] = os.path.dirname(os.path.abspath(__file__))

	general_information["Execution time"] = datetime.datetime.now().ctime()

	system_info = os.uname()

	general_information["Operating system"] = OrderedDict()

	general_information["Operating system"]["System"] = system_info.sysname

	general_information["Operating system"]["Server"] = system_info.nodename

	general_information["Operating system"]["Operating version"] = system_info.version

	general_information["Operating system"]["Release"] = system_info.release

	general_information["Operating system"]["Machine"] = system_info.machine

	general_information["Tool versions"] = OrderedDict()

	general_information["Tool versions"]["Python"] = re.sub("\n", "", sys.version)

	general_information["Tool versions"]["Snakemake"] = snakemake.__version__

	general_information["Tool versions"]["FastQC"] = get_version("fastqc --version")

	if not arguments["notrimming"]:

		general_information["Tool versions"]["Trimmomatic"] = get_version("trimmomatic -version", "trimmomatic")

		if arguments["technology"] == "Illumina" and not arguments["adapter"]:

			general_information["adapter"] = os.path.join(os.path.dirname(os.path.abspath(__file__)), "config",

			                                              "adapters.fa")

		elif arguments["technology"] == "Illumina" and arguments["adapter"]:

			general_information["adapter"] = os.path.join(config["DEFAULT"]["adapter"], arguments["adapter"] + ".fa")

	if arguments["reference"] or arguments["index"]:

		general_information["Tool versions"]["Bowtie2"] = get_version("bowtie2 --version")

	if not arguments["nokraken"]:

		general_information["Tool versions"]["Kraken"] = get_version("kraken --version")

	general_information["Sample information"] = OrderedDict()

	general_information["Sample information"]["type"] = type

	general_information["Sample information"]["samples"] = samples

	general_information["Sample information"]["join_reads"] = join_reads

	general_information["Sample information"]["join_lanes"] = joined_samples

	general_information["Sample information"]["rename"] = name_dict

	os.makedirs(os.path.join(arguments["output"], "QCResults", "_data"), exist_ok=True)

	general_information_file = os.path.join(arguments["output"], "QCResults", "_data",

	                                        "general_information.json")

	json.dump(general_information, open(general_information_file, "w"))

	process = subprocess.Popen(

		"snakemake --snakefile {snakefile} {additional_commands} --directory {workdir} --configfile {configfile} --cores {cores} {targets} ".format(

			additional_commands = " ".join(unknown_args),

			snakefile=os.path.join(os.path.dirname(__file__), "Snakefile"),

			workdir=arguments["output"],

			configfile=general_information_file,

			cores=arguments["threads"],

			targets="QCResults/batch_report.html"),

		shell=True)

	process.communicate()

#!/bin/bash

#conda-env update -f environment/packages.yaml

source activate qcumber

python3 /home/lieuv/Pipelines/QCumber/QCumber-2 $*

source deactivate qcumber