summaryrefslogtreecommitdiff
path: root/Code/make_parameter_rnaseq.py
blob: 18ef568d1d16100ffd3491c6944c742614183d05 (plain)
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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
# Usage: python make_parameter_rnaseq.py [id-list.txt] > parameter_for_buildRmatrix.txt
#        Edit QUANT_PATH, set NON_ZERO_RATIO.
#
# Purpose: automatically generate parameter_for_buildRmatrix.txt
#
# Update: 26 Feb 2017, slcu, hui
# Update: 19 Sep 2019, slcu, hui [add  read_ena_data_info_json]

import sys, os, glob, json
import fnmatch, re
from datetime import datetime
from configure import RNA_SEQ_INFO_FILE, UPDATE_NETWORK_LOG_FILE

NON_ZERO_RATIO = 0.2 # omit *_quant.txt files with too many zeros.
QUANT_PATH = ['../Data/R/Mapped/public', '../Data/R/Mapped/inhouse', '../Data/R/Mapped/other']  # places where all _quant.txt reside. _quant.txt in sub-directories will also be used.

def get_quant_file_list(fname):
    f = open(fname)
    result = []
    for line in f:
        line = line.strip()
        lst = line.split()
        result.append(lst[0])
    f.close()
    return result

def extract_id(s, src):
    if src == 'inhouse' or src == 'other':
        dirname = os.path.dirname(s)
        lst = dirname.split('/')
        parent_dir = lst[-1]
        libno = filter(str.isdigit, parent_dir) # extrac library number
        libno = libno.zfill(3)
        sample_no = os.path.basename(s) # extract sample number
        first_match = re.findall('_S\d+_', sample_no)[0]
        sample_no = filter(str.isdigit, first_match).zfill(2)
        return '0000' + libno + sample_no
    if src == 'sra':
        index = s.find('_quant')
        if index > 0:
            return s[:index]
        else:
            return 'NA'
    return 'NA'

def zfill2(s, n):
    return s + 'X' * (n-len(s))
    
def make_id(success_flag, lab_id, myid): # should also work for non-SRA id
    if lab_id == 'SRR' or lab_id == 'ERR' or lab_id == 'DRR':
        result = 'R' + success_flag + lab_id + zfill2(myid, 9)
    else: # inhouse or other
        result = 'R' + success_flag + lab_id + myid
    return result
    
def glob_files_include_path(directory, pattern):
    ''' return all file names (with paths) given directory and pattern '''
    result = []
    for root, dirnames, filenames in os.walk(directory):
        for filename in fnmatch.filter(filenames, pattern):
            result.append(os.path.join(root, filename))
    return result

def glob_files_include_paths(directory_lst, pattern):
    ''' return all file names (with paths) given directory and pattern '''
    result = []
    for directory in directory_lst:
        result.extend(glob_files_include_path(os.path.abspath(directory), pattern))
    return result

def non_zero_ratio(fname):
    non_zero_count = 0
    total_count = 0
    f = open(fname)
    lines = f.readlines()
    f.close()
    for line in lines[1:]:
        line = line.strip()
        lst = line.split()
        if len(lst) < 4: # this should not occur. Report error if occurred.
            return -1
        tpm = lst[3]
        if not tpm == '0' and not 'nan' in tpm:
            non_zero_count += 1
        total_count += 1
    return 1.0 * non_zero_count / total_count

def get_rna_word_count(s):
    ''' If s looks like description of an RNA-seq experiment, return 1 otherwise return 0. '''
    count = 0
    s = s.lower()
    if 'rna-seq' in s or 'rnaseq' in s or 'transcriptome' in s or 'transcript' in s or 'mrna' in s:
        count = 1
    return count

def read_ena_data_info_json(fname):
    d = {}
    with open(fname) as json_data:
        json_dict = json.load(json_data)
        for run_id in json_dict:
           d[run_id] = 1
    return d


def write_log_file(s, fname):
    f = open(fname, 'a')
    curr_time = datetime.now().strftime('%Y-%m-%d %H:%M')
    s = '[' + curr_time + ']: ' + s
    if not '\n' in s:
        s += '\n'
    f.write(s)
    f.close()

### main
if not os.path.exists(RNA_SEQ_INFO_FILE):
    print('make_parameter_rnaseq.py: you must provide %s. See parse_ena_xml.py on how to make it.' % (RNA_SEQ_INFO_FILE))
    sys.exit()
    
rna_data_info_dict = read_ena_data_info_json(RNA_SEQ_INFO_FILE)


print('%%GENE_LIST=../Data/information/gene-list.txt\n%%HOLDON=NO\n') # head

if len(sys.argv) > 1: # only get these samples specified in id-list.txt
    quant_files = get_quant_file_list(sys.argv[1])
else:
    quant_files = glob_files_include_paths(QUANT_PATH, '*_quant.txt')


include_count = 0
total = 0
already_added_dict = {}
for fn in sorted(quant_files):
    total += 1

    nzr = non_zero_ratio(fn)
    if nzr > NON_ZERO_RATIO: # files with too many zeros are ignored
        fn2 = os.path.basename(fn)
        data_name = 'None'
        if 'inhouse' in fn:
            myid = extract_id(fn, 'inhouse')
            myid2 = make_id('0','001', myid) # lab id is 001, phil's lab
            data_name = 'LAB001' + '_LIB' + myid[-5:]
            desc = 'in-house'
            include_me = True
        elif 'public' in fn:
            myid = extract_id(fn2, 'sra')                
            myid2 = make_id('0',myid[0:3], myid[3:])
            data_name = myid
            desc = 'SRA'
            include_me = True if myid in rna_data_info_dict and rna_data_info_dict[myid] >= 0 else False # IMPORTANT
        elif 'other' in fn:
            myid = extract_id(fn, 'other')
            # get lab id
            dirname = os.path.dirname(fn)
            lst = dirname.split('/')
            lab_id = lst[-2]
            myid2 = make_id('0', lab_id, myid) # lab id is 001, phil's lab
            data_name = 'LAB' + lab_id + '_LIB' + myid[-5:]
            desc = 'Other lab'
            include_me = True
        if include_me and not myid2 in already_added_dict:
            print('@%s' % (myid2))
            print('DATA_NAME:%s' % (data_name))
            print('DATA_FORMAT:%s'% ('txt'))
            print('DESCRIPTION:%s' % (desc))
            print('LOCATION:%s' % (fn))
            print('NOTE: non zero ratio is %4.2f' % (nzr))
            print('')
            include_count += 1
            already_added_dict[myid2] = 'yes'
    elif nzr < 0:
        write_log_file('[make_parameter_rnaseq.py] Warning: incomplete line in file %s' % (fn), UPDATE_NETWORK_LOG_FILE)        
    else:
        #print('%s has too many zeroes. ignore.' % (fn))
        pass

print('#Done.  Processed %d files.  Included %d files.' % (total, include_count))