## ＰＴＴ大神相助 1 #########################
##########################################
# andrew43 建議approx()解決
##########################################


library(ggplot2)

loci <- read.table("loci list", header=TRUE)
ci <- read.table("ci_add0", header=TRUE)


dt <- data.frame(id = loci$Locus, x = loci$Het, y = loci$Fst)
line.upper <- data.frame(x = ci$Het, y = ci$X0.5.1.pval.quantile)
line.lower <- data.frame(x = ci$Het, y = ci$Fst)

dt$upper <- approx(line.upper$x, line.upper$y, dt$x)$y
dt$lower <- approx(line.lower$x, line.lower$y, dt$x)$y
dt$code <- dt$y > dt$lower & dt$y < dt$upper

print(dt)

ggplot(dt, aes(x, y)) +
  geom_point(aes(color = code)) +
  geom_line(data = line.upper, color = "red", aes(x = x, y = y)) +
  geom_line(data = line.lower, color = "green", aes(x = x, y = y))
		#原作者andrew43提醒
		# 紅點表示在界外（dt$code值為FALSE），
		# 綠點表示在界內（dt$code值為TRUE），
		# 灰點表示不能判斷（上界或下界不存在；dt$code值為NA）

# 抓取neutral loci，以yes為抓取基準
neutral <- subset(dt, code=="TRUE")
outlier <- subset(dt, code=="FALSE")

# 計算數量(多種方式)
dim(outlier)[1]
nrow(outlier)
length(outlier[,1])

# 輸出
write.table(outlier$id, file="/Users/path/list.txt", quote = FALSE, row.names = FALSE, col.names = FALSE)




## ＰＴＴ大神相助 2 #########################
##########################################
# showfeb 建議 sp::point.in.polygon 
##########################################
	
require(ggplot2)
require(sp)

loci <- read.table("loci list edit", header=TRUE)
ci <- read.table("ci_add0", header=TRUE)


df=data.frame(id=loci$Locus, x=loci$Het, y=loci$Fst)          

Line1=data.frame(x=ci$Het, y=ci$Fst, g="Lower")
Line2=data.frame(x=rev(ci$Het),y=rev(ci$X0.5.1.pval.quantile),g="Upper")
Linedf=rbind(Line1,Line2)

df$in.pg = ifelse(as.logical(point.in.polygon(df$x, df$y, pol.x=Linedf$x, pol.y=Linedf$y))
                  ,"Yes","No")
head(df)

ggplot()+
  geom_point(aes(x=x,y=y,colour=in.pg),data=df,inherit.aes = F)+
  geom_path(aes(x=x,y=y,group=g,colour=g),size=1.2,data=Linedf,inherit.aes = F)
		# 注意多邊形繪製的順序，尤其Line2當中的x與y有"rev"，倒序後才可以繪製成完整的多邊形
		#  可使用下列方式輸出多邊形確認是否正確 
		#	plot.new()
		#	polygon(Linedf$x, Linedf$y)

# 抓取neutral loci，以yes為抓取基準
neutral <- subset(df, in.pg=="Yes")
outlier <- subset(df, in.pg=="No")

# 計算數量(多種方式)
dim(outlier)[1]
nrow(outlier)
length(outlier[,1])

# 輸出
write.table(outlier$id, file="/Users/path/list.txt", quote = FALSE, row.names = FALSE, col.names = FALSE)