本篇文章給大家分享的是有關(guān)怎樣使用R語(yǔ)言利用vcf格式文件計(jì)算核苷酸多樣性，小編覺(jué)得挺實(shí)用的，因此分享給大家學(xué)習(xí)，希望大家閱讀完這篇文章后可以有所收獲，話不多說(shuō)，跟著小編一起來(lái)看看吧。
創(chuàng)新互聯(lián)公司成立以來(lái)不斷整合自身及行業(yè)資源、不斷突破觀念以使企業(yè)策略得到完善和成熟，建立了一套“以技術(shù)為基點(diǎn)，以客戶需求中心、市場(chǎng)為導(dǎo)向”的快速反應(yīng)體系。對(duì)公司的主營(yíng)項(xiàng)目，如中高端企業(yè)網(wǎng)站企劃 / 設(shè)計(jì)、行業(yè) / 企業(yè)門戶設(shè)計(jì)推廣、行業(yè)門戶平臺(tái)運(yùn)營(yíng)、app軟件開發(fā)、手機(jī)網(wǎng)站制作設(shè)計(jì)、微信網(wǎng)站制作、軟件開發(fā)、四川雅安服務(wù)器托管等實(shí)行標(biāo)準(zhǔn)化操作，讓客戶可以直觀的預(yù)知到從創(chuàng)新互聯(lián)公司可以獲得的服務(wù)效果。
重慶分公司
重慶分公司
bcftools view snp.vcf.gz scaffold_1 > popgenome-vcf/scaffold_1
bcftools view snp.vcf.gz scaffold_2 > popgenome-vcf/scaffold_2
如果當(dāng)前目錄下只有vcf格式文件，會(huì)遇到報(bào)錯(cuò)Failed to open .vcf.gz: could not load index，可以參考 https://www.cnblogs.com/chenwenyan/p/11945445.html
tabix -p vcf snp.vcf.gz
如果當(dāng)前目錄下沒(méi)有popgenome-vcf這個(gè)目錄，還需要新建目錄
mkdir popgenome-vcf
今天參考的文章里寫道 In theory, the r PopGenome can read VCF files directly, using the readVCF function. However, because our samples are haploid, we need to use a different function, r readData, which requires a folder with a separate VCF for each scaffold. 這個(gè)是為什么呢？
#install.packages("PopGenome")
library(PopGenome)
getwd()
setwd("VCF/")
snp<-readData("popgenome-vcf",format = "VCF")
get.sum.data(snp)
這里可以直接統(tǒng)計(jì) 轉(zhuǎn)換和顛換的比例
pops<-get.individuals(snp)[[1]]
pop1<-pops[grep("B\\.bam",pops)]
pop2<-pops[grep("b\\.bam",pops)]
pop1
pop2
snp<-set.populations(snp,list(pop1,pop2))
snp@populations
snp<-F_ST.stats(snp)
get.F_ST(snp)
這里的指標(biāo)都是什么意思呢？
get.diversity(snp)[[1]]
這里的指標(biāo)也看不懂是什么意思呀
win_snp<-sliding.window.transform(snp,
width = 10000,
jump = 2000,type = 2)
win_snp<-F_ST.stats(win_snp)
win_snp@nucleotide.F_ST
win_snp@nuc.diversity.within
library(ggplot2)
win_fst <- data.frame(x=1:dim(win_snp@nucleotide.F_ST)[1],
y=win_snp@nucleotide.F_ST[,1])
head(win_fst)
p1<-ggplot(win_fst,aes(x=x,y=y))+
geom_point()+
geom_line()+
theme_bw()+
theme(panel.grid = element_blank())+
scale_x_continuous(breaks = win_fst$x,
labels = win_fst$x)+
labs(x=NULL,y=NULL,title = "FST")
ggsave("FST.pdf",p1,width = 15,height = 4)
bb_div <- win_snp@nuc.diversity.within[,1] # diversity among B (bb = "big B")
lb_div <- win_snp@nuc.diversity.within[,2] # diversity among B (lb = "little b")
bb_div
df1<-data.frame(x=1:length(bb_div),y=bb_div)
df2<-data.frame(x=1:length(lb_div),y=lb_div)
p2<-ggplot()+
geom_line(data=df1,aes(x=x,y=y),color="red")+
geom_point(data=df1,aes(x=x,y=y),size=2,color="red")+
geom_line(data=df2,aes(x=x,y=y),color="blue")+
geom_point(data=df2,aes(x=x,y=y),size=2,color="blue")+
theme_bw()+labs(x=NULL,y=NULL)
ggsave("diversity.pdf",p2,width = 15,height = 4)
以上就是怎樣使用R語(yǔ)言利用vcf格式文件計(jì)算核苷酸多樣性，小編相信有部分知識(shí)點(diǎn)可能是我們?nèi)粘９ぷ鲿?huì)見(jiàn)到或用到的。希望你能通過(guò)這篇文章學(xué)到更多知識(shí)。更多詳情敬請(qǐng)關(guān)注創(chuàng)新互聯(lián)行業(yè)資訊頻道。
