- 讀取文件
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è)置每一個位點對應(yīng)的reference position
GENOME.class <- readData('VCF_chr18',format = 'VCF')
#### 或者這種方式讀取VCF
### 文件都要放在一個文件夾里床玻，這里的FASTA_chr18是一個文件夾固蚤，里面有一個fasta文件
- 定義sliding window联贩，計算
GENOME.class@n.sites
GENOME.class@region.names
get.sum.data(GENOME.class)
### 簡單看一些統(tǒng)計
GENOME.class <- set.populations(GENOME.class,list(
c(1:9),
c(10:17)
))
### 根據(jù)文件中的順序定義種群烈钞，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的意思是定義的長度包含非多態(tài)性位點，1000就是1000個核苷酸長度晋涣，甭管里面幾個SNPs仪媒；
#### type=1的意思是定義的長度只考慮多態(tài)性位點，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算吩，中性檢驗統(tǒng)計量
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，分別把種群內(nèi)的中性檢驗結(jié)果輸入對應(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
#### 重新標定一下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單獨拿出來構(gòu)建數(shù)據(jù)框佃扼，因為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)
#### 寫入數(shù)據(jù)
- 畫圖
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
