7093643877

Committed 04 Dec 2023 11:00PM UTC coverage: 73.729% (+0.3%) from 73.466%

Build # 7093643877

Build Type

push

github

Committed by

cokelaer

Commit Message

rename plot legend

Run Details

0 of 1 new or added line in 1 file covered. (0.0%)

271 existing lines in 19 files now uncovered.

13398 of 18172 relevant lines covered (73.73%)

2.21 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

from snakemake import shell as shellcmd

import rich_click as click
import colorlog

from sequana.fastq import FastQ
from sequana.fasta import FastA

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

cokelaer / sequana / 7093643877

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