7117302237

Committed 06 Dec 2023 04:23PM UTC coverage: 75.482% (+1.8%) from 73.729%

Build # 7117302237

Build Type

push

github

Committed by

cokelaer

Commit Message

Update version

Run Details

13709 of 18162 relevant lines covered (75.48%)

2.26 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

84.09

/sequana/scripts/main/mapping.py

#
#  This file is part of Sequana software
#
#  Copyright (c) 2016-2022 - Sequana Development Team
#
#  Distributed under the terms of the 3-clause BSD license.
#  The full license is in the LICENSE file, distributed with this software.
#
#  website: https://github.com/sequana/sequana
#  documentation: http://sequana.readthedocs.io
#
##############################################################################
import os

import colorlog
import rich_click as click
from snakemake import shell as shellcmd

from sequana.fasta import FastA
from sequana.fastq import FastQ
from sequana.scripts.utils import CONTEXT_SETTINGS

logger = colorlog.getLogger(__name__)


@click.command(context_settings=CONTEXT_SETTINGS)
@click.option("-1", "--file1", show_default=True, required=True, help="R1 Fastq file ; zipped file expected")
@click.option("-2", "--file2", show_default=True, required=False, help="R2 Fastq file")
@click.option("-r", "--reference", show_default=True, required=True, help="Reference where to map data")
@click.option("-t", "--threads", show_default=True, default=1, help="number of threads to use")
@click.option("-p", "--pacbio", show_default=True, default=False, help="", is_flag=True)
@click.option("-s", "--use-sambamba", is_flag=True, help="use sambamba instad of samtools for sorting")
def mapping(**kwargs):
    """map FastQ data onto a reference

This is a simple mapper for quick test. The commands are as follows:

# Indexing
bwa index REFERENCE\
samtools faidx REFERENCE

# mapping
bwa mem -t 4 -R @RG\\tID:1\\tSM:1\\tPL:illumina -T 30 REFERENCE FASTQ_FILES  | samtools
view -Sbh -> REFERENCE.bam

samtools sort -o REFERENCE.sorted.bam  REFERENCE.bam
    """
    reference = kwargs["reference"]
    file1 = kwargs["file1"]
    file2 = kwargs["file2"]
    threads = kwargs["threads"]

    if file1 and file2:
        fastq = f"{file1} {file2}"
    elif file1 and not file2:
        fastq = f"{file1}"
    elif file1 is None:
        raise ValueError("--file1 must be used")

    S = sum([len(this["sequence"]) for this in FastQ(file1)])

    if file2:
        S += sum([len(this["sequence"]) for this in FastQ(file2)])

    ref = FastA(reference)
    coverage = float(S) / len(ref.sequences[0])
    print("Theoretical Depth of Coverage : %s" % coverage)

    # indexing
    shellcmd(f"bwa index {reference} ")
    cmd = f"samtools faidx {reference} "

    # mapping
    cmd = "bwa mem -M "  # mark shorter split read as secondary; -M is not compulsary but recommended
    if kwargs["pacbio"]:
        cmd += "-x pacbio "
    cmd += f" -t {threads} " + r"-R @RG\\tID:1\\tSM:1\\tPL:illumina -T 30 " + f"{reference} {fastq}  "

    # Samtools options:
    #   S:ignore input format
    #   h:include header
    #   b:bam output
    if kwargs["use_sambamba"] is False:
        cmd += "| samtools view -Sbh | "
        # sorting BAM
        cmd += f"samtools sort -@ {threads} -o {reference}.sorted.bam -"
        shellcmd(cmd)
    else:
        # FIXME use sambamba for the view as well
        cmd += f"| samtools view -Sbu - | sambamba sort /dev/stdin -o {reference}.sorted.bam -t {threads}  --tmpdir=./tmp  "
        shellcmd(cmd)

