-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy path6_position_encoding.py
More file actions
executable file
·145 lines (108 loc) · 3.98 KB
/
Copy path6_position_encoding.py
File metadata and controls
executable file
·145 lines (108 loc) · 3.98 KB
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
import os
import re
import argparse
import pandas as pd
import numpy as np
from sklearn.preprocessing import OneHotEncoder
def main(args):
# Load the different species
SPECIES = pd.read_csv('../data/species.txt')
# Load the genome length
genome = pd.read_excel('../data/PATHOgenex_strains_genome_sizes.xlsx')
lessthan15 = np.loadtxt('../data/processed_files/splits/lessThan15genes.txt', dtype=str)
# Extract and save the expression values of interest
# from the raw file
c = 0
for species in SPECIES.values:
#for species in ["Shigella flexneri 5a str. M90T_new.csv"]:
species = species[0]
name = return_name(species)
if (name not in lessthan15):
c+=1
print(species.split('.csv')[0])
# Loading the cleaned expression data
expression_data = pd.read_csv('{}{}.csv'.format(args.inexp, name), index_col = 0)
# we need this for the position on the genome
raw_data = pd.read_csv(args.inraw + species, sep = ';', low_memory=False, index_col = 0)
raw_data = raw_data.set_index(expression_data.index.name)
# obtain polar coordinates
if("_new" in species):
gen_sp = genome[genome['Strain'] == species.split('_new.csv')[0]]
else:
gen_sp = genome[genome['Strain'] == species.split('.csv')[0]]
sub_raw_data = raw_data.loc[expression_data.index]
coord = polar(sub_raw_data, gen_sp)
# Unify the information
all_features = pd.concat([expression_data, coord], axis=1)
all_features.to_csv('{}{}_genomic_info.csv'.format(args.inexp, name))
print(all_features)
print('')
print(c)
def polar(raw_data, genome):
strand, starts, ends = get_info(raw_data)
chromosome = raw_data['Chromosome'].values
check(chromosome, genome)
chrom_info = np.zeros(len(strand)) # this is a binary vector that tells you
# whether the genetic information is located on a chromosome (0) or on a plasmide (1)
for idx, row in genome.iterrows():
gen_id_idx = (raw_data['Chromosome'] == row['Accession']).values
starts[gen_id_idx] = starts[gen_id_idx]/row['Size (bp)']
ends[gen_id_idx] = ends[gen_id_idx]/row['Size (bp)']
if(row['Genome'] == 'Plasmid'):
chrom_info[gen_id_idx] = 1
print(len(np.unique(chrom_info)))
# convert starts and end to radiant
starts_rad = starts*2*np.pi
ends_rad = ends*2*np.pi
data = np.c_[np.cos(starts_rad), np.sin(starts_rad), np.cos(ends_rad), np.sin(ends_rad), strand, chrom_info,]
info_df = pd.DataFrame(data = data, columns = ['cos-start', 'sin-start','cos-end', 'sin-end', 'strand', 'genome'], index = raw_data.index)
print(info_df)
return info_df
def check(chromosome, genome):
'''
This function checks whether we miss the information on
the length of the chromosomes/plasmides in some cases
'''
c = np.unique(chromosome)
g = genome['Accession'].values
inter = np.intersect1d(c, g)
if (len(c) > len(g)):
print('Something wrong!')
def get_info(genome):
region = genome['Region'].values
strand, starts, ends = [], [], []
for r in region:
if('complement' in r):
strand.append(0)
else:
strand.append(1)
parts = r.split('..')
start = re.sub(r'[^0-9]', '', parts[0])
end = re.sub(r'[^0-9]', '', parts[1])
if (start < end):
starts.append(float(start))
ends.append(float(end))
elif(start > end):
starts.append(float(end))
ends.append(float(start))
else:
print('Region of length 0 bps')
return np.array(strand), np.array(starts), np.array(ends)
def return_name(infile):
parts = infile.split('.csv')
name = parts[0]
return name.replace(' ', '_')
def parse_arguments():
'''
Definition of the command line arguments
'''
parser = argparse.ArgumentParser()
parser.add_argument('--inraw', required = False,
help = 'Raw expression file path', default = '../data/OneDrive_1_2-1-2022/Gene expression and DEG files/')
parser.add_argument('--inexp', required = False,
help = 'output dir processed expression data', default = '../data/processed_files/features/')
args = parser.parse_args()
return args
if __name__ == '__main__':
arguments = parse_arguments()
main(arguments)