working commit of the fmx-assign-to-gt functionality

single_cell_RNAseq/bin/examine_gtcheck.R

@@ -0,0 +1,62 @@
+library(tidyverse)
+library(gplots)
+
+args = commandArgs(trailingOnly=T)
+POOL = args[1]
+gtcheck_out = args[2]
+
+tab_in = read_tsv(gtcheck_out, comment="#", col_names=F) 
+if (ncol(tab_in) == 6 & "CN" %in% tab_in$X1){ ### for 1.3.1
+  my_tab =  tab_in %>%
+    filter(X1=="CN") %>%
+    filter(!str_detect(X5, "CLUST"), str_detect(X6, "CLUST")) %>%
+    select(-X1, -X3, -X4) %>%
+    mutate(X5=str_replace_all(X5, "_[0-9]+", "")) %>%
+    dplyr::rename(individual=X5, cluster=X6) %>%
+    mutate(err=as.numeric(X2)) %>%
+    select(-X2)
+} else if (ncol(tab_in)==6 & "DC" %in% tab_in$X1){ # 1.18
+  my_tab = tab_in %>% 
+    select(-X1) %>%
+    filter(X6!=0, str_detect(X3, "CLUST"), !str_detect(X2, "CLUST")) %>%
+    rename(individual=X2, cluster=X3, err=X4) %>%
+    select(-X5, -X6)
+} else { # 1.10
+  my_tab = tab_in %>%
+    filter(!str_detect(X4, "CLUST"), str_detect(X5, "CLUST"), X1=="ERR") %>%
+    select(-X1, -X3) %>%
+    mutate(X4=str_replace_all(X4, "_[0-9]+", "")) %>%
+    dplyr::rename(individual=X4, cluster=X5) %>%
+    mutate(err=as.numeric(X2)) %>%
+    select(-X2)
+}
+
+
+my_mat = my_tab %>% pivot_wider(names_from="cluster", values_from="err", values_fill=NA) %>%
+  column_to_rownames("individual") %>%
+  as.matrix() %>%
+  t() 
+if(ncol(my_mat)==1 | nrow(my_mat) == 1 | any(is.na(my_mat))){
+  print("Warning - unable to plot because too few values in heatmap.")
+
+} else {
+  pdf(sprintf("%s_gtcheck_out.pdf", POOL), height=5, width=5)
+  heatmap.2(my_mat, trace="none", scale="none", cexRow = 0.5, cexCol = 0.5)
+  dev.off()
+}
+
+
+cutoff = mean(my_tab$err) - sd(my_tab$err)
+
+
+
+filt_tab = my_tab %>% 
+  filter(err < cutoff) 
+filt_tab %>%
+  write_tsv(sprintf("%s_gtcheck_map_to_ind.tsv", POOL))
+
+if (length(unique(filt_tab$individual)) != nrow(filt_tab) | 
+    length(unique(filt_tab$cluster)) != nrow(filt_tab) | 
+    nrow(filt_tab) != nrow(my_mat) ){
+  print("Error, the number of assigned samples and clusters does not match the total")
+}

single_cell_RNAseq/bin/run_gtcheck.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+
+POOL=$1
+REF_VCF_FILE=$2
+FMX_VCF_FILE=$3
+
+cp ${REF_VCF_FILE} ground_truth.vcf.gz
+bcftools view -O b -o ground_truth.bcf ground_truth.vcf.gz
+bcftools view -O b -o ${POOL}.bcf ${FMX_VCF_FILE}
+bcftools index ground_truth.bcf
+bcftools index ${POOL}.bcf
+
+bcftools filter --exclude "GQ<20" -O b -o ${POOL}_filtered.bcf ${POOL}.bcf
+bcftools index ${POOL}_filtered.bcf
+
+bcftools isec -O b -p ${POOL}_with_ground_truth -n =2 ${POOL}_filtered.bcf ground_truth.bcf
+
+bcftools merge --merge snps -o pool_with_ground_truth.bcf -O b ${POOL}_with_ground_truth/0000.bcf \
+        ${POOL}_with_ground_truth/0001.bcf
+bcftools index pool_with_ground_truth.bcf
+
+bcftools gtcheck pool_with_ground_truth.bcf | grep -v "^INFO" > ${POOL}_gtcheck.out

single_cell_RNAseq/config/container.config

@@ -4,6 +4,7 @@ params {
         cellranger = "${params.container.singularity_dir}/cellranger/6.0.2/cellranger.sif"
         popscle = "${params.container.singularity_dir}/popscle/da70fc78da385ef049e0e890342acfd62842cae0/popscle.sif"
         rsinglecell = "${params.container.singularity_dir}/RSingleCell/v4/RSingleCell.sif"
+        rplus_bcftools = "${params.container.singularity_dir}/rplus_bcftools/rsinglecell_bcftools.sif"
         python = "${params.container.singularity_dir}/cytof/v3/cytof.sif"
     }
 }