duplicate vcf pos deleted [run ci]

Author: amatur

.github/workflows/linux.yml

@@ -26,7 +26,7 @@ jobs:
           
     - name: Build selscan
       working-directory: src/
-      run: make -f ../unused/static_makefiles/Makefile_linux -j
+      run: make -f Makefile_linux -j
 
     - name: Test selscan installation
       working-directory: src/

README.md

@@ -238,7 +238,15 @@ in the construction of your haplotypes please use the --keep-low-freq flag.
 
 ## 📝 Change Log
 ```
-
+26SEP2025 - selscan v2.1.1 - Bug fixes for the fast and memory-efficient version introduced in v2.1:
+
+    - Refined cutoff handling for edge cases, improving correlation with v2.0 outputs.
+    - Corrected reporting of sites with low minor allele counts (avoids zero-area and infinite nSL/iHS).
+    - Fixed nSL distance cutoff in multi-parameter mode.
+    - XP-EHH now uses genetic distance from map files instead of defaulting to physical distance.
+    - Binaries updated for compatibility with older GLIBC.
+    - Improvements in manual: added a section describing norm usage.
+	
 09APR2025 - selscan v2.1.0 - Introducing fast and memory-efficient versions for all statistics. See Rahman, et al. (2025) Biorxiv for details. 
 
 	There is a new batch option for efficient processing of multiple statistics or parameters at once. See the manual for full details.

Selscan_Manual_v2_1.pdf

@@ -0,0 +1,8 @@
+1	100	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1000	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1000	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1000	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1000	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1000	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1000	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1
+1	1100	rs10	G	A	.	PASS	.	GT	0|1	0|0	0|1	1|1	0|0	1|1	0|1	0|0	0|1	1|1

Selscan_Manual_v2_1_1.pdf

@@ -0,0 +1,154 @@
+plotHaplotypes <- function(filename)
+{
+  #filename <- "outfile.rs7215219.ihs.out";
+  
+  # filename.der <- paste(filename,".der.colormap",sep="")
+  # filename.anc <- paste(filename,".anc.colormap",sep="")
+  # data<-as.matrix(read.table(filename.der))
+  rawData<-read.table(filename)
+
+
+  print("Plotting EHH decay...")
+  setEPS()
+  postscript(file=paste(filename,".ehh.eps",sep=""))
+  plot(rawData$V1/1000000,rawData$V2,pch=" ",xlab="Distance from locus (Mb)",ylab="EHH",ylim=c(0,1))
+  lines(rawData$V1/1000000,rawData$V3,col="red")
+  lines(rawData$V1/1000000,rawData$V4,col="blue")
+  legend("topright",legend=c("Derived","Ancestral"),col=c("red","blue"),lty=1)
+  dev.off()
+  
+  numHaps <- dim(data)[1]
+  numSites <- dim(data)[2]
+  numCols<-range(data)[2]-range(data)[1]
+  
+  pos <- rawData$V1
+  hap <- seq(1,numHaps)
+
+  numSortCols <- min(numCols,5)
+  
+  ordering <- as.data.frame(matrix(rep(0,numSortCols*numHaps),nrow=numHaps,ncol=numSortCols))
+
+  print("Sorting derived colors...")
+  
+  for (h in 1:numSortCols)
+    {
+      for (i in 1:numHaps)
+        {
+          for (j in 1:numSites)
+            {
+              if ( data[i,j] == h-1 )
+                {
+                  ordering[i,h] <- ordering[i,h]+1
+                }
+            }
+        }
+    }
+  
+  order.string <- "sorted.order<-order("
+  for (h in 1:numSortCols)
+    {
+      order.string <- paste(order.string,"ordering$V",h,",",sep="")
+    }
+  substr(order.string,nchar(order.string),nchar(order.string)) <- ")"
+  eval(parse(text = order.string))
+  
+  blank<-rep("",length(pos))
+  
+  
+  data2<-as.matrix(read.table(filename.anc))
+  numHaps2 <- dim(data2)[1]
+  numSites2 <- dim(data2)[2]
+  numCols2 <- range(data2)[2]-range(data2)[1]
+  
+  pos2 <- rawData$V1
+  hap2 <- seq(1,numHaps2)
+
+  numSortCols2 <- min(numCols2,5)
+  
+  ordering2 <- as.data.frame(matrix(rep(0,numSortCols2*numHaps2),nrow=numHaps2,ncol=numSortCols2))
+
+  print("Sorting ancestral colors...")
+  
+  for (h in 1:numSortCols2)
+    {
+      for (i in 1:numHaps2)
+        {
+          for (j in 1:numSites2)
+            {
+              if ( data2[i,j] == h-1 )
+                {
+                  ordering2[i,h] <- ordering2[i,h]+1
+                }
+            }
+        }
+    }
+  
+  order.string2 <- "sorted.order2<-order("
+  for (h in 1:numSortCols2)
+    {
+      order.string2 <- paste(order.string2,"ordering2$V",h,",",sep="")
+    }
+  substr(order.string2,nchar(order.string2),nchar(order.string2)) <- ")"
+  eval(parse(text = order.string2))
+  
+  blank2<-rep("",length(pos2))
+  space <- rep(-1,numSites)
+
+  padding <- as.integer((numHaps+numHaps2) * 0.05)
+  if(padding < 2)
+    {
+      padding <- 2
+    } else if(padding %% 2 == 1)
+    {
+      padding <- padding-1
+    }
+  
+  combined.data <- rbind(data[sorted.order,],space);
+
+  for (i in 1:padding)
+    {
+      combined.data <- rbind(combined.data,space);
+    }
+  
+  combined.data <- rbind(combined.data,data2[sorted.order2,])
+  combined.numHaps <- dim(combined.data)[1]
+  combined.hap <- seq(1,combined.numHaps)
+  combined.numCols <- max(numCols,numCols2)
+  
+  col.der <- c(rgb(227/255, 26/255, 28/255),rgb(51/255, 160/255, 44/255),rgb(31/255, 120/255, 180/255),rgb(255/255, 127/255, 0/255),rgb(106/255, 61/255, 154/255),rgb(251/255, 154/255, 153/255),rgb(178/255, 223/255, 138/255),rgb(166/255, 206/255, 227/255),rgb(253/255, 191/255, 111/255),rgb(202/255, 178/255, 214/255))
+  colors.der <- rep("",combined.numCols)
+  
+  col.anc <- c(rgb(31/255, 120/255, 180/255),rgb(255/255, 127/255, 0/255),rgb(106/255, 61/255, 154/255),rgb(227/255, 26/255, 28/255),rgb(51/255, 160/255, 44/255),rgb(166/255, 206/255, 227/255),rgb(253/255, 191/255, 111/255),rgb(202/255, 178/255, 214/255),rgb(251/255, 154/255, 153/255),rgb(178/255, 223/255, 138/255))
+  colors.anc <- rep("",combined.numCols)
+  
+  index <- 1;
+  for(i in 0:(combined.numCols-1))
+    {
+      colors.der[i+1] <- col.der[(i %% 10) + 1]
+      colors.anc[i+1] <- col.anc[(i %% 10) + 1]
+    }
+
+  pos <- pos/1000000
+  
+  print("Plotting haplotype colors...")
+  setEPS()
+  postscript(file=paste(filename,".hapcolor.eps",sep=""))
+  
+  image(pos,combined.hap,t(combined.data),zlim=c(-0.1,-0.01),yaxt="n",xaxt="n",ylab="",xlab="Distance from locus (Mb)",ylim=c(1-(padding/2),combined.numHaps+(padding/2)))
+  axis(3,at=pos,lab=blank,tck=(1/(combined.numHaps) * padding/4))
+  axis(3)
+  axis(1,at=pos,lab=blank,tck=(1/(combined.numHaps) * padding/4))
+  axis(1)
+  abline(h=numHaps+(padding/2)+1)
+  
+  image(pos,combined.hap[1:numHaps],t(combined.data[1:numHaps,]),zlim=c(0,max(combined.data)),col=colors.der,add=TRUE)#,yaxt="n",xaxt="n",ylab="",ylim=c(1-(padding/10),combined.numHaps+(padding/10)),xlab="")
+  mtext("Derived",2,at=(numHaps/2 + padding/2))  
+  image(pos,combined.hap[(numHaps+1):(combined.numHaps)],t(combined.data[(numHaps+1):(combined.numHaps),]),zlim=c(0,max(combined.data)),col=colors.anc,add=TRUE)#,yaxt="n",xaxt="n",ylab="",ylim=c(-1,combined.numHaps+2),xlab="")
+  mtext("Ancestral",2,at= (numHaps + (padding+1) + numHaps2/2 ) )
+
+  dev.off()
+
+  return
+}
+
+plotHaplotypes("/Users/amatur/code/selscan/src/out/outfile.ihs.out")

example/multiallele2.vcf

@@ -0,0 +1,92 @@
+"""
+get_map.py
+
+Generate a 4-column space-delimited map file accepted by selscan from ms or VCF input.
+
+Output columns:
+1) Chromosome name (default: "chr1")
+2) SNP serial number starting at 1
+3) Genetic position (coordinate used for integrating iHH decay in iHS)
+4) Physical position (bp)
+
+Notes:
+- Column 3: genetic position used for iHH integration (here equal to physical position).
+- Column 4: physical position used by selscan to determine gap thresholds for EHH truncation.
+"""
+
+import gzip
+
+def open_maybe_gz(file_path):
+    if file_path.endswith(".gz"):
+        return gzip.open(file_path, "rt")
+    else:
+        return open(file_path, "r")
+    
+
+import argparse
+import sys
+
+def parse_ms_positions(ms_file, genome_length):
+    physical_positions = []
+    with open_maybe_gz(ms_file) as f:
+        for line in f:
+            line = line.strip()
+            if line.startswith("positions:"):
+                parts = line.split()
+                rel_pos = parts[1:]  # skip 'positions:'
+                physical_positions = [int(float(p) * genome_length) for p in rel_pos]
+                break
+    if not physical_positions:
+        raise ValueError("No positions line found in ms file.")
+    return physical_positions
+
+def parse_vcf_positions(vcf_file):
+    physical_positions = []
+    with open_maybe_gz(vcf_file) as f:
+        for line in f:
+            if line.startswith("#"):
+                continue
+            parts = line.strip().split('\t')
+            pos = int(parts[1])
+            physical_positions.append(pos)
+    if not physical_positions:
+        raise ValueError("No variant positions found in VCF file.")
+    return physical_positions
+
+def write_map(positions, chrom="chr1", out_file="output.map"):
+    with open(out_file, "w") as f:
+        ## Write header line describing columns (space-delimited)
+        # f.write("chrom snp_index genetic_pos physical_pos\n")
+        for i, pos in enumerate(positions, start=0):
+            f.write(f"{chrom} rs{i} {pos} {pos}\n")
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Generate a 4-column (<chrom> <snp_index> <genetic_pos> <physical_pos>) space-delimited map file from ms or VCF input.",
+        epilog="Note: Col3=genetic pos for iHH integration (here equals physical pos, assuming constant recombination rate), Col4=physical pos for maximum allowable EHH gap detection."
+    )
+    group = parser.add_mutually_exclusive_group(required=True)
+    group.add_argument("--ms", help="Input ms file with positions line")
+    group.add_argument("--vcf", help="Input VCF file")
+
+    parser.add_argument("--len", type=int,
+                        help="Genome length for scaling relative positions in ms file (required if --ms is used)")
+    parser.add_argument("--chrom", default="chr1", help="Chromosome name to use in output (default: chr1)")
+    parser.add_argument("--out", default="output.map", help="Output map file name (default: output.map)")
+
+    args = parser.parse_args()
+
+    if args.ms and args.len is None:
+        print("Error: --len is required when using --ms", file=sys.stderr)
+        sys.exit(1)
+
+    if args.ms:
+        positions = parse_ms_positions(args.ms, args.len)
+    else:
+        positions = parse_vcf_positions(args.vcf)
+
+    write_map(positions, chrom=args.chrom, out_file=args.out)
+    print(f"Map file written to {args.out}")
+
+if __name__ == "__main__":
+    main()

scripts/colormap.plotting.R

@@ -0,0 +1,31 @@
+library(ggplot2)
+library(reshape2)
+
+# --- Read the data ---
+ehh_data <- read.table("outfile.ehh.1943.out", header = TRUE, sep = "\t", stringsAsFactors = FALSE)
+
+# Rename derEHH to derived for plotting
+colnames(ehh_data)[colnames(ehh_data) == "derEHH"] <- "derived"
+
+# Choose x-axis: pdist or gdist
+x_axis <- "pdist"  # Change to "gdist" if needed
+
+ehh_data$pdist <- ehh_data$pdist / 1e6
+
+# Reshape data to long format
+ehh_long <- melt(ehh_data, id.vars = x_axis, 
+                 measure.vars = c("derived", "ancEHH", "EHH"),
+                 variable.name = "Type", value.name = "EHH_Value")
+
+# Define custom colors
+ehh_colors <- c("derived" = "red", "ancEHH" = "blue", "EHH" = "gray50")
+
+# Plot
+ggplot(ehh_long, aes_string(x = x_axis, y = "EHH_Value", color = "Type")) +
+  geom_line(size = 1) +
+  scale_color_manual(values = ehh_colors) +
+  labs(x = ifelse(x_axis == "pdist", "Physical Distance (bp)", "Genetic Distance (cM)"),
+       y = "EHH",
+       title = "EHH Decay Plot") +
+  theme_minimal() +
+  theme(text = element_text(size = 14)) 

scripts/get_map.py

@@ -0,0 +1,82 @@
+# Makefile for selscan (macOS ARM) -
+
+# Compiler
+CC = g++ 	
+#clang++ and g++ are both supported
+
+# Optimization and architecture flags  and language standard setup
+G++FLAG = -O3 -m64 -ftree-vectorize -w -Wno-psabi -std=c++17 
+#G++FLAG = -DDEBUG -O0 -g -w -Wno-psabi -std=c++17 
+# -O3              : Maximum optimization level
+# -m64             : Generate 64-bit code
+# -arch arm64      : Target Apple Silicon (ARM64) architecture
+# -ftree-vectorize : Enable automatic vectorization of loops
+# -w               : Suppress all warnings
+# -Wno-psabi       : Disable warnings about ABI differences (useful for ARM cross-compilation)
+# -std=c++17       : Use the C++17 language standard
+
+
+# Include and library paths
+I_PATH = -I../include                     # Header file search path
+L_PATH = ../lib/macos-arm                 # Path to GSL libraries
+
+# Linking options
+LINK_OPTS_SELSCAN = -lpthread -lz         # Link with pthread and zlib for selscan
+LINK_OPTS_NORM = -L$(L_PATH) -lgsl -lgslcblas # Link with GSL and CBLAS for norm #LINK_OPTS_NORM = $(L_PATH)/libgsl.a $(L_PATH)/libgslcblas.a   #For static linking of norm program to gsl libs
+
+# Object file groups
+OBJ_STATS = xpihh.o ehh.o ihs.o ihh12.o pi.o ehh12.o
+OBJ_HAPMAP = hapdata.o mapdata.o bitset.o
+
+# Build all targets
+all: selscan norm2
+
+# Build selscan binary
+selscan: selscan-main.o $(OBJ_STATS) selscan-data.o $(OBJ_HAPMAP) selscan-cli.o binom.o param_t.o gzstream.o thread_pool.o
+	$(CC) $(G++FLAG) -o selscan $^ $(LINK_OPTS_SELSCAN)
+
+# Build norm binary
+norm2: norm2.o param_t.o
+	$(CC) $(G++FLAG) -o norm2 $^ $(LINK_OPTS_NORM)
+
+# Compilation rules
+selscan-main.o: selscan-main.cpp
+	$(CC) $(G++FLAG) -c $< $(I_PATH)
+
+selscan-data.o: selscan-data.cpp selscan-data.h
+	$(CC) $(G++FLAG) -c $< $(I_PATH)
+
+selscan-cli.o: selscan-cli.cpp selscan-cli.h
+	$(CC) $(G++FLAG) -c $< $(I_PATH)
+
+binom.o: binom.cpp binom.h
+	$(CC) $(G++FLAG) -c $<
+
+norm2.o: norm2.cpp
+	$(CC) $(G++FLAG) -c $< $(I_PATH)
+
+param_t.o: param_t.cpp param_t.h
+	$(CC) $(G++FLAG) -c $<
+
+gzstream.o: gzstream.cpp gzstream.h
+	$(CC) $(G++FLAG) -c $< $(I_PATH)
+
+thread_pool.o: thread_pool.cpp thread_pool.h
+	$(CC) $(G++FLAG) -c $< $(LINK_OPTS_SELSCAN)
+
+# Pattern rules for hapmap and stats modules
+%.o: hapmap/%.cpp hapmap/%.h
+	$(CC) $(G++FLAG) -c $< -o $@ $(I_PATH)
+
+%.o: stats/%.cpp stats/%.h
+	$(CC) $(G++FLAG) -c $< -o $@ $(I_PATH)
+
+# Optional: filetype modules (currently disabled)
+# %.o: filetype/%.cpp filetype/%.h
+# 	$(CC) $(G++FLAG) -c $< -o $@ $(I_PATH)
+
+# Cleanup rules
+.PHONY: clean
+clean:
+	rm -rf *.o selscan norm2 norm
+# rm -rf outfile*