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
|
"""
Test the qtlfiles export of traits files
Run with:
env SQL_URI="mysql://<user>:<password>@<host>:<port>/db_webqtl" python3 qtlfilesexport.py
replacing the variables in the angled brackets with the appropriate values
"""
from gn3.computations.slink import slink
from gn3.db_utils import database_connector
from gn3.computations.qtlreaper import run_reaper
from gn3.computations.heatmap import export_trait_data
from gn3.db.traits import retrieve_trait_data, retrieve_trait_info
from gn3.db.genotypes import build_genotype_file, load_genotype_samples
from gn3.computations.qtlreaper import random_string, generate_traits_file
from gn3.computations.heatmap import (
cluster_traits,
compute_heatmap_order,
retrieve_strains_and_values)
TMPDIR = "tmp/qtltests"
def trait_fullnames():
"""Return sample names for traits"""
return [
"UCLA_BXDBXH_CARTILAGE_V2::ILM103710672",
"UCLA_BXDBXH_CARTILAGE_V2::ILM2260338",
"UCLA_BXDBXH_CARTILAGE_V2::ILM3140576",
"UCLA_BXDBXH_CARTILAGE_V2::ILM5670577",
"UCLA_BXDBXH_CARTILAGE_V2::ILM2070121",
"UCLA_BXDBXH_CARTILAGE_V2::ILM103990541",
"UCLA_BXDBXH_CARTILAGE_V2::ILM1190722",
"UCLA_BXDBXH_CARTILAGE_V2::ILM6590722",
"UCLA_BXDBXH_CARTILAGE_V2::ILM4200064",
"UCLA_BXDBXH_CARTILAGE_V2::ILM3140463"]
def main():
"""entrypoint function"""
conn = database_connector()[0]
threshold = 0
traits = [
retrieve_trait_info(threshold, fullname, conn)
for fullname in trait_fullnames()]
traits_data_list = [retrieve_trait_data(t, conn) for t in traits]
genotype_filename = build_genotype_file(traits[0]["riset"])
strains = load_genotype_samples(genotype_filename)
exported_traits_data_list = [
export_trait_data(td, strains) for td in traits_data_list]
slinked = slink(cluster_traits(exported_traits_data_list))
orders = compute_heatmap_order(slinked)
strains_and_values = retrieve_strains_and_values(
orders, strains, exported_traits_data_list)
strains_values = strains_and_values[0][1]
trait_values = [t[2] for t in strains_and_values]
traits_filename = "{}/traits_test_file_{}.txt".format(
TMPDIR, random_string(10))
generate_traits_file(strains_values, trait_values, traits_filename)
print("Generated file: {}".format(traits_filename))
main_output, permutations_output = run_reaper(
genotype_filename, traits_filename, separate_nperm_output=True)
print("Main output: {}, Permutation output: {}".format(
main_output, permutations_output))
if __name__ == "__main__":
main()
|
