rehh计算xpehh-CSDN博客

1.结果说明

selscan-2.0.3计算xpehh的在另外一篇文章里面，计算的结果和rehh的结果比对，selscan标准化之后和rehh出的xpehh是一样的，然后ihs的结果是相反的，应该问题不大，反正我们需要的超过阈值的位点是一样的。

这个补充一点的是selscan计算ihs的时候出现的文件是这个，标准化的文件类似xpehh

2.rehh计算ihs和xpehh的流程

#########################
#install.packages("rehh")
#install.packages("ggplot2")
#install.packages("magrittr")
#install.packages("tidyverse")
#install.packages("CMplot")
library(rehh)
library(magrittr)
library(dplyr)
library(CMplot)
#install.packages("tidyverse")
library(tidyverse)

hh_map1 <- data2haplohh(
  hap_file = "2g_gwas.phased17.vcf",
  polarize_vcf = FALSE,#设置极化标记方向，确定等位基因的方向性
  vcf_reader = "data.table",
  chr.name = 17,
  map_file = "2g_new17.map",
  allele_coding = "map")


hh_map2 <- data2haplohh(
  hap_file = "8g_gwas.phased17.vcf",
  vcf_reader = "data.table",
  polarize_vcf = FALSE,
  map_file = "CHR17/8g_new17.map",
  allele_coding = "map")
#函数scan_hh（）为haplohh对象中的所有标记计算iHH\iES值。
#如果标记未极化，则polarized选项应设置为FALSE，从而导致该功能使用频率最高和频率第二高等位基因

res.scan1 <- scan_hh(hh_map1)
res.scan2 <- scan_hh(hh_map2)

res.ihs1<- ihh2ihs(res.scan1)
res.ihs2<- ihh2ihs(res.scan2)

write.table(res.ihs1$ihs,"g2.chr17.ihs.txt",quote=F,row.names=F)
write.table(res.ihs2$ihs,"g8.chr17.ihs.txt",quote=F,row.names=F)

### plot the iHS statistics
ggplot(res.ihs1$ihs, aes(POSITION, IHS)) + geom_point()+  ggtitle("g2_POSITION_ihs_r包")
ggplot(res.ihs2$ihs, aes(POSITION, IHS)) + geom_point()+  ggtitle("g8_POSITION_ihs_r包")
### plot the log P-value 
ggplot(res.ihs1$ihs, aes(POSITION, LOGPVALUE)) + geom_point()


#### 这里同时跑了rehh的xpehh和rsb两种方法，但其实差不多，后者对core标准化，而前者不
res.xpehh <- ies2xpehh(res.scan1, res.scan2, "g2", "g8")
res.rsb<-ines2rsb(res.scan1, res.scan2, "g2", "g8")

#PLOT the xpEHH values
ggplot(res.xpehh, aes(POSITION, XPEHH_g2_g8)) + geom_point() + ggtitle("g2_g8_POSITION_XPEHH_r包")


### 把结果输入到文件
write.table(res.xpehh,"g2_g8.chr17.xpehh.txt",quote=F,row.names=F)
write.table(res.rsb,"g2_g8.chr17.rsb.txt",quote=F,row.names=F)
############################这个是画出selscan的图

input1=read.table("2g_ihs.ihs.out",header = T)
input2=read.table("8g_ihs.ihs.out",header = T)

### plot the iHS statistics
ggplot(input1, aes(Pos, unstandardized_iHS)) + geom_point()+  ggtitle("g2_POSITION_ihs")
ggplot(input2, aes(Pos, unstandardized_iHS)) + geom_point()+  ggtitle("g8_POSITION_ihs")

input3=read.table(“chr17.8_2.xpehh.out",header = T)
#PLOT the xpEHH values
ggplot(input3, aes(pos,xpehh)) + geom_point() + ggtitle("g2_g8_POSITION_XPEHH")