IBD_3500=$1
NormRMSE=$2
P=$3

paste -d" " <(awk '{print $0}' $IBD_3500) <(awk '{print $4}' $NormRMSE) > temp
awk '{print $1, $2, $3, $4, $5, $6, $7, $8, (1-$6 + $7 + 1-$8)/3}' temp > temp2
sed -r '1s/\S+//9' temp2 > temp
sed -i '1s/$/score/' temp
# Look only at haploid
mv temp ${P}_GL_IBD_3500kb_onlyHap.txt; rm temp2

head -n 1 ${P}_GL_IBD_3500kb_onlyHap.txt > header.txt
sed '1d' ${P}_GL_IBD_3500kb_onlyHap.txt > data.txt

sort -nk9 data.txt > sorted_data.txt ; rm data.txt
cat header.txt sorted_data.txt > ${P}_GL_IBD_3500kb_onlyHap.txt ; rm header.txt sorted_data.txt

cat ${P}_GL_IBD_3500kb_onlyHap.txt

PARAMS=`tail -n 1 ${P}_GL_IBD_3500kb_onlyHap.txt | cut -f 2-9 -d" "`
PARAMS=(${PARAMS//,/ })
echo $PARAMS
BITS=${PARAMS[1]}
ERRHOM=${PARAMS[2]}
PERC_DIFF=${PARAMS[4]}
PERC_MATCH=${PARAMS[5]}
KIN_RESIDS=${PARAMS[6]}
SCORE=${PARAMS[7]}

echo -e $P GERMLINE2 parameters"\n"Type: ${TYPE}"\n"Bits: ${BITS}"\n"Err_hom: ${ERRHOM}"\n"Percent difference: ${PERC_DIFF}"\n"Percent match: ${PERC_MATCH}"\n"Kin residuals: ${KIN_RESIDS}"\n"Score: ${SCORE} > results/other/${P}_GL_parameters.txt
# awk -v Pop="$P" -v IGNORECASE=1 '$1 ~ Pop && $3 ~ Pop {print}' RoH_param_combs/${BITS}/IBD_${BITS}_${ERRHOM}_hap_segsJoined_QCed.ibd | cut -f1-7,9 > ${P}_only_QCed.ibd


cat results/RoH_param_combs/${BITS}/IBD_${BITS}_${ERRHOM}_hap_segsJoined_QCed.ibd > results/other/${P}_only_QCed.ibd

args = commandArgs(trailingOnly=TRUE)

# Defaults
perc.overlap <- 0.95

# Error handling and updating default values
if (length(args)<2) {
  stop("At least two arguments must be supplied", call.=FALSE)
} else if (length(args)==3) {
  perc.overlap <- as.numeric(args[3])
} else if (length(args)>3) {
  stop("Too many arguments", call.=FALSE)
}

options(stringsAsFactors=F)
fileNameSplit <- strsplit(args[1], split=".", fixed=T)[[1]]
ext.ind <- which(fileNameSplit %in% c("match", "hom", "ibd"))
stem <- paste(fileNameSplit[-ext.ind], collapse=".")
ext <- fileNameSplit[ext.ind]

if (ext=="hom") {
  seg.file <- read.table(args[1], header=T)
  colnames(seg.file)[7:8] <- c("start", "end")
} else if (ext=="match") {
  seg.file <- read.table(args[1])
  if (ncol(seg.file)==13) {
    colnames(seg.file)[3:5] <- c("CHR", "start", "end")
  } else if (ncol(seg.file)==15) {
    colnames(seg.file)[5:7] <- c("CHR", "start", "end")
  }
} else if (ext=="ibd") {
  seg.file <- read.table(args[1])
  colnames(seg.file)[5:7] <- c("CHR", "start", "end")
}

excl.regs <- read.table(args[2])
  colnames(excl.regs) <- c("CHR", "start", "end", "type")

suppressMessages(library("GenomicRanges"))
Excl.Regs <- makeGRangesFromDataFrame(excl.regs[,c("CHR", "start", "end")], seqnames.field="CHR")
Seg.File <- makeGRangesFromDataFrame(seg.file[ ,c("CHR", "start", "end")], seqnames.field="CHR")

Hits <- findOverlaps(Seg.File, Excl.Regs) # find all hits
Overlaps <- pintersect(Seg.File[queryHits(Hits)], Excl.Regs[subjectHits(Hits)])
percentOverlap <- width(Overlaps)/width(Seg.File[queryHits(Hits)])
HighOverlapHits <- Hits[percentOverlap > perc.overlap] # Which have high overlap? Some will be removed

new.seg.file <- seg.file[-queryHits(HighOverlapHits),]

if (ext=="hom") {
  COLNAMES <- c(colnames(new.seg.file)[1:6], "POS1", "POS2", colnames(new.seg.file)[9:13])
  write.table(new.seg.file, "", row.names=F, col.names=COLNAMES, quote=F, sep="\t")
} else {
  write.table(new.seg.file, "", row.names=F, col.names=F, quote=F, sep="\t")
}

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!requireNamespace("GenomicRanges", quietly = TRUE))
    BiocManager::install("GenomicRanges")


#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)

# Defaults
perc.thresh <- 0.05
perc.overlap <- 0.8

# Error handling and updating default values
if (length(args)<3) {
  stop("At least three arguments must be supplied", call.=FALSE)
} else if (length(args)==4) {
  perc.thresh <- as.numeric(args[4])
} else if (length(args)==5) {
  perc.thresh <- as.numeric(args[4])
  perc.overlap <- as.numeric(args[5])
} else if (length(args)>5) {
  stop("Too many arguments", call.=FALSE)
}

options(stringsAsFactors=F)
Bim <- read.table(args[1])
WinSize <- as.numeric(args[3])

# Determine type of input segment file
cat("Determining input file type... ")
fileNameSplit <- strsplit(args[2], split=".", fixed=T)[[1]]
ext.ind <- match("hom", fileNameSplit)

if (length(fileNameSplit) > 1 & !is.na(ext.ind)) {
  stem <- paste(fileNameSplit[-ext.ind], collapse=".")
  ext <- "hom"
} else {
  stem <- args[2]
  ext="garlic"
}

if (ext=="hom") {
  Seg <- read.table(args[2], header=T)
  if (ncol(Seg)==13) {
    fileType <- "PLINK"
    colnames(Seg)[9] <- "LENGTH"
  }
  cat(fileType, "file type detected\n")
} else {
  Seg <- read.table(args[2], skip=1)
  Seg <- cbind(Seg, rep("placeholder_ID", nrow(Seg)))
  if (ncol(Seg)==10) {
    fileType <- "GARLIC"
    cat("WARNING: GARLIC file type detected - header line dropped\n")
    colnames(Seg)[c(1,2,3,4,5,10)] <- c("CHR","POS1","POS2","CLASS","LENGTH","IID")
    for (i in 1:22) {
      Seg[which(Seg$CHR==paste("chr", i, sep="")), "CHR"] <- i
    }
    Seg$CHR <- as.numeric(Seg$CHR)
    class.ranges <- list(); length(class.ranges) <- 3; names(class.ranges) <- c("A","B","C")
    for (class in c("A","B","C")) {
      class.ranges[[class]] <- range(as.numeric(Seg[which(Seg$CLASS==class), "LENGTH"]))
    }
    class.ranges[["A"]][2] <- mean(class.ranges[["A"]][2], class.ranges[["B"]][1])
    class.ranges[["B"]][1] <- mean(class.ranges[["A"]][2], class.ranges[["B"]][1]) + 0.5
    class.ranges[["B"]][2] <- mean(class.ranges[["B"]][2], class.ranges[["C"]][1]) 
    class.ranges[["C"]][1] <- mean(class.ranges[["B"]][2], class.ranges[["C"]][1]) + 0.5
  }
}

# Create copy of input file - eventual output file
Seg.new <- Seg

genome.SNP.count <- c()
chr.hists <- list()
  length(chr.hists) <- length(unique(Bim$V1))
  names(chr.hists) <- unique(Bim$V1)

for (chr in unique(Bim$V1)) {
  chr.data <- Bim[Bim$V1==chr,]
  start <- chr.data[1, "V4"]
  end <- chr.data[nrow(chr.data), "V4"]
  chr.hists[[as.character(chr)]] <- hist(chr.data$V4, breaks=seq(start, end+WinSize, WinSize), plot=F)
  genome.SNP.count <- c(genome.SNP.count, chr.hists[[as.character(chr)]]$counts)
}
min.SNPs <- max(sort(genome.SNP.count)[length(genome.SNP.count)*perc.thresh], 1/5e4*WinSize)

suppressMessages(library("GenomicRanges"))
for (chr in sort(unique(Seg$CHR))) {
  cat("Starting chr", chr, "\n")
  # Calculate SNP density across the chromosome in windows
  Hist <- chr.hists[[as.character(chr)]]
  Counts <- Hist$counts
  Breaks <- Hist$breaks
  # Store count information per region
  if (fileType=="PLINK") {
    chrWinInfo <- as.data.frame(matrix(ncol=4, nrow=length(Counts)))
      colnames(chrWinInfo) <- c("chr", "start", "end", "nSNPs")
  } else if (fileType=="GARLIC") {
    chrWinInfo <- as.data.frame(matrix(ncol=3, nrow=length(Counts)))
      colnames(chrWinInfo) <- c("chr", "start", "end")
  }
  chrWinInfo$chr <- chr
  chrWinInfo$start <- Breaks[1:length(Counts)]
  chrWinInfo$end <- Breaks[2:length(Breaks)]
  chrWinInfo$nSNPs <- Counts
  # Find low density regions
  low.dens.wins <- chrWinInfo[which(chrWinInfo$nSNPs < min.SNPs), 1:3]
  # If low density regions are found join continguous windows into regions
  if (nrow(low.dens.wins) > 0) {
    Low.Density.Regions <- low.dens.wins[1, ]
    if (nrow(low.dens.wins) > 1) {
      for (x in 2:nrow(low.dens.wins)) {
        if (low.dens.wins[x, "start"] == Low.Density.Regions[nrow(Low.Density.Regions), "end"]) {
          Low.Density.Regions[nrow(Low.Density.Regions), "end"] <- low.dens.wins[x, "end"]
        } else {
          Low.Density.Regions[nrow(Low.Density.Regions)+1, ] <- low.dens.wins[x, ]
        }
      }
    }
    Low.Density.Ranges <- makeGRangesFromDataFrame(Low.Density.Regions)
    # Join segments in input file
    for (Ind in unique(Seg$IID)) {
      chrSegInd <- Seg[Seg$CHR==chr&Seg$IID==Ind,]
        chrSegInd <- chrSegInd[order(chrSegInd$POS1),]
      if (nrow(chrSegInd)>1) {
        RoH.gaps <- as.data.frame(matrix(ncol=3,nrow=0))
          colnames(RoH.gaps) <- c("chr","start","end")
        for (y in 2:max(2, nrow(chrSegInd))) {
          RoH.gaps[y-1, ] <- c(chr, chrSegInd[y-1, "POS2"], chrSegInd[y, "POS1"])
        }
        RoH.gaps <- RoH.gaps[which(RoH.gaps$end > RoH.gaps$start),]
        if (nrow(RoH.gaps) > 0) {
          RoH.gap.Ranges <- makeGRangesFromDataFrame(RoH.gaps)
          # Do any gaps overlap a low density region?
          Hits <- findOverlaps(RoH.gap.Ranges, Low.Density.Ranges) # find all hits
          Overlaps <- pintersect(Low.Density.Ranges[subjectHits(Hits)], RoH.gap.Ranges[queryHits(Hits)])
          percentOverlap <- width(Overlaps) / width(RoH.gap.Ranges[queryHits(Hits)]) # calculate overlap
          Hits <- Hits[percentOverlap > perc.overlap] # sufficient overlap?
          joinOverGaps <- RoH.gaps[queryHits(Hits),]
          if (nrow(joinOverGaps)>0) {
            for (z in 1:nrow(joinOverGaps)) {
              seg1 <- rownames(chrSegInd[which(chrSegInd$POS2 == joinOverGaps[z, "start"]),])[1]
              seg2 <- rownames(chrSegInd[which(chrSegInd$POS1 == joinOverGaps[z, "end"]),])[1]
              joinedSeg <- chrSegInd[as.character(seg1),]
              joinedSeg$POS2 <- chrSegInd[as.character(seg2),"POS2"]
              joinedSeg$LENGTH <- chrSegInd[as.character(seg1),"LENGTH"] + chrSegInd[as.character(seg2),"LENGTH"]
              if (fileType=="PLINK") {
                joinedSeg$SNP2 <- chrSegInd[as.character(seg2),"SNP2"]
                joinedSeg$NSNP <- chrSegInd[as.character(seg1),"NSNP"] + chrSegInd[as.character(seg2),"NSNP"]
              } else if (fileType=="GARLIC") {
                if (joinedSeg$LENGTH < class.ranges[["A"]][1]) {
                  joinedSeg$CLASS <- "A"
                } else if (joinedSeg$LENGTH >= class.ranges[["A"]][1] & joinedSeg$LENGTH <= class.ranges[["A"]][2]) {
                  joinedSeg$CLASS <- "A"
                } else if (joinedSeg$LENGTH >= class.ranges[["B"]][1] & joinedSeg$LENGTH <= class.ranges[["B"]][2]) {
                  joinedSeg$CLASS <- "B"
                } else if (joinedSeg$LENGTH >= class.ranges[["C"]][1] & joinedSeg$LENGTH <= class.ranges[["C"]][2]) {
                  joinedSeg$CLASS <- "C"
                } else if (joinedSeg$LENGTH > class.ranges[["C"]][2]) {
                  joinedSeg$CLASS <- "C"
                }
              }
              Seg.new <- rbind(Seg.new, joinedSeg)
              # Remove the two split segments from the Seg.joined object
              Seg.new <- Seg.new[-match(as.character(seg1), rownames(Seg.new)),]
              Seg.new <- Seg.new[-match(as.character(seg2), rownames(Seg.new)),]
            }
          }
        }
      }
    }
  }
}

if (fileType=="PLINK") {
  colnames(Seg.new)[9] <- "KB"
  write.table(Seg.new, paste(stem, "_segsJoined.hom", sep=""), quote=F, row.names=F, col.names=T)
} else if (fileType=="GARLIC") {
  Seg.new <- Seg.new[,-ncol(Seg.new)]
  Seg.new$CHR <- paste("chr", Seg.new$CHR, sep="")
  write.table(Seg.new, paste(stem, "_segsJoined.garlic", sep=""), quote=F, row.names=F, col.names=F, sep="\t")
}

INPUT=$1
MAP=$2
LOG=$3
OUTPUT=$4
REF_HAP=$5
REF_LEG=$6
REF_SAMP=$7

# ShapeIT phasing and conversion to plink files
#----------------------------------------------
# Divide plink file by chromosome, trim SNPs that lie outside the boundaries of the genetic map
# Divide plink file by chromosome, trim SNPs that lie outside the boundaries of the genetic map

Rscript scripts/trim_bimfile.R ${INPUT} ${MAP}
plink --bfile ${INPUT} --exclude ${INPUT}_removeSNPs --make-bed --out ${INPUT}_endsTrimmed
#cat ${FILE}.chr${chr}_endsTrimmed.bim >> dat_QCed_endsTrimmed.bim

# Run SHAPEIT
shapeit -B ${INPUT}_endsTrimmed -M ${MAP} --input-ref ${REF_HAP} ${REF_LEG} ${REF_SAMP} --duohmm -W 5 -O ${OUTPUT} --output-log ${LOG}.log -T 4

files=( ${LOG}'.snp.strand' )
if (( ${#files[@]} )); then
  N=`ls -l ${LOG}.snp.strand | wc -l`
  while [ $N -gt 0 ]
  do
      grep Missing ${LOG}.snp.strand | cut -f4 > ${INPUT}_missing_from_ref.txt
      grep Strand ${LOG}.snp.strand | awk '$11==1 {print $4}' > ${INPUT}_flip_in_dataset.txt
      grep Strand ${LOG}.snp.strand | awk '$11==0 {print $4}' > ${INPUT}_cannot_flip_in_dataset.txt
      cat ${INPUT}_missing_from_ref.txt ${INPUT}_cannot_flip_in_dataset.txt > ${INPUT}_remove_from_dataset.txt
      rm -f ${LOG}.snp.strand*
      rm -f ${INPUT}_missing_from_ref.txt ${INPUT}_cannot_flip_in_dataset.txt
      plink --bfile ${INPUT}_endsTrimmed --flip ${INPUT}_flip_in_dataset.txt --exclude ${INPUT}_remove_from_dataset.txt --make-bed --out ${INPUT}_endsTrimmed_misalignRepaired
      shapeit -B ${INPUT}_endsTrimmed_misalignRepaired -M ${MAP} --input-ref ${REF_HAP} ${REF_LEG} ${REF_SAMP} --duohmm -W 5 -O ${OUTPUT} --output-log ${LOG}.log -T 4
      N=`ls -l ${LOG}.snp.strand | wc -l`
  done
fi

from numpy import *
from scipy.interpolate import interp1d
from sys import argv

# use, e.g.:
# python shapeit_to_germline.py shapeitoutputprefix geneticmapfile
# python shapeit_to_germline.py Ethiopians_chr21_SHAPEITphased_duoHMM genetic_map_chr21_combined_b37.txt

prefix = argv[-2]
hapsfile = prefix + '.haps'
samplefile = prefix + '.sample'

hapmapfile = argv[-1]

famdata = [line.strip().split() for line in file(samplefile)][2:]

# read in haplotypes
print 'reading in haps data...'
hapsdata = [line.strip().split() for line in file(hapsfile)]
print 'splitting haps data'
mapinfo = array([line[:3] for line in hapsdata])
variants = array([line[3:5] for line in hapsdata])
alleles = array([line[5:] for line in hapsdata],dtype=int)


#hapsdata = array([line.strip().split()[:3] for line in file(hapsfile)])
bps = array( mapinfo[:,2] , dtype=int)

# read in hapmap file
mapdata = array([line.strip().split() for line in file(hapmapfile)][2:],dtype=float)
happos = mapdata[:,0]
hapcms = mapdata[:,2]

# now create interpolation function
f = interp1d(happos,hapcms,fill_value='extrapolate')

bpcms = f(bps)

# make map file
outmap = file(prefix + '.map' , 'w')
print 'writing mapfile to %s' % (outmap)
for i, line in enumerate(mapinfo):
	outmap.write('%s\t%s\t%s\t%s\n' % (line[0],line[1],bpcms[i],line[2]))


# make ped file. More annoying for sure!
# brute force for now, hopefully it runs okay (not huge files)
outped = file(prefix + '.ped' , 'w')
print 'writing pedfile to %s' % (outped)
alleles_acgt = array([variants[i][alleles[i]] for i in range(len(variants))])
for i, line in enumerate(famdata):
	if i % 20 == 0:
		print 'on individual %s' % (i)
	# zip together haplotypes 2*i and 2*i+1 so they're on the same line
	plinkalleles = array(ndarray.flatten(array(zip(alleles_acgt[:,2*i],alleles_acgt[:,2*i+1] ) ) ), dtype='S1')
	outped.write('%s %s\n' % ( ' '.join(line[:2] + line[3:]) , ' '.join(plinkalleles)))

outmap.close()
outped.close()

print 'done'