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))
          ##### 并設(shè)置每一個(gè)位點(diǎn)對(duì)應(yīng)的reference position

GENOME.class <- readData('VCF_chr18',format = 'VCF')
             #### 或者這種方式讀取VCF


### 文件都要放在一個(gè)文件夾里，這里的FASTA_chr18是一個(gè)文件夾，里面有一個(gè)fasta文件

2. 定義sliding window，計(jì)算

GENOME.class@n.sites
GENOME.class@region.names
get.sum.data(GENOME.class)
### 簡(jiǎn)單看一些統(tǒng)計(jì)

GENOME.class <- set.populations(GENOME.class,list(
             c(1:9),
            c(10:17)
           ))
### 根據(jù)文件中的順序定義種群，1-9為一個(gè)，10-17為一個(gè)

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


GENOME.class.slide  <- sliding.window.transform(GENOME.class,width=1000,jump=1,type=2,whole.data=TRUE)
#### 使用window 長(zhǎng)度1000，每次step為1，
#### type2的意思是定義的長(zhǎng)度包含非多態(tài)性位點(diǎn)，1000就是1000個(gè)核苷酸長(zhǎng)度，甭管里面幾個(gè)SNPs；
#### type=1的意思是定義的長(zhǎng)度只考慮多態(tài)性位點(diǎn)，1000就是1000個(gè)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
### 計(jì)算Fst，中性檢驗(yàn)統(tǒng)計(jì)量

pop1_out_df<-as.data.frame(get.neutrality(GENOME.class.slide)[[1]])
pop2_out_df<-as.data.frame(get.neutrality(GENOME.class.slide)[[2]])
### 我這里兩個(gè)population，分別把種群內(nèi)的中性檢驗(yàn)結(jié)果輸入對(duì)應(yīng)數(shù)據(jù)框

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
#### 重新標(biāo)定一下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
}
#### 重新標(biāo)定一下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單獨(dú)拿出來(lái)構(gòu)建數(shù)據(jù)框，因?yàn)镕st只有種群間的

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)
#### 寫入數(shù)據(jù)

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