mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-10-28 02:58:18 +08:00
239 lines
7.8 KiB
Python
239 lines
7.8 KiB
Python
#!/usr/bin/env python
|
|
from optparse import OptionParser, SUPPRESS_HELP
|
|
import os, random, quake
|
|
|
|
############################################################
|
|
# cov_model.py
|
|
#
|
|
# Given a file of kmer counts, reports the cutoff to use
|
|
# to separate trusted/untrusted kmers.
|
|
############################################################
|
|
|
|
############################################################
|
|
# main
|
|
############################################################
|
|
def main():
|
|
usage = 'usage: %prog [options] <counts file>'
|
|
parser = OptionParser(usage)
|
|
parser.add_option('--int', dest='count_kmers', action='store_true', default=False, help='Kmers were counted as integers w/o the use of quality values [default: %default]')
|
|
parser.add_option('--ratio', dest='ratio', type='int', default=200, help='Likelihood ratio to set trusted/untrusted cutoff [default: %default]')
|
|
parser.add_option('--no_sample', dest='no_sample', action='store_true', default=False, help='Do not sample kmer coverages into kmers.txt because its already done [default: %default]')
|
|
# help='Model kmer coverage as a function of GC content of kmers [default: %default]'
|
|
parser.add_option('--gc', dest='model_gc', action='store_true', default=False, help=SUPPRESS_HELP)
|
|
(options, args) = parser.parse_args()
|
|
|
|
if len(args) != 1:
|
|
parser.error('Must provide kmers counts file')
|
|
else:
|
|
ctsf = args[0]
|
|
|
|
if options.count_kmers:
|
|
model_cutoff(ctsf, options.ratio)
|
|
print 'Cutoff: %s' % open('cutoff.txt').readline().rstrip()
|
|
|
|
else:
|
|
if options.model_gc:
|
|
model_q_gc_cutoffs(ctsf, 25000, options.ratio)
|
|
else:
|
|
model_q_cutoff(ctsf, 50000, options.ratio, options.no_sample)
|
|
print 'Cutoff: %s' % open('cutoff.txt').readline().rstrip()
|
|
|
|
|
|
############################################################
|
|
# model_cutoff
|
|
#
|
|
# Make a histogram of kmers to give to R to learn the cutoff
|
|
############################################################
|
|
def model_cutoff(ctsf, ratio):
|
|
# make kmer histogram
|
|
cov_max = 0
|
|
for line in open(ctsf):
|
|
cov = int(line.split()[1])
|
|
if cov > cov_max:
|
|
cov_max = cov
|
|
|
|
kmer_hist = [0]*cov_max
|
|
for line in open(ctsf):
|
|
cov = int(line.split()[1])
|
|
kmer_hist[cov-1] += 1
|
|
|
|
cov_out = open('kmers.hist', 'w')
|
|
for cov in range(0,cov_max):
|
|
if kmer_hist[cov]:
|
|
print >> cov_out, '%d\t%d' % (cov+1,kmer_hist[cov])
|
|
cov_out.close()
|
|
|
|
os.system('R --slave --args %d < %s/cov_model.r 2> r.log' % (ratio,quake.quake_dir))
|
|
|
|
|
|
############################################################
|
|
# model_q_cutoff
|
|
#
|
|
# Sample kmers to give to R to learn the cutoff
|
|
# 'div100' is necessary when the number of kmers is too
|
|
# large for random.sample, so we only consider every 100th
|
|
# kmer.
|
|
############################################################
|
|
def model_q_cutoff(ctsf, sample, ratio, no_sample=False):
|
|
if not no_sample:
|
|
# count number of kmer coverages
|
|
num_covs = 0
|
|
for line in open(ctsf):
|
|
num_covs += 1
|
|
|
|
# choose random kmer coverages
|
|
div100 = False
|
|
if sample >= num_covs:
|
|
rand_covs = range(num_covs)
|
|
else:
|
|
if num_covs > 1000000000:
|
|
div100 = True
|
|
rand_covs = random.sample(xrange(num_covs/100), sample)
|
|
else:
|
|
rand_covs = random.sample(xrange(num_covs), sample)
|
|
rand_covs.sort()
|
|
|
|
# print to file
|
|
out = open('kmers.txt', 'w')
|
|
kmer_i = 0
|
|
rand_i = 0
|
|
for line in open(ctsf):
|
|
if div100:
|
|
if kmer_i % 100 == 0 and kmer_i/100 == rand_covs[rand_i]:
|
|
print >> out, line.split()[1]
|
|
rand_i += 1
|
|
if rand_i >= sample:
|
|
break
|
|
else:
|
|
if kmer_i == rand_covs[rand_i]:
|
|
print >> out, line.split()[1]
|
|
rand_i += 1
|
|
if rand_i >= sample:
|
|
break
|
|
kmer_i += 1
|
|
out.close()
|
|
|
|
os.system('R --slave --args %d < %s/cov_model_qmer.r 2> r.log' % (ratio,quake.quake_dir))
|
|
|
|
|
|
############################################################
|
|
# model_q_gc_cutoffs
|
|
#
|
|
# Sample kmers to give to R to learn the cutoff for each
|
|
# GC value
|
|
############################################################
|
|
def model_q_gc_cutoffs(ctsf, sample, ratio):
|
|
# count number of kmer coverages at each at
|
|
k = len(open(ctsf).readline().split()[0])
|
|
num_covs_at = [0]*(k+1)
|
|
for line in open(ctsf):
|
|
kmer = line.split()[0]
|
|
num_covs_at[count_at(kmer)] += 1
|
|
|
|
# for each AT bin
|
|
at_cutoffs = []
|
|
for at in range(1,k):
|
|
# sample covs
|
|
if sample >= num_covs_at[at]:
|
|
rand_covs = range(num_covs_at[at])
|
|
else:
|
|
rand_covs = random.sample(xrange(num_covs_at[at]), sample)
|
|
rand_covs.sort()
|
|
|
|
# print to file
|
|
out = open('kmers.txt', 'w')
|
|
kmer_i = 0
|
|
rand_i = 0
|
|
for line in open(ctsf):
|
|
(kmer,cov) = line.split()
|
|
if count_at(kmer) == at:
|
|
if kmer_i == rand_covs[rand_i]:
|
|
print >> out, cov
|
|
rand_i += 1
|
|
if rand_i >= sample:
|
|
break
|
|
kmer_i += 1
|
|
out.close()
|
|
|
|
os.system('R --slave --args %d < %s/cov_model_qmer.r 2> r%d.log' % (ratio,quake.quake_dir,at))
|
|
|
|
at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
|
|
if at in [1,k-1]: # setting extremes to next closests
|
|
at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
|
|
|
|
os.system('mv kmers.txt kmers.at%d.txt' % at)
|
|
os.system('mv cutoff.txt cutoff.at%d.txt' % at)
|
|
|
|
out = open('cutoffs.gc.txt','w')
|
|
print >> out, '\n'.join(at_cutoffs)
|
|
out.close()
|
|
|
|
|
|
############################################################
|
|
# model_q_gc_cutoffs_bigmem
|
|
#
|
|
# Sample kmers to give to R to learn the cutoff for each
|
|
# GC value
|
|
############################################################
|
|
def model_q_gc_cutoffs_bigmem(ctsf, sample, ratio):
|
|
# input coverages
|
|
k = 0
|
|
for line in open(ctsf):
|
|
(kmer,cov) = line.split()
|
|
if k == 0:
|
|
k = len(kmer)
|
|
at_covs = ['']*(k+1)
|
|
else:
|
|
at = count_at(kmer)
|
|
if at_covs[at]:
|
|
at_covs[at].append(cov)
|
|
else:
|
|
at_covs[at] = [cov]
|
|
|
|
for at in range(1,k):
|
|
print '%d %d' % (at,len(at_covs[at]))
|
|
|
|
# for each AT bin
|
|
at_cutoffs = []
|
|
for at in range(1,k):
|
|
# sample covs
|
|
if sample >= len(at_covs[at]):
|
|
rand_covs = at_covs[at]
|
|
else:
|
|
rand_covs = random.sample(at_covs[at], sample)
|
|
|
|
# print to file
|
|
out = open('kmers.txt', 'w')
|
|
for rc in rand_covs:
|
|
print >> out, rc
|
|
out.close()
|
|
|
|
os.system('R --slave --args %d < %s/cov_model_qmer.r 2> r%d.log' % (ratio,quake.quake_dir,at))
|
|
|
|
at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
|
|
if at in [1,k-1]: # setting extremes to next closests
|
|
at_cutoffs.append( open('cutoff.txt').readline().rstrip() )
|
|
|
|
os.system('mv kmers.txt kmers.at%d.txt' % at)
|
|
os.system('mv cutoff.txt cutoff.at%d.txt' % at)
|
|
|
|
out = open('cutoffs.gc.txt','w')
|
|
print >> out, '\n'.join(at_cutoffs)
|
|
out.close()
|
|
|
|
|
|
############################################################
|
|
# count_at
|
|
#
|
|
# Count A's and T's in the given sequence
|
|
############################################################
|
|
def count_at(seq):
|
|
return len([nt for nt in seq if nt in ['A','T']])
|
|
|
|
|
|
############################################################
|
|
# __main__
|
|
############################################################
|
|
if __name__ == '__main__':
|
|
main()
|