Mercurial > repos > mvdbeek > damid_deseq2_to_peaks

comparison damid_deseq2_to_peaks.py @ 0:3fd7995da4fd draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damid_deseq2_to_peaks commit f37f4b741fd81f663d10523e1636039578c5bb55
author mvdbeek
date Mon, 07 Jan 2019 12:58:55 -0500
parents
children edca422b6cd6
comparison
equal deleted inserted replaced
-1:000000000000 0:3fd7995da4fd
1 import click
2 import pandas as pd
3 import numpy as np
4
5
6 def order_index(df):
7 """
8 Split chr_start_stop in df index and order by chrom and start.
9 """
10 idx = df.index.str.split('_')
11 idx = pd.DataFrame.from_records(list(idx))
12
13 idx.columns = ['chr', 'start', 'stop']
14 idx = idx.astype(dtype={"chr": "object",
15 "start": "int32",
16 "stop": "int32"})
17 coordinates = idx.sort_values(['chr', 'start'])
18 df.index = np.arange(len(df.index))
19 df = df.loc[coordinates.index]
20 df = coordinates.join(df)
21 # index is center of GATC site
22 df.index = df['start'] + 2
23 return df
24
25
26 def significant_gatcs_to_peaks(df, p_value_cutoff):
27 # Add `pass` column for sig. GATCs
28 df['pass'] = 0
29 df.loc[df[6] < p_value_cutoff, 'pass'] = 1
30 # Create pass_id column for consecutive pass or no-pass GATCs
31 # True whenever there is a value change (from previous value):
32 df['pass_id'] = df.groupby('chr')['pass'].diff().ne(0).cumsum()
33 gb = df.groupby('pass_id')
34 # aggregate
35 consecutive_gatcs = gb.aggregate({'chr': np.min, 'start': np.min, 'stop': np.max, 'pass': np.max})
36 # keep only groups with 2 or more GATCS
37 s = gb.size() > 1
38 consecutive_only = consecutive_gatcs[s]
39 # drop GATC groups that are not significant
40 peaks = consecutive_only[consecutive_only['pass'] == 1][['chr', 'start', 'stop']]
41 # calculate region that is not covered.
42 no_peaks = consecutive_only[consecutive_only['pass'] == 0][['chr', 'start', 'stop']]
43 s = no_peaks['stop'] - no_peaks['start']
44 print("%s nt not covered by peaks" % s.sum())
45 s = peaks['stop'] - peaks['start']
46 print("%s nt covered by peaks" % s.sum())
47 return peaks
48
49
50 @click.command()
51 @click.argument('input_path', type=click.Path(exists=True))
52 @click.argument('output_path', type=click.Path())
53 @click.option('--p_value_cutoff', type=float, default=0.01, help="Minimum adjusted p-value for a significant GATC site")
54 def deseq2_gatc_to_peak(input_path, output_path, p_value_cutoff):
55 df = pd.read_csv(input_path, sep='\t', header=None, index_col=0)
56 df = order_index(df)
57 peaks = significant_gatcs_to_peaks(df, p_value_cutoff)
58 peaks.to_csv(output_path, sep='\t', header=None, index=None)
59
60
61 if __name__ == '__main__':
62 deseq2_gatc_to_peak()