PopGenome包计算Tajima's D, Fst, Fu

1. 读取文件

library(ape)
library(PopGenome)
# dir_of_gzvcffile=
chr='chr18'
GENOME.class <- readData('FASTA_chr18',format = 'fasta',
                         SNP.DATA=T) ### 这种方式读取fasta格式的SNPs
pos_data<-read.csv('chr17_W_data.fasta.pos')$Position
GENOME.class <-set.ref.positions(GENOME.class,list(pos_data))
          ##### 并设置每一个位点对应的reference position

GENOME.class <- readData('VCF_chr18',format = 'VCF')
             #### 或者这种方式读取VCF


### 文件都要放在一个文件夹里，这里的FASTA_chr18是一个文件夹，里面有一个fasta文件

2. 定义sliding window，计算

GENOME.class@n.sites
GENOME.class@region.names
get.sum.data(GENOME.class)
### 简单看一些统计

GENOME.class <- set.populations(GENOME.class,list(
             c(1:9),
            c(10:17)
           ))
### 根据文件中的顺序定义种群，1-9为一个，10-17为一个

GENOME.class@region.data@populations
### 看下种群信息


GENOME.class.slide  <- sliding.window.transform(GENOME.class,width=1000,jump=1,type=2,whole.data=TRUE)
#### 使用window 长度1000，每次step为1，
#### type2的意思是定义的长度包含非多态性位点，1000就是1000个核苷酸长度，甭管里面几个SNPs；
#### type=1的意思是定义的长度只考虑多态性位点，1000就是1000个SNPs。


GENOME.class.slide@region.names
GENOME.class.slide   <- neutrality.stats(GENOME.class.slide)
GENOME.class.slide  <- F_ST.stats(GENOME.class.slide)
GENOME.class.slide  <- diversity.stats(GENOME.class.slide)
GENOME.class.slide  <- detail.stats(GENOME.class.slide)
GENOME.class.slide@n.sites
### 计算Fst，中性检验统计量

pop1_out_df<-as.data.frame(get.neutrality(GENOME.class.slide)[[1]])
pop2_out_df<-as.data.frame(get.neutrality(GENOME.class.slide)[[2]])
### 我这里两个population，分别把种群内的中性检验结果输入对应数据框

row<-rownames(pop1_out_df)
pop1_out_df$start<-NA
pop1_out_df$end<-NA
pop1_out_df$mid_pos<-NA
pop2_out_df$start<-NA
pop2_out_df$end<-NA
pop2_out_df$mid_pos<-NA
#### 重新标定一下window的位置


rowcount=0
for (i in row){
    rowcount=rowcount+1
    start=as.numeric(strsplit(i,split = " - ")[[1]][1])
    end=as.numeric(strsplit(strsplit(i,split = " - ")[[1]][2],split=' :')[[1]][1])
    print(end)
    mid=(start+end)/2
    pop1_out_df$start[rowcount]<-start
    pop1_out_df$end[rowcount]<-end
    pop1_out_df$mid_pos[rowcount]<-mid
    pop2_out_df$start[rowcount]<-start
    pop2_out_df$end[rowcount]<-end
    pop2_out_df$mid_pos[rowcount]<-mid
}
#### 重新标定一下window的位置


fst<-data.frame(chr=chr,start=pop1_out_df$start,end=pop1_out_df$end,mid_pos=pop1_out_df$mid_pos,fst=as.data.frame(get.F_ST(GENOME.class.slide))[[2]])
pop1_out_df$chr=chr
pop2_out_df$chr=chr
#### fst单独拿出来构建数据框，因为Fst只有种群间的

write.csv(pop1_out_df,paste0(chr,'.1.stats.csv'),row.names=F,quote=F)
write.csv(pop2_out_df,paste0(chr,'.2.stats.csv'),row.names=F,quote=F)
#### 写入数据

3. 画图

pdf(paste0(chr,'.selection.stats.pdf'))
library(ggplot2)
p<-ggplot()+geom_point(data=pop1_out_df,aes(x=mid_pos,y=Tajima.D,color="pop_1"))+
    geom_point(data=pop2_out_df,aes(x=mid_pos,y=Tajima.D,color="pop_2"))
print(p)
p<-ggplot()+geom_point(data=pop1_out_df,aes(x=mid_pos,y=Fu.Li.F,color="pop_1"))+
    geom_point(data=pop2_out_df,aes(x=mid_pos,y=Fu.Li.F,color="pop_2"))
print(p)
p<-ggplot()+geom_point(data=pop1_out_df,aes(x=mid_pos,y=Fu.Li.D,color="pop_1"))+
    geom_point(data=pop2_out_df,aes(x=mid_pos,y=Fu.Li.D,color="pop_2"))
print(p)
p<-ggplot()+geom_point(data=fst,aes(x=mid_pos,y=fst))
print(p)
dev.off()
#### 画图

image.png image.png image.png