Make the read GC plot for Fastp
499928
def fastp_read_gc_plot(self):
        
        data_labels, pdata = self.filter_pconfig_pdata_subplots(self.fastp_gc_content_data, 'Base Content Percent')
        pconfig = {
            'id': 'fastp-seq-content-gc-plot',
            'title': 'Fastp: Read GC Content',
            'xlab': 'Read Position',
            'ylab': 'R1 Before filtering: Base Content Percent',
            'ymax': 100,
            'ymin': 0,
            'xDecimals': False,
            'yLabelFormat': '{value}%',
            'tt_label': '{point.x}: {point.y:.2f}%',
            'data_labels': data_labels
        }
        return linegraph.plot(pdata, pconfig)
1243616
def gp_rapp():
  """"""
  inDir, outDir = getWorkDirs()
  # prepare data
  yields = {}
  for infile in os.listdir(inDir):
    energy = re.compile('\d+').search(infile).group()
    medium = np.loadtxt(open(os.path.join(inDir, infile), 'rb'))
    getMassRangesSums(energy, medium, yields)
  data = dict( # sort by energy
    (k, np.array(sorted(v)))
    for k, v in yields.iteritems()
  )
  for k in data: data[k][:,1] /= data[k][-1,1] # divide by 200
  # make plot
  nSets = len(data)
  make_plot(
    data = [ data[k][:-1] for k in data ],
    properties = [
      'with linespoints lt 1 lw 4 ps 1.5 lc %s pt 18' % default_colors[i]
      for i in xrange(nSets)
    ],
    titles = data.keys(), # TODO: titles order correct?
    name = os.path.join(outDir, 'rapp'),
    xlabel = '{/Symbol \326}s_{NN} (GeV)', ylabel = 'Rapp Ratio to 200 GeV',
    lmargin = 0.1, key = ['left'], yr = [0.1, 0.8]
  )
  return 'done'
