comparison prepare_population_structure.py @ 27:8997f2ca8c7a

Update to Miller Lab devshed revision bae0d3306d3b

26:91e835060ad2

#!/usr/bin/env python

import errno
import os
import shutil
import subprocess
import sys
from Population import Population
import gd_composite

################################################################################

def do_import(filename, files_path, min_reads, min_qual, min_spacing, tags, using_info, population_list):
    info_page = gd_composite.InfoPage()
    info_page.set_title('Prepare to look for population structure Galaxy Composite Dataset')

    display_file = gd_composite.DisplayFile()
    display_value = gd_composite.DisplayValue()

    info_page.add_misc(misc_populations)

    with open(filename, 'w') as ofh:
        print >> ofh, info_page.render()

def mkdir_p(path):
    try:
        os.makedirs(path)
    except OSError, e:
        if e.errno <> errno.EEXIST:
            raise

def die(message, exit=True):
    print >> sys.stderr, message
    if exit:
        sys.exit(1)

################################################################################

if len(sys.argv) < 10:
    die("Usage")

# parse command line
input_snp_filename, input_type, min_reads, min_qual, min_spacing, output_filename, output_files_path = sys.argv[1:8]
args = sys.argv[8:]

individual_metadata = []
population_files = []
population_names = []
all_individuals = False

for arg in args:
    if arg == 'all_individuals':
        all_individuals = True
    elif len(arg) > 11:
        tag = arg[:11]
        value = arg[11:]
        if tag == 'individual:':
            individual_metadata.append(value)
        elif tag == 'population:':
            filename, name = value.split(':', 1)
            population_files.append(filename)
            population_names.append(name)

p_total = Population()
p_total.from_tag_list(individual_metadata)

individual_population = {}

population_list = []

if all_individuals:
    p1 = p_total
    p1.name = 'All Individuals'
    population_list.append(p1)
else:
    p1 = Population()
    for idx in range(len(population_files)):
        population_file = population_files[idx]
        population_name = population_names[idx]
        this_pop = Population(population_name)
        this_pop.from_population_file(population_file)
        population_list.append(this_pop)
        p1.from_population_file(population_file)
        tags = p1.tag_list()
        for tag in tags:
            if tag not in individual_population:
                individual_population[tag] = population_name

if not p_total.is_superset(p1):
    print >> sys.stderr, 'There is an individual in the population that is not in the SNP table'
    sys.exit(1)

# run tool
prog = 'admix_prep'

args = []
args.append(prog)
args.append(input_snp_filename)
args.append(min_reads)
args.append(min_qual)
args.append(min_spacing)

tags = p1.tag_list()
for tag in tags:
    if input_type == 'gd_genotype':
        column, name = tag.split(':')
        tag = '{0}:{1}'.format(int(column) - 2, name)
    args.append(tag)

#print "args:", ' '.join(args)
p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=subprocess.PIPE, stderr=sys.stderr)
(stdoutdata, stderrdata) = p.communicate()
rc = p.returncode

if rc != 0:
    die('admix_prep failed: rc={0}'.format(rc))

using_info = stdoutdata.rstrip('\r\n')
mkdir_p(output_files_path)
output_ped_filename = os.path.join(output_files_path, 'admix.ped')
output_map_filename = os.path.join(output_files_path, 'admix.map')
shutil.copy2('admix.ped', output_ped_filename)
shutil.copy2('admix.map', output_map_filename)
do_import(output_filename, output_files_path, min_reads, min_qual, min_spacing, tags, using_info, population_list)

os.unlink('admix.ped')
os.unlink('admix.map')

sys.exit(0)

27:8997f2ca8c7a

#!/usr/bin/env python

import gd_util
import os
import shutil
import sys
from Population import Population
import gd_composite

################################################################################

def do_import(filename, files_path, min_reads, min_qual, min_spacing, using_info, population_list):
    info_page = gd_composite.InfoPage()
    info_page.set_title('Prepare to look for population structure Galaxy Composite Dataset')

    display_file = gd_composite.DisplayFile()
    display_value = gd_composite.DisplayValue()

    info_page.add_misc(misc_populations)

    with open(filename, 'w') as ofh:
        print >> ofh, info_page.render()

################################################################################

if len(sys.argv) < 10:
    gd_util.die('Usage')

# parse command line
input_snp_filename, input_type, min_reads, min_qual, min_spacing, output_filename, output_files_path, ind_arg = sys.argv[1:9]
args = sys.argv[9:]

population_files = []
all_individuals = False

for arg in args:
    if arg == 'all_individuals':
        all_individuals = True
    elif len(arg) > 11 and arg[:11] == 'population:':
        file, name = arg[11:].split(':', 1)
        population_files.append((file, name))

p_total = Population()
p_total.from_wrapped_dict(ind_arg)

individual_population = {}
population_list = []

if all_individuals:
    p1 = p_total
    p1.name = 'All Individuals'
    population_list.append(p1)
else:
    p1 = Population()
    for file, name in population_files:
        this_pop = Population(name)
        this_pop.from_population_file(file)
        population_list.append(this_pop)

        for tag in this_pop.tag_list():
            if tag not in individual_population:
                individual_population[tag] = name

        # add individuals from this file to p1
        p1.from_population_file(file)


if not p_total.is_superset(p1):
    gd_util.die('There is an individual in the population that is not in the SNP table')

################################################################################

prog = 'admix_prep'

args = [ prog ]
args.append(input_snp_filename)
args.append(min_reads)
args.append(min_qual)
args.append(min_spacing)

for tag in p1.tag_list():
    if input_type == 'gd_genotype':
        column, name = tag.split(':', 1)
        tag = '{0}:{1}'.format(int(column) - 2, name)
    args.append(tag)

stdoutdata, stderrdata = gd_util.run_program(prog, args)
using_info = stdoutdata.rstrip('\r\n')

################################################################################

gd_util.mkdir_p(output_files_path)

output_ped_filename = os.path.join(output_files_path, 'admix.ped')
output_map_filename = os.path.join(output_files_path, 'admix.map')
shutil.copy2('admix.ped', output_ped_filename)
shutil.copy2('admix.map', output_map_filename)

do_import(output_filename, output_files_path, min_reads, min_qual, min_spacing, using_info, population_list)
sys.exit(0)