參考
命令
vcftools --vcf snp.bialles.vcf --SNPdensity 100000 --out StatResults/SNPdensity
100000 是指定窗口長(zhǎng)度
--out 是輸出文件的前綴
使用R語(yǔ)言中的circlize包畫(huà)圖
參考
代碼
df<-read.table("SNPdensity.snpden",sep="\t",header=T)
head(df)
df<-df[,c(1,2,4)]
colnames(df)<-c("Chr","X","Y")
head(df)
df$X<-df$X/1000000
options(scipen=999)
library(circlize)
library(RColorBrewer)
col<-RColorBrewer::brewer.pal(8,"Paired")
circos.initialize(factors=df$Chr,x=df$X)
circos.trackPlotRegion(factors=df$Chr,y=df$Y,
panel.fun=function(x,y){
circos.axis()
},track.height = 0.05)
for(i in 1:8){
highlight.sector(sector.index = paste0("LG",i),col=col[i])
circos.text(CELL_META$xcenter, CELL_META$ycenter,
labels = paste0("LG",i),
sector.index = paste0("LG",i),cex=0.5)
}
circos.trackPlotRegion(factors=df$Chr,y=df$Y)
circos.trackLines(df$Chr,df$X,df$Y,col=col)
circos.trackPlotRegion(factors=df$Chr,y=df$Y)
circos.trackPoints(df$Chr,df$X,df$Y,col=col,cex=0.5)
circos.trackPlotRegion(factors=df$Chr,y=df$Y)
circos.trackHist(df$Chr,df$X,col=col)
circos.clear()
結(jié)果
image.png
還想實(shí)現(xiàn)的其他功能：
想開(kāi)一個(gè)口子用倆給每一圈添加文字標(biāo)簽，如何用代碼實(shí)現(xiàn)還不太清楚录择，想到一個(gè)解決辦法是多加一條染色體，然后出圖后手工編輯將這條染色體刪除掉，然后添加文字標(biāo)簽
如何填加一圈柱形圖呢践樱？
vcftools計(jì)算SNP密度的時(shí)候是否可以設(shè)置滑動(dòng)窗口？
歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本
公眾號(hào)二維碼.jpg