677722
def gb(args):
    
    from Bio.Alphabet import generic_dna
    try:
        from BCBio import GFF
    except ImportError:
        print("You need to install dep first: $ easy_install bcbio-gff", file=sys.stderr)

    p = OptionParser(gb.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    gff_file, fasta_file = args
    pf = op.splitext(gff_file)[0]
    out_file = pf + ".gb"
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")
1366056
def gcn(args):
    
    p = OptionParser(gcn.__doc__)
    p.set_cpus()
    p.set_tmpdir(tmpdir="tmp")
    p.set_outfile()
    opts, args = p.parse_args(args)

    if len(args) < 2:
        sys.exit(not p.print_help())

    exonbed = args[0]
    canvasvcfs = args[1:]
    tsvfile = opts.outfile
    tmpdir = opts.tmpdir

    mkdir(tmpdir)
    set_tempdir(tmpdir)

    df = vcf_to_df(canvasvcfs, exonbed, opts.cpus)
    for suffix in (".avgcn", ".medcn"):
        df_to_tsv(df, tsvfile, suffix)
682287
def parse_gc_dist(self, f):
    
    # Get the sample name from the parent parent directory
    # Typical path: <sample name>/raw_data_qualimapReport/mapped_reads_gc-content_distribution.txt
    s_name = self.get_s_name(f)

    d = dict()
    reference_species = None
    reference_d = dict()
    avg_gc = 0
    for l in f['f']:
        if l.startswith('#'):
            sections = l.strip("\n").split("\t", 3)
            if len(sections) > 2:
                reference_species = sections[2]
            continue
        sections = l.strip("\n").split("\t", 3)
        gc = int(round(float(sections[0])))
        content = float(sections[1])
        avg_gc += gc * content
        d[gc] = content
        if len(sections) > 2:
            reference_content = float(sections[2])
            reference_d[gc] = reference_content

    # Add average GC to the general stats table
    self.general_stats_data[s_name]['avg_gc'] = avg_gc

    # Save results
    if s_name in self.qualimap_bamqc_gc_content_dist:
        log.debug("Duplicate Mapped Reads GC content distribution sample name found! Overwriting: {}".format(s_name))
    self.qualimap_bamqc_gc_content_dist[s_name] = d
    if reference_species and reference_species not in self.qualimap_bamqc_gc_by_species:
        self.qualimap_bamqc_gc_by_species[reference_species] = reference_d
    self.add_data_source(f, s_name=s_name, section='mapped_gc_distribution')
1348160
def wgsim(args):
    
    p = OptionParser(wgsim.__doc__)
    p.add_option("--erate", default=.01, type="float",
                 help="Base error rate of the read [default: %default]")
    p.add_option("--noerrors", default=False, action="store_true",
                 help="Simulate reads with no errors [default: %default]")
    p.add_option("--genomesize", type="int",
                 help="Genome size in Mb [default: estimate from data]")
    add_sim_options(p)
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    fastafile, = args
    pf = op.basename(fastafile).split(".")[0]

    genomesize = opts.genomesize
    size = genomesize * 1000000 if genomesize else Fasta(fastafile).totalsize
    depth = opts.depth
    readlen = opts.readlen
    readnum = int(math.ceil(size * depth / (2 * readlen)))

    distance = opts.distance
    stdev = distance / 10

    outpf = opts.outfile or "{0}.{1}bp.{2}x".format(pf, distance, depth)

    logging.debug("Total genome size: {0} bp".format(size))
    logging.debug("Target depth: {0}x".format(depth))
    logging.debug("Number of read pairs (2x{0}): {1}".format(readlen, readnum))

    if opts.noerrors:
        opts.erate = 0

    cmd = "dwgsim -e {0} -E {0}".format(opts.erate)
    if opts.noerrors:
        cmd += " -r 0 -R 0 -X 0 -y 0"

    cmd += " -d {0} -s {1}".format(distance, stdev)
    cmd += " -N {0} -1 {1} -2 {1}".format(readnum, readlen)
    cmd += " {0} {1}".format(fastafile, outpf)
    sh(cmd)
16464
def read_gcvs(filename):
    
    with open(filename, 'r') as fp:
        parser = GcvsParser(fp)
        for star in parser:
            yield star
571152
def parse_fastp_log(self, f):
        
        try:
            parsed_json = json.load(f['f'])
        except:
            log.warn("Could not parse fastp JSON: '{}'".format(f['fn']))
            return None

        # Fetch a sample name from the command
        s_name = f['s_name']
        cmd = parsed_json['command'].split()
        for i,v in enumerate(cmd):
            if v == '-i':
                s_name = self.clean_s_name(cmd[i+1], f['root'])
        if s_name == 'fastp':
            log.warn('Could not parse sample name from fastp command: {}'.format(f['fn']))

        self.add_data_source(f, s_name)
        self.fastp_data[s_name] = {}
        self.fastp_duplication_plotdata[s_name] = {}
        self.fastp_insert_size_data[s_name] = {}
        self.fastp_all_data[s_name] = parsed_json
        for k in ['read1_before_filtering','read2_before_filtering','read1_after_filtering','read2_after_filtering']:
            self.fastp_qual_plotdata[k][s_name] = {}
            self.fastp_gc_content_data[k][s_name] = {}
            self.fastp_n_content_data[k][s_name] = {}

        # Parse filtering_result
        try:
            for k in parsed_json['filtering_result']:
                self.fastp_data[s_name]['filtering_result_{}'.format(k)] = float(parsed_json['filtering_result'][k])
        except KeyError:
            log.debug("fastp JSON did not have 'filtering_result' key: '{}'".format(f['fn']))

        # Parse duplication
        try:
            self.fastp_data[s_name]['pct_duplication'] = float(parsed_json['duplication']['rate'] * 100.0)
        except KeyError:
            log.debug("fastp JSON did not have a 'duplication' key: '{}'".format(f['fn']))

        # Parse after_filtering
        try:
            for k in parsed_json['summary']['after_filtering']:
                self.fastp_data[s_name]['after_filtering_{}'.format(k)] = float(parsed_json['summary']['after_filtering'][k])
        except KeyError:
            log.debug("fastp JSON did not have a 'summary'-'after_filtering' keys: '{}'".format(f['fn']))


        # Parse data required to calculate Pct reads surviving
        try:
            self.fastp_data[s_name]['before_filtering_total_reads'] = float(parsed_json['summary']['before_filtering']['total_reads'])
        except KeyError:
            log.debug("Could not find pre-filtering # reads: '{}'".format(f['fn']))

        try:
            self.fastp_data[s_name]['pct_surviving'] = (self.fastp_data[s_name]['after_filtering_total_reads'] / self.fastp_data[s_name]['before_filtering_total_reads']) * 100.0
        except KeyError:
            log.debug("Could not calculate 'pct_surviving': {}".format(f['fn']))

        # Parse adapter_cutting
        try:
            for k in parsed_json['adapter_cutting']:
                try:
                    self.fastp_data[s_name]['adapter_cutting_{}'.format(k)] = float(parsed_json['adapter_cutting'][k])
                except (ValueError, TypeError):
                    pass
        except KeyError:
            log.debug("fastp JSON did not have a 'adapter_cutting' key, skipping: '{}'".format(f['fn']))

        try:
            self.fastp_data[s_name]['pct_adapter'] = (self.fastp_data[s_name]['adapter_cutting_adapter_trimmed_reads'] / self.fastp_data[s_name]['before_filtering_total_reads']) * 100.0
        except KeyError:
            log.debug("Could not calculate 'pct_adapter': {}".format(f['fn']))

        # Duplication rate plot data
        try:
            # First count the total read count in the dup analysis
            total_reads = 0
            for v in parsed_json['duplication']['histogram']:
                total_reads += v
            if total_reads == 0:
                raise KeyError
            # Calculate percentages
            for i, v in enumerate(parsed_json['duplication']['histogram']):
                self.fastp_duplication_plotdata[s_name][i+1] = (float(v) / float(total_reads)) * 100.0
        except KeyError:
            log.debug("No duplication rate plot data: {}".format(f['fn']))

        # Insert size plot data
        try:
            # First count the total read count in the insert size analysis
            total_reads = 0
            max_i = 0
            for i, v in enumerate(parsed_json['insert_size']['histogram']):
                total_reads += v
                if float(v) > 0:
                    max_i = i
            if total_reads == 0:
                raise KeyError
            # Calculate percentages
            for i, v in enumerate(parsed_json['insert_size']['histogram']):
                if i <= max_i:
                    self.fastp_insert_size_data[s_name][i+1] = (float(v) / float(total_reads)) * 100.0
        except KeyError:
            log.debug("No insert size plot data: {}".format(f['fn']))

        for k in ['read1_before_filtering','read2_before_filtering','read1_after_filtering','read2_after_filtering']:
            # Read quality data
            try:
                for i, v in enumerate(parsed_json[k]['quality_curves']['mean']):
                    self.fastp_qual_plotdata[k][s_name][i+1] = float(v)
            except KeyError:
                log.debug("Read quality {} not found: {}".format(k, f['fn']))

            # GC and N content plots
            try:
                for i, v in enumerate(parsed_json[k]['content_curves']['GC']):
                    self.fastp_gc_content_data[k][s_name][i+1] = float(v) * 100.0
                for i, v in enumerate(parsed_json[k]['content_curves']['N']):
                    self.fastp_n_content_data[k][s_name][i+1] = float(v) * 100.0
            except KeyError:
                log.debug("Content curve data {} not found: {}".format(k, f['fn']))


        # Remove empty dicts
        if len(self.fastp_data[s_name]) == 0:
            del self.fastp_data[s_name]
        if len(self.fastp_duplication_plotdata[s_name]) == 0:
            del self.fastp_duplication_plotdata[s_name]
        if len(self.fastp_insert_size_data[s_name]) == 0:
            del self.fastp_insert_size_data[s_name]
        if len(self.fastp_all_data[s_name]) == 0:
            del self.fastp_all_data[s_name]
349926
def fastqc_general_stats(self):
        

        # Prep the data
        data = dict()
        for s_name in self.fastqc_data:
            bs = self.fastqc_data[s_name]['basic_statistics']
            data[s_name] = {
                'percent_gc': bs['%GC'],
                'avg_sequence_length': bs['avg_sequence_length'],
                'total_sequences': bs['Total Sequences'],
            }
            try:
                data[s_name]['percent_duplicates'] = 100 - bs['total_deduplicated_percentage']
            except KeyError:
                pass # Older versions of FastQC don't have this
            # Add count of fail statuses
            num_statuses = 0
            num_fails = 0
            for s in self.fastqc_data[s_name]['statuses'].values():
                num_statuses += 1
                if s == 'fail':
                    num_fails += 1
            data[s_name]['percent_fails'] = (float(num_fails)/float(num_statuses))*100.0

        # Are sequence lengths interesting?
        seq_lengths = [x['avg_sequence_length'] for x in data.values()]
        hide_seq_length = False if max(seq_lengths) - min(seq_lengths) > 10 else True

        headers = OrderedDict()
        headers['percent_duplicates'] = {
            'title': '% Dups',
            'description': '% Duplicate Reads',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'RdYlGn-rev'
        }
        headers['percent_gc'] = {
            'title': '% GC',
            'description': 'Average % GC Content',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'Set1',
            'format': '{:,.0f}'
        }
        headers['avg_sequence_length'] = {
            'title': 'Length',
            'description': 'Average Sequence Length (bp)',
            'min': 0,
            'suffix': ' bp',
            'scale': 'RdYlGn',
            'format': '{:,.0f}',
            'hidden': hide_seq_length
        }
        headers['percent_fails'] = {
            'title': '% Failed',
            'description': 'Percentage of modules failed in FastQC report (includes those not plotted here)',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'Reds',
            'format': '{:,.0f}',
            'hidden': True
        }
        headers['total_sequences'] = {
            'title': '{} Seqs'.format(config.read_count_prefix),
            'description': 'Total Sequences ({})'.format(config.read_count_desc),
            'min': 0,
            'scale': 'Blues',
            'modify': lambda x: x * config.read_count_multiplier,
            'shared_key': 'read_count'
        }
        self.general_stats_addcols(data, headers)
3718
def gff(args):
    
    align_choices = ("left", "center", "right")
    p = OptionParser(gff.__doc__)
    p.add_option("--align", default="left", choices=align_choices,
                 help="Horizontal alignment [default: %default]")
    p.add_option("--noUTR", default=False, action="store_true",
                 help="Do not plot UTRs [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(not p.print_help())

    fig = plt.figure(1, (8, 5))
    root = fig.add_axes([0, 0, 1, 1])

    gffiles = args
    ngenes = len(gffiles)

    canvas = .6
    setups, ratio = get_setups(gffiles, canvas=canvas, noUTR=opts.noUTR)
    align = opts.align
    xs = .2 if align == "left" else .8
    yinterval = canvas / ngenes
    ys = .8
    tip = .01
    for genename, mrnabed, cdsbeds in setups:
        ExonGlyph(root, xs, ys, mrnabed, cdsbeds, ratio=ratio, align=align)
        if align == "left":
            root.text(xs - tip, ys, genename, ha="right", va="center")
        elif align == "right":
            root.text(xs + tip, ys, genename, ha="left", va="center")
        ys -= yinterval

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    figname = "exons.pdf"
    savefig(figname, dpi=300)
508546
def _add_umis_with_fastp(read_fq, umi_fq, out_fq, cores):
    
    with utils.open_gzipsafe(umi_fq) as in_handle:
        in_handle.readline()  # name
        umi_size = len(in_handle.readline().strip())
    cmd = ("fastp -Q -A -L -G -w 1 --in1 {read_fq} --in2 {umi_fq} "
           "--umi --umi_prefix UMI --umi_loc read2 --umi_len {umi_size} "
           "--out1 >(bgzip --threads {cores} -c > {out_fq}) --out2 /dev/null "
           "-j /dev/null -h /dev/null")
    do.run(cmd.format(**locals()), "Add UMIs to fastq file with fastp")
720076
def read_fastas(input_files):
    
    tumor_file = [y for x, y in input_files.items() if x.startswith('T')][0]
    normal_file = [y for x, y in input_files.items() if x.startswith('N')][0]
    output_files = defaultdict(list)
    output_files = _read_fasta(tumor_file, output_files)
    num_entries = len(output_files)
    output_files = _read_fasta(normal_file, output_files)
    assert len(output_files) == num_entries
    return output_files
296685
def fastp_filtered_reads_chart(self):
        
        # Specify the order of the different possible categories
        keys = OrderedDict()
        keys['filtering_result_passed_filter_reads'] =  { 'name': 'Passed Filter' }
        keys['filtering_result_low_quality_reads'] =    { 'name': 'Low Quality' }
        keys['filtering_result_too_many_N_reads'] =     { 'name': 'Too Many N' }
        keys['filtering_result_too_short_reads'] =      { 'name': 'Too short' }

        # Config for the plot
        pconfig = {
            'id': 'fastp_filtered_reads_plot',
            'title': 'Fastp: Filtered Reads',
            'ylab': '# Reads',
            'cpswitch_counts_label': 'Number of Reads',
            'hide_zero_cats': False,
        }
        return bargraph.plot(self.fastp_data, keys, pconfig)
857845
def plotCurve(self):
        """"""
        try:
            attenuations, freqs = self.datafile.get_calibration(str(self.ui.calChoiceCmbbx.currentText()), self.calf)
            self.pw = SimplePlotWidget(freqs, attenuations, parent=self)
            self.pw.setWindowFlags(QtCore.Qt.Window)
            self.pw.setLabels('Frequency', 'Attenuation', 'Calibration Curve')
            self.pw.show()
        except IOError:
            QtGui.QMessageBox.warning(self, "File Read Error", "Unable to read calibration file")
        except KeyError:
            QtGui.QMessageBox.warning(self, "File Data Error", "Unable to find data in file")
345429
def write_prepareref_fasta_file(outfile, gene_coords, genes_need_upstream, genes_non_upstream, upstream_before=100, upstream_after=100):
    ''''''
    tmp_dict = {}
    fasta_in = os.path.join(data_dir, 'NC_000962.3.fa.gz')
    pyfastaq.tasks.file_to_dict(fasta_in, tmp_dict)
    ref_seq = tmp_dict['NC_000962.3']

    with open(outfile, 'w') as f:
        for gene in genes_non_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start:end+1])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[end:start+1])
                gene_fa.revcomp()

            print(gene_fa, file=f)

        for gene in genes_need_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_before:start + upstream_after])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_after + 1:start + upstream_before + 1])
                gene_fa.revcomp()

            gene_fa.id += '_upstream'
            print(gene_fa, file=f)
