<reponame>PavriLab/classifyIS-nf
#!/usr/bin/env python

import pandas as pd
import argparse as ap
import numpy as np
import logging

if __name__ == '__main__':
    logging.basicConfig(format='%(asctime)s - %(message)s', level=logging.INFO)
    parser = ap.ArgumentParser(description='''converts readcounts of IS to log2(RPM) and assigns
                                              each peak a class according to the respective signal
                                              intensities in each condition''')
    parser.add_argument('-ct', '--countTable', required = True,
                        help = '''tab-separated table holding the readquantification information
                                  as generated by deeptools multiBamSummary''')
    parser.add_argument('-b', '--bed', required = True,
                        help = 'BEDfile used to generate counttable')
    parser.add_argument('-o', '--outputFile', required = True,
                        help = 'name of the output file')
    parser.add_argument('-FC', '--foldChange', default = 0.585, type = float,
                        help = 'foldchange threshold for differential regulation determination in log2(RPM)')
    parser.add_argument('-wt', '--wtcolumn', required = True,
                        help = 'name of the column in --counttable holding the WT counts ')
    parser.add_argument('-kd', '--kdcolumn', required = True,
                        help = 'name of the column in --counttable holding the KD counts ')
    parser.add_argument('-t', '--thresholdDA', default = 2, type = float,
                        help = 'cutoff to use to define dormant/absent origins in log2(RPM)')
    args = parser.parse_args()

    logging.info('reading counts %s' % args.countTable)
    counts = pd.read_csv(args.countTable, sep = '\t', header = None, skiprows = 1,
                         names = ['chr', 'start', 'end', args.wtcolumn, args.kdcolumn])
    counts.sort_values(by=['chr', 'start'], inplace = True)
    counts.reset_index(drop = True, inplace = True)

    bed = pd.read_csv(args.bed, sep = '\t', header = None, usecols = [0, 1, 2, 3],
                      names = ['chr', 'start', 'end', 'name'])

    counts = counts.merge(bed, on = ['chr', 'start', 'end'], how = 'left')

    logging.info('computing log2(RPM) values')
    for col in [args.wtcolumn, args.kdcolumn]:
        counts.loc[counts[col] == 0, col] = 1
        counts.loc[:, col] = np.log2((counts[col]*1000000)/counts[col].sum())

    logging.info('assigning regulation class')
    counts['class'] = 0
    fc = args.foldChange

    # dormant
    counts.loc[counts[args.wtcolumn] < counts[counts[args.kdcolumn] < args.thresholdDA][args.wtcolumn].min(), 'class'] = 1
    # absent
    counts.loc[counts[args.kdcolumn] < counts[counts[args.wtcolumn] < args.thresholdDA][args.kdcolumn].min(), 'class'] = 5
    # upregulated
    counts.loc[(counts['class'] == 0) & (counts[args.kdcolumn] > fc + counts[args.wtcolumn]), 'class'] = 2
    # downregulated
    counts.loc[(counts['class'] == 0) & (counts[args.wtcolumn] > fc + counts[args.kdcolumn]), 'class'] = 4
    # unchanged
    counts.loc[counts['class'] == 0, 'class'] = 3

    counts = counts.loc[:, ['chr', 'start', 'end', 'name', args.wtcolumn, args.kdcolumn, 'class']]
    counts.to_csv(args.outputFile, sep = '\t', index = False)

#!/usr/bin/env python

import pandas as pd
import argparse as ap
import numpy as np
import logging

if __name__ == '__main__':
    logging.basicConfig(format='%(asctime)s - %(message)s', level=logging.INFO)
    parser = ap.ArgumentParser(description='''converts readcounts of IS to log2(RPM) and assigns
                                              each peak a class according to the respective signal
                                              intensities in each condition''')
    parser.add_argument('-ct', '--countTable', required = True,
                        help = '''tab-separated table holding the readquantification information
                                  as generated by deeptools multiBamSummary''')
    parser.add_argument('-b', '--bed', required = True,
                        help = 'BEDfile used to generate counttable')
    parser.add_argument('-o', '--outputFile', required = True,
                        help = 'name of the output file')
    parser.add_argument('-FC', '--foldChange', default = 0.585, type = float,
                        help = 'foldchange threshold for differential regulation determination in log2(RPM)')
    parser.add_argument('-wt', '--wtcolumn', required = True,
                        help = 'name of the column in --counttable holding the WT counts ')
    parser.add_argument('-kd', '--kdcolumn', required = True,
                        help = 'name of the column in --counttable holding the KD counts ')
    parser.add_argument('-t', '--thresholdDA', default = 2, type = float,
                        help = 'cutoff to use to define dormant/absent origins in log2(RPM)')
    args = parser.parse_args()

    logging.info('reading counts %s' % args.countTable)
    counts = pd.read_csv(args.countTable, sep = '\t', header = None, skiprows = 1,
                         names = ['chr', 'start', 'end', args.wtcolumn, args.kdcolumn])
    counts.sort_values(by=['chr', 'start'], inplace = True)
    counts.reset_index(drop = True, inplace = True)

    bed = pd.read_csv(args.bed, sep = '\t', header = None, usecols = [0, 1, 2, 3],
                      names = ['chr', 'start', 'end', 'name'])

    counts = counts.merge(bed, on = ['chr', 'start', 'end'], how = 'left')

    logging.info('computing log2(RPM) values')
    for col in [args.wtcolumn, args.kdcolumn]:
        counts.loc[counts[col] == 0, col] = 1
        counts.loc[:, col] = np.log2((counts[col]*1000000)/counts[col].sum())

    logging.info('assigning regulation class')
    counts['class'] = 0
    fc = args.foldChange

    # dormant
    counts.loc[counts[args.wtcolumn] < counts[counts[args.kdcolumn] < args.thresholdDA][args.wtcolumn].min(), 'class'] = 1
    # absent
    counts.loc[counts[args.kdcolumn] < counts[counts[args.wtcolumn] < args.thresholdDA][args.kdcolumn].min(), 'class'] = 5
    # upregulated
    counts.loc[(counts['class'] == 0) & (counts[args.kdcolumn] > fc + counts[args.wtcolumn]), 'class'] = 2
    # downregulated
    counts.loc[(counts['class'] == 0) & (counts[args.wtcolumn] > fc + counts[args.kdcolumn]), 'class'] = 4
    # unchanged
    counts.loc[counts['class'] == 0, 'class'] = 3

    counts = counts.loc[:, ['chr', 'start', 'end', 'name', args.wtcolumn, args.kdcolumn, 'class']]
    counts.to_csv(args.outputFile, sep = '\t', index = False)
#!/usr/bin/env python converts readcounts of IS to log2(RPM) and assigns                                               each peak a class according to the respective signal                                               intensities in each condition tab-separated table holding the readquantification information                                   as generated by deeptools multiBamSummary # dormant # absent # upregulated # downregulated # unchanged