add many more tests and R ports

.github/workflows/ci.yml

@@ -26,9 +26,14 @@ jobs:
 
       - name: Set up pixi environment
         uses: prefix-dev/setup-pixi@v0.10.0
+        with:
+          run-install: false
 
       - name: Script tests
-        run: pixi run test_scripts
+        run: |
+          pixi config set --local run-post-link-scripts insecure
+          pixi install
+          pixi run test_scripts
 
       - name: Notebook tests
         run: pixi run test_notebooks

.gitignore

@@ -25,6 +25,8 @@ results/
 
 # Data files (too large for git)
 *.parquet
+# committed test fixtures are the exception
+!tests/fixtures/**/*.parquet
 data/
 !renovated_code/snakemake/modular_sos/tests/data/
 !renovated_code/snakemake/modular_sos/tests/data/*.tar.gz

code/SoS/data_preprocessing/covariate/covariate_hidden_factor.ipynb

@@ -495,23 +495,14 @@
    "outputs": [],
    "source": [
     "[BiCV_2]\n",
-    "# For cluster jobs, number commands to run per job\n",
-    "import time\n",
     "output: f'{cwd:a}/{_input:bn}.fake.vcf.gz'\n",
-    "#task: trunk_workers = 1,trunk_size = job_size , walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
-    "R: container=container, expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', entrypoint = entrypoint\n",
-    "    library(\"dplyr\")\n",
-    "    library(\"readr\")\n",
-    "    ## Add fake header\n",
-    "    cat(paste(\"##fileformat=VCFv4.2\\n##fileDate=$[time.strftime(\"%Y%m%d\",time.localtime())]\\n##source=FAKE\\n\"), file=$[_output:nr], append=FALSE)\n",
-    "    ## Add colnames based on bed\n",
-    "    pheno = read_delim(\"$[_input]\", delim = \"\\t\",n_max = 1)\n",
-    "    colnames(pheno)[1:3] = c(\"#CHROM\",\"POS\",\"ID\") \n",
-    "    pheno = cbind(pheno[1:3]%>%mutate(REF = \"A\", ALT = \"C\", QUAL = \".\",FILTER = \".\", INFO = \"PR\", FORMAT = \"GT\"),pheno[,5:ncol(pheno)])\n",
-    "    pheno%>%write_delim($[_output:nr],delim = \"\\t\",col_names = T, append = T)\n",
-    "bash: container=container, expand= \"$[ ]\", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', entrypoint = entrypoint\n",
-    "    bgzip -f $[_output:n]\n",
-    "    tabix -p vcf $[_output]"
+    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
+    "bash: expand= \"${ }\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container, entrypoint = entrypoint\n",
+    "    Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \\\n",
+    "        --step BiCV_2 \\\n",
+    "        --cwd \"${cwd}\" \\\n",
+    "        --residFile \"${_input}\" \\\n",
+    "        --numThreads ${numThreads}"
    ]
   },
   {
@@ -530,29 +521,20 @@
     "parameter: N = 0\n",
     "# The number of iteration\n",
     "parameter: iteration = 10\n",
-    "\n",
-    "import pandas as pd\n",
-    "data = pd.read_csv(phenoFile, sep=\"\t\", index_col = 3).drop([\"#chr\",\"start\",\"end\"],axis = 1)\n",
-    "if N == 0:\n",
-    "    if len(data.columns) < 150:\n",
-    "        N = 15\n",
-    "    elif len(data.columns) < 250:\n",
-    "        N = 30\n",
-    "    elif len(data.columns) < 350:\n",
-    "        N = 45\n",
-    "    else:\n",
-    "        N = 60\n",
     "input: output_from(1), output_from(2)\n",
     "output: f'{cwd:a}/{_input[0]:bnn}.BiCV.gz'\n",
     "task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
-    "bash: container=container, expand= \"$[ ]\", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', entrypoint = entrypoint\n",
-    "    apex factor \\\n",
-    "        --out $[_output[0]:nn] \\\n",
-    "        --iter $[iteration] \\\n",
-    "        --factors $[N] \\\n",
-    "        --bed $[_input[0]] \\\n",
-    "        --vcf $[_input[1]] \\\n",
-    "        --threads $[numThreads]  $[ f'--cov {covFile}' if covFile.is_file() else '']"
+    "bash: expand= \"${ }\", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint\n",
+    "    Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \\\n",
+    "        --step BiCV_3 \\\n",
+    "        --cwd \"${cwd}\" \\\n",
+    "        --residFile \"${_input[0]}\" \\\n",
+    "        --vcfFile \"${_input[1]}\" \\\n",
+    "        --covFile \"${covFile}\" \\\n",
+    "        --N ${N} \\\n",
+    "        --iteration ${iteration} \\\n",
+    "        --output \"${_output[0]}\" \\\n",
+    "        --numThreads ${numThreads}"
    ]
   },
   {

code/SoS/data_preprocessing/genotype/GRM.ipynb

@@ -98,9 +98,7 @@
    },
    "outputs": [],
    "source": [
-    "sos run pipeline/GRM.ipynb grm \\\n",
-    "    --cwd output/grm_uf \\\n",
-    "    --genoFile input/genotype/protocol_example.genotype.list"
+    "sos run pipeline/GRM.ipynb grm \\\n    --cwd output/grm_uf \\\n    --genoFile input/genotype/protocol_example.genotype.list"
    ]
   },
   {
@@ -132,72 +130,7 @@
    },
    "outputs": [],
    "source": [
-    "[global]\n",
-    "# Work directory & output directory\n",
-    "parameter: cwd = path(\"output\")\n",
-    "# The filename name for output data\n",
-    "parameter: container = ''\n",
-    "# For cluster jobs, number commands to run per job\n",
-    "parameter: job_size = 1\n",
-    "# Wall clock time expected\n",
-    "parameter: walltime = \"5h\"\n",
-    "# Memory expected\n",
-    "parameter: mem = \"16G\"\n",
-    "# Number of threads\n",
-    "parameter: numThreads = 20\n",
-    "parameter: genoFile = paths\n",
-    "cwd = f\"{cwd:a}\"\n",
-    "\n",
-    "## Code below are copied from genotype_formatting.ipynb. Change should be made to codes in that document and sync with this chunk.\n",
-    "from sos.utils import expand_size\n",
-    "import os\n",
-    "\n",
-    "def get_genotype_file(geno_file_paths):\n",
-    "    #\n",
-    "    def valid_geno_file(x):\n",
-    "        suffixes = path(x).suffixes\n",
-    "        if suffixes[-1] == '.bed':\n",
-    "            return True\n",
-    "        if len(suffixes)>1 and ''.join(suffixes[-2:]) == \".vcf.gz\":\n",
-    "            return True\n",
-    "        return False\n",
-    "    #\n",
-    "    def complete_geno_path(x, geno_file):\n",
-    "        if not valid_geno_file(x):\n",
-    "            raise ValueError(f\"Genotype file {x} should be VCF (end with .vcf.gz) or PLINK bed file (end with .bed)\")\n",
-    "        if not os.path.isfile(x):\n",
-    "            # relative path\n",
-    "            if not os.path.isfile(f'{geno_file:ad}/' + x):\n",
-    "                raise ValueError(f\"Cannot find genotype file {x}\")\n",
-    "            else:\n",
-    "                x = f'{geno_file:ad}/' + x\n",
-    "        return x\n",
-    "    # \n",
-    "    def format_chrom(chrom):\n",
-    "        if chrom.startswith('chr'):\n",
-    "            chrom = chrom[3:]\n",
-    "        return chrom\n",
-    "    # Inputs are either VCF or bed, or a vector of them \n",
-    "    if len(geno_file_paths) > 1:\n",
-    "        if all([valid_geno_file(x) for x in geno_file_paths]):\n",
-    "            return paths(geno_file_paths)\n",
-    "        else: \n",
-    "            raise ValueError(f\"Invalid input {geno_file_paths}\")\n",
-    "    # Input is one genotype file or text list of genotype files\n",
-    "    geno_file = geno_file_paths[0]\n",
-    "    if valid_geno_file(geno_file):\n",
-    "        return path(geno_file)\n",
-    "    else: \n",
-    "        units = [x.strip().split() for x in open(geno_file).readlines() if x.strip() and not x.strip().startswith('#')]\n",
-    "        if all([len(x) == 1 for x in units]):\n",
-    "            return paths([complete_geno_path(x[0], geno_file) for x in units])\n",
-    "        elif all([len(x) == 2 for x in units]):\n",
-    "            genos = dict([(format_chrom(x[0]), path(complete_geno_path(x[1], geno_file))) for x in units])\n",
-    "        else:\n",
-    "            raise ValueError(f\"{geno_file} should contain one column of file names, or two columns of chrom number and corresponding file name\")\n",
-    "        return genos\n",
-    "                        \n",
-    "genotypes = get_genotype_file(genoFile)"
+    "[global]\n# Work directory & output directory\nparameter: cwd = path(\"output\")\n# The filename name for output data\nparameter: modular_script_dir = path('code/script')\n# For cluster jobs, number commands to run per job\nparameter: job_size = 1\n# Wall clock time expected\nparameter: walltime = \"5h\"\n# Memory expected\nparameter: mem = \"16G\"\n# Number of threads\nparameter: numThreads = 20\nparameter: genoFile = paths\ncwd = f\"{cwd:a}\"\n\n## Code below are copied from genotype_formatting.ipynb. Change should be made to codes in that document and sync with this chunk.\nfrom sos.utils import expand_size\nimport os\n\ndef get_genotype_file(geno_file_paths):\n    #\n    def valid_geno_file(x):\n        suffixes = path(x).suffixes\n        if suffixes[-1] == '.bed':\n            return True\n        if len(suffixes)>1 and ''.join(suffixes[-2:]) == \".vcf.gz\":\n            return True\n        return False\n    #\n    def complete_geno_path(x, geno_file):\n        if not valid_geno_file(x):\n            raise ValueError(f\"Genotype file {x} should be VCF (end with .vcf.gz) or PLINK bed file (end with .bed)\")\n        if not os.path.isfile(x):\n            # relative path\n            if not os.path.isfile(f'{geno_file:ad}/' + x):\n                raise ValueError(f\"Cannot find genotype file {x}\")\n            else:\n                x = f'{geno_file:ad}/' + x\n        return x\n    # \n    def format_chrom(chrom):\n        if chrom.startswith('chr'):\n            chrom = chrom[3:]\n        return chrom\n    # Inputs are either VCF or bed, or a vector of them \n    if len(geno_file_paths) > 1:\n        if all([valid_geno_file(x) for x in geno_file_paths]):\n            return paths(geno_file_paths)\n        else: \n            raise ValueError(f\"Invalid input {geno_file_paths}\")\n    # Input is one genotype file or text list of genotype files\n    geno_file = geno_file_paths[0]\n    if valid_geno_file(geno_file):\n        return path(geno_file)\n    else: \n        units = [x.strip().split() for x in open(geno_file).readlines() if x.strip() and not x.strip().startswith('#')]\n        if all([len(x) == 1 for x in units]):\n            return paths([complete_geno_path(x[0], geno_file) for x in units])\n        elif all([len(x) == 2 for x in units]):\n            genos = dict([(format_chrom(x[0]), path(complete_geno_path(x[1], geno_file))) for x in units])\n        else:\n            raise ValueError(f\"{geno_file} should contain one column of file names, or two columns of chrom number and corresponding file name\")\n        return genos\n                        \ngenotypes = get_genotype_file(genoFile)"
    ]
   },
   {
@@ -218,13 +151,7 @@
    },
    "outputs": [],
    "source": [
-    "# Generate LOCO file list\n",
-    "[grm_1]\n",
-    "chrom = list(genotypes.keys())\n",
-    "input: for_each = 'chrom'\n",
-    "output: f'{cwd}/{genotypes[_chrom]:bn}.loco.txt'\n",
-    "with open(_output, 'w') as f:\n",
-    "    f.write('\\n'.join([str(genotypes[x].with_suffix('')) for x in genotypes if x != _chrom]))"
+    "# Generate LOCO file list\n[grm_1]\nchrom = list(genotypes.keys())\ninput: for_each = 'chrom'\noutput: f'{cwd}/{genotypes[_chrom]:bn}.loco.txt'\nwith open(_output, 'w') as f:\n    f.write('\\n'.join([str(genotypes[x].with_suffix('')) for x in genotypes if x != _chrom]))"
    ]
   },
   {
@@ -245,15 +172,7 @@
    },
    "outputs": [],
    "source": [
-    "# Compute LOCO-GRM\n",
-    "[grm_2]\n",
-    "output: f'{_input:nn}.grm.gz'\n",
-    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
-    "bash: expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container\n",
-    "   gcta64 \\\n",
-    "   --mbfile $[_input] \\\n",
-    "   --make-grm-gz \\\n",
-    "   --out $[_output:nn]"
+    "# Compute LOCO-GRM\n[grm_2]\noutput: f'{_input:nn}.grm.gz'\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\nbash: expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'\n   gcta \\\n   --mbfile $[_input] \\\n   --make-grm-gz \\\n   --out $[_output:nn]\n"
    ]
   },
   {
@@ -274,33 +193,7 @@
    },
    "outputs": [],
    "source": [
-    "# Format output for APEX\n",
-    "[grm_3]\n",
-    "output: f'{_input:nn}.apex.grm'\n",
-    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
-    "R: expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container\n",
-    "\n",
-    "    # Load necessary libraries\n",
-    "    library(dplyr)\n",
-    "    library(readr)\n",
-    "    library(purrr)\n",
-    "\n",
-    "    # Read data\n",
-    "    grm <- read_delim(gzfile(\"$[_input]\"), '\\t', col_names = FALSE) %>%\n",
-    "      select(1, 2, 4) %>%\n",
-    "      set_names(c(\"#id1\", \"id2\", \"kinship\"))\n",
-    "\n",
-    "    id <- read_delim(\"$[_input:n].id\", '\\t', col_names = FALSE)[, 2] %>%\n",
-    "      unlist()\n",
-    "\n",
-    "    # Modify grm data\n",
-    "    grm_updated <- grm %>%\n",
-    "      mutate(`#id1` = map_chr(`#id1`, ~id[.x]),\n",
-    "             `id2` = map_chr(`id2`, ~id[.x]))\n",
-    "\n",
-    "    # Write to output\n",
-    "    grm_updated %>%\n",
-    "      write_delim(\"$[_output]\", '\\t')"
+    "# Format output for APEX\n[grm_3]\noutput: f'{_input:nn}.apex.grm'\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\nbash: expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'\n    Rscript $[modular_script_dir:a]/data_preprocessing/genotype/GRM.R --step format_apex \\\n        --grm $[_input] \\\n        --id $[_input:n].id \\\n        --output $[_output]\n"
    ]
   },
   {
@@ -328,30 +221,7 @@
    },
    "outputs": [],
    "source": [
-    "# Generate APEX input list\n",
-    "[grm_4]\n",
-    "input: group_by = \"all\"\n",
-    "output: f'{cwd}/{path(list(genotypes.values())[1]):bnn}.loco_grm_list.txt'\n",
-    "R: expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container\n",
-    "\n",
-    "    # Load necessary libraries\n",
-    "    library(dplyr)\n",
-    "    library(readr)\n",
-    "    library(purrr)\n",
-    "\n",
-    "    # Define chromosomal names\n",
-    "    chromosomes <- paste0(\"chr\", c($[paths(genotypes.keys()):r,]), collapse = \"\")\n",
-    "\n",
-    "    # Define input files and directory\n",
-    "    input_files <- c($[_input:r,])\n",
-    "\n",
-    "    # Create and sort the tibble\n",
-    "    grm_tibble <- tibble(`#chr` = chromosomes, dir = input_files) %>%\n",
-    "      arrange(`#chr`)\n",
-    "\n",
-    "    # Write to output\n",
-    "    grm_tibble %>%\n",
-    "      write_delim(\"$[_output]\", \"\\t\")"
+    "# Generate APEX input list\n[grm_4]\ninput: group_by = \"all\"\noutput: f'{cwd}/{path(list(genotypes.values())[1]):bnn}.loco_grm_list.txt'\nbash: expand= \"$[ ]\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'\n    Rscript $[modular_script_dir:a]/data_preprocessing/genotype/GRM.R --step apex_list \\\n        --input $[_input:r,] \\\n        --chroms $[\" \".join(genotypes.keys())] \\\n        --output $[_output]\n"
    ]
   }
  ],
@@ -391,4 +261,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}