1	#
2	# This file is part of Sequana software
3	#
4	# Copyright (c) 2016-2022 - Sequana Development Team
5	#
6	# Distributed under the terms of the 3-clause BSD license.
7	# The full license is in the LICENSE file, distributed with this software.
8	#
9	# website: https://github.com/sequana/sequana
10	# documentation: http://sequana.readthedocs.io
11	#
12	##############################################################################
13	import os	3✔
14
15	import colorlog	3✔
16	import rich_click as click	3✔
17	from snakemake import shell as shellcmd	3✔
18
19	from sequana.fasta import FastA	3✔
20	from sequana.fastq import FastQ	3✔
21	from sequana.scripts.utils import CONTEXT_SETTINGS	3✔
22
23	logger = colorlog.getLogger(__name__)	3✔
24
25
26	@click.command(context_settings=CONTEXT_SETTINGS)	3✔
27	@click.option("-1", "--file1", show_default=True, required=True, help="R1 Fastq file ; zipped file expected")	3✔
28	@click.option("-2", "--file2", show_default=True, required=False, help="R2 Fastq file")	3✔
29	@click.option("-r", "--reference", show_default=True, required=True, help="Reference where to map data")	3✔
30	@click.option("-t", "--threads", show_default=True, default=1, help="number of threads to use")	3✔
31	@click.option("-p", "--pacbio", show_default=True, default=False, help="", is_flag=True)	3✔
32	@click.option("-s", "--use-sambamba", is_flag=True, help="use sambamba instad of samtools for sorting")	3✔
33	def mapping(**kwargs):	3✔
34	"""map FastQ data onto a reference
35
36	This is a simple mapper for quick test. The commands are as follows:
37
38	# Indexing
39	bwa index REFERENCE\
40	samtools faidx REFERENCE
41
42	# mapping
43	bwa mem -t 4 -R @RG\\tID:1\\tSM:1\\tPL:illumina -T 30 REFERENCE FASTQ_FILES \| samtools
44	view -Sbh -> REFERENCE.bam
45
46	samtools sort -o REFERENCE.sorted.bam REFERENCE.bam
47	"""
48	reference = kwargs["reference"]	3✔
49	file1 = kwargs["file1"]	3✔
50	file2 = kwargs["file2"]	3✔
51	threads = kwargs["threads"]	3✔
52
53	if file1 and file2:	3✔
54	fastq = f"{file1} {file2}"	3✔
55	elif file1 and not file2:	×
56	fastq = f"{file1}"	×
57	elif file1 is None:	×
58	raise ValueError("--file1 must be used")	×
59
60	S = sum([len(this["sequence"]) for this in FastQ(file1)])	3✔
61
62	if file2:	3✔
63	S += sum([len(this["sequence"]) for this in FastQ(file2)])	3✔
64
65	ref = FastA(reference)	3✔
66	coverage = float(S) / len(ref.sequences[0])	3✔
67	print("Theoretical Depth of Coverage : %s" % coverage)	3✔
68
69	# indexing
70	shellcmd(f"bwa index {reference} ")	3✔
71	cmd = f"samtools faidx {reference} "	3✔
72
73	# mapping
74	cmd = "bwa mem -M " # mark shorter split read as secondary; -M is not compulsary but recommended	3✔
75	if kwargs["pacbio"]:	3✔
76	cmd += "-x pacbio "	×
77	cmd += f" -t {threads} " + r"-R @RG\\tID:1\\tSM:1\\tPL:illumina -T 30 " + f"{reference} {fastq} "	3✔
78
79	# Samtools options:
80	# S:ignore input format
81	# h:include header
82	# b:bam output
83	if kwargs["use_sambamba"] is False:	3✔
84	cmd += "\| samtools view -Sbh \| "	3✔
85	# sorting BAM
86	cmd += f"samtools sort -@ {threads} -o {reference}.sorted.bam -"	3✔
87	shellcmd(cmd)	3✔
88	else:
89	# FIXME use sambamba for the view as well
90	cmd += f"\| samtools view -Sbu - \| sambamba sort /dev/stdin -o {reference}.sorted.bam -t {threads} --tmpdir=./tmp "	×
91	shellcmd(cmd)	×

cokelaer / sequana / 7117302237

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous