fork download
copy 
  1. ## ＰＴＴ大神相助 1 #########################
  2. ##########################################
  3. # andrew43 建議approx()解決
  4. ##########################################
  5.  
  6.  
  7. library(ggplot2)
  8.  
  9. loci <- read.table("loci list", header=TRUE)
  10. ci <- read.table("ci_add0", header=TRUE)
  11.  
  12.  
  13. dt <- data.frame(id = loci$Locus, x = loci$Het, y = loci$Fst)
  14. line.upper <- data.frame(x = ci$Het, y = ci$X0.5.1.pval.quantile)
  15. line.lower <- data.frame(x = ci$Het, y = ci$Fst)
  16.  
  17. dt$upper <- approx(line.upper$x, line.upper$y, dt$x)$y
  18. dt$lower <- approx(line.lower$x, line.lower$y, dt$x)$y
  19. dt$code <- dt$y > dt$lower & dt$y < dt$upper
  20.  
  21. print(dt)
  22.  
  23. ggplot(dt, aes(x, y)) +
  24.   geom_point(aes(color = code)) +
  25.   geom_line(data = line.upper, color = "red", aes(x = x, y = y)) +
  26.   geom_line(data = line.lower, color = "green", aes(x = x, y = y))
  27. 		#原作者andrew43提醒
  28. 		# 紅點表示在界外（dt$code值為FALSE），
  29. 		# 綠點表示在界內（dt$code值為TRUE），
  30. 		# 灰點表示不能判斷（上界或下界不存在；dt$code值為NA）
  31.  
  32. # 抓取neutral loci，以yes為抓取基準
  33. neutral <- subset(dt, code=="TRUE")
  34. outlier <- subset(dt, code=="FALSE")
  35.  
  36. # 計算數量(多種方式)
  37. dim(outlier)[1]
  38. nrow(outlier)
  39. length(outlier[,1])
  40.  
  41. # 輸出
  42. write.table(outlier$id, file="/Users/path/list.txt", quote = FALSE, row.names = FALSE, col.names = FALSE)
  43.  
  44.  
  45.  
  46.  
  47. ## ＰＴＴ大神相助 2 #########################
  48. ##########################################
  49. # showfeb 建議 sp::point.in.polygon 
  50. ##########################################
  51.  
  52. require(ggplot2)
  53. require(sp)
  54.  
  55. loci <- read.table("loci list edit", header=TRUE)
  56. ci <- read.table("ci_add0", header=TRUE)
  57.  
  58.  
  59. df=data.frame(id=loci$Locus, x=loci$Het, y=loci$Fst)          
  60.  
  61. Line1=data.frame(x=ci$Het, y=ci$Fst, g="Lower")
  62. Line2=data.frame(x=rev(ci$Het),y=rev(ci$X0.5.1.pval.quantile),g="Upper")
  63. Linedf=rbind(Line1,Line2)
  64.  
  65. df$in.pg = ifelse(as.logical(point.in.polygon(df$x, df$y, pol.x=Linedf$x, pol.y=Linedf$y))
  66.                   ,"Yes","No")
  67. head(df)
  68.  
  69. ggplot()+
  70.   geom_point(aes(x=x,y=y,colour=in.pg),data=df,inherit.aes = F)+
  71.   geom_path(aes(x=x,y=y,group=g,colour=g),size=1.2,data=Linedf,inherit.aes = F)
  72. 		# 注意多邊形繪製的順序，尤其Line2當中的x與y有"rev"，倒序後才可以繪製成完整的多邊形
  73. 		#  可使用下列方式輸出多邊形確認是否正確 
  74. 		#	plot.new()
  75. 		#	polygon(Linedf$x, Linedf$y)
  76.  
  77. # 抓取neutral loci，以yes為抓取基準
  78. neutral <- subset(df, in.pg=="Yes")
  79. outlier <- subset(df, in.pg=="No")
  80.  
  81. # 計算數量(多種方式)
  82. dim(outlier)[1]
  83. nrow(outlier)
  84. length(outlier[,1])
  85.  
  86. # 輸出
  87. write.table(outlier$id, file="/Users/path/list.txt", quote = FALSE, row.names = FALSE, col.names = FALSE)
Success #stdin #stdout #stderr 0.41s 50564KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
Error in file(file, "rt") : cannot open the connection
Calls: read.table -> file
In addition: Warning message:
In file(file, "rt") :
  cannot open file 'loci list': No such file or directory
Execution halted
https://ideone.com/SmNNEv
language:
R (R 3.5.2)
created:
2 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!