747883
def set_plot_CC_T_rho_new(self,fig='CC evol',linestyle=['-'],burn_limit=0.997,color=['r'],marker=['o'],markevery=500): 
		'''
			
		'''
		if len(linestyle)==0:
			linestyle=200*['-']
		plt.figure(fig)
                for i in range(len(self.runs_H5_surf)):
                        sefiles=se(self.runs_H5_out[i])
			t1_model=-1
                        mini=sefiles.get('mini')
			zini=sefiles.get('zini')
			label=str(mini)+'$M_{\odot}$, Z='+str(zini)
                        model=sefiles.se.cycles
                        model_list=[]
                        for k in range(0,len(model),1):
                                model_list.append(model[k])
			print 'REad Rho,T, this might take a while...'
                        rho1=sefiles.get(model_list,'rho')   #[:(t1_model-t0_model)]
                        T1=sefiles.get(model_list,'temperature')#[:(t1_model-t0_model)]
			print 'finished'
                        rho=[]
                        T=[]
			T_unit=sefiles.get('temperature_unit')
			labeldone=False
			for k in range(len(model_list)):
				#print 'test model ',model_list[k]
                                t9max=max(np.array(T1[k])*T_unit/1.e9)
                                #T.append(max(t9))
				rho1max=max(rho1[k])
				print 'model',model_list[k]
				print 'maxT, maxrho'
				print t9max,rho1max
				if k==0:
					t9_prev=t9max
					rho1_prev=rho1max
					idx_T=0
					idx_rho=0
					continue
				if t9max>t9_prev:
					idx_T=k
					t9_prev=t9max
				if rho1max>rho1_prev:
					idx_rho=k
					rho1_prev=rho1max
			print 'found highest rho',idx_rho,max(np.array(T1[idx_rho])*T_unit/1.0e9),max(rho1[idx_rho]),model_list[idx_rho]
			print 'found highest T',idx_T,max(np.array(T1[idx_T])*T_unit/1.0e9),max(rho1[idx_T]),model_list[idx_T]
			if idx_T==idx_rho:
				x=np.array(T1[idx_T])*T_unit/1e9
				y=rho1[idx_T]
				rho1=[]	
				T1=[]
				for k in range(len(x)):
				    if not y[k]==1.0:
			        	rho1.append(y[k])
        				T1.append(x[k])
				x=T1
				y=rho1
				plt.plot(x,y,label=label,color=color[i],marker=marker[i],linestyle=linestyle[i],markevery=markevery)
                                #rhoi.append(max(rho1[k]))	       
			else:
				#for max T
				x=np.array(T1[idx_T])*T_unit/1e9
				y=rho1[idx_T]
                                rho_temp=[]
                                T_temp=[]
                                for k in range(len(x)):
                                    if not y[k]==1.0:
                                        rho_temp.append(y[k])
                                        T_temp.append(x[k])
                                x=T_temp
                                y=rho_temp
				plt.plot(x,y,label=label,color=color[i],marker=marker[i],linestyle=linestyle[i],markevery=markevery)
				#for max rho
        			x=np.array(T1[idx_rho])*T_unit/1e9
				y=rho1[idx_rho]
                                rho_temp=[]
                                T_temp=[]
                                for k in range(len(x)):
                                    if not y[k]==1.0:
                                        rho_temp.append(y[k])
                                        T_temp.append(x[k])
                                x=T_temp
                                y=rho_temp
				plt.plot(x,y,color=color[i],marker=marker[i],linestyle=linestyle[i],markevery=markevery)
29482
def lgdistplot(self,dict_to_use,orientation):
        """"""

        data = dict()
        for s_name in dict_to_use:
            try:
                data[s_name] = {int(d): int (dict_to_use[s_name][d]) for d in dict_to_use[s_name]}
            except KeyError:
                pass
        if len(data) == 0:
            log.debug('No valid data for forward read lgdist input!')
            return None

        config = {
            'id': 'length-distribution-{}'.format(orientation),
            'title': 'DamageProfiler: Read length distribution: {} '.format(orientation),
            'ylab': 'Number of reads',
            'xlab': 'Readlength (bp)',
            'xDecimals': False,
            'tt_label': '{point.y} reads of length {point.x}',
            'ymin': 0,
            'xmin': 0
        }
        return linegraph.plot(data,config)