pythonmisc/AssessMappingQuality.py at master · cacampbell/pythonmisc · GitHub

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
#!/usr/bin/env python3
from PairedEndCommand import PairedEndCommand


class Assesser(PairedEndCommand):
    def __init__(self, *args, **kwargs):
        self.input_regex = ".*"
        self.read_regex = ".*"
        self.extension = r".bam$"
        super(Assesser, self).__init__(*args, **kwargs)
        self.set_default("mode", None)
        self.set_default("reference", "reference.fa")
        self.set_default("ref_flat", "ref_flat")
        self.set_default("tmp_dir", "~/tmp")
        self.set_default("picard", "~/.prog/picard-tools-2.5.0/picard.jar")

    def make_command(self, read):
        # Call CollectAlignmentSummaryMetrics, CollectQualityYieldMetrics, and
        # CollectWGSMetrics if the input mode is not RNA
        output = self.rebase_file(read)
        alnsum = self.replace_extension_with(".aln_s.txt", output)
        qualyield = self.replace_extension_with(".qual_yld.txt", output)
        wgsmetrics = self.replace_extension_with(".wgs_met.txt", output)
        rnamet = self.replace_extension_with(".rna_met.txt")
        command = ("java -Xms{xms} -Xmx{xmx} -Djava.io.tmpdir={tmp} -jar "
                   "{picard} CollectAlignmentSummaryMetrics R={ref} I={input} "
                   "O={alnsum} && java -Xms{xms} -Xmx{xmx} "
                   "-Djava.io.tmpdir={tmp} -jar {picard} "
                   "CollectQualityYieldMetrics I={input} O={qualyield}").format(
            xms=self.get_mem(0.95),
            xmx=self.get_mem(0.98),
            tmp=self.tmp_dir,
            picard=self.picard,
            ref=self.reference,
            input=read,
            alnsum=alnsum,
            qualyield=qualyield
        )

        if self.mode.upper().strip() in ["DNA", "RNA"]:
            command += (" && java -Xms{xms} -Xmx{xmx} -Djava.io.tmpdir={tmp} "
                        "-jar {picard} ").format(
                xms=self.get_mem(0.95),
                xmx=self.get_mem(0.98),
                tmp=self.tmp_dir,
                picard=self.picard
            )
            if self.mode.upper().strip() == "DNA":
                command += (" CollectWgsMetrics I={input} O={wgsmet} R={ref} "
                            "INCLUDE_BQ_HISTOGRAM").format(
                    input=read,
                    wgsmet=wgsmetrics,
                    ref=self.reference
                )
            elif self.mode.upper().strip() == "RNA":
                command += (" CollectRnaSeqMetrics I={input} O={rnamet} "
                            "REF_FLAT={refflat}").format(
                    input=read,
                    rnamet=rnamet,
                    refflat=self.ref_flat
                )

        return (command)