用popart構建常染色體單倍型網絡（Autosomal haplotypes network construction with popart）

1）將vcf轉化為plink格式，假定輸入的vcf文件名為：17893893-17898893.vcf,也可以參考鏈接：將vcf文件轉化為plink格式並且保持phasing狀態

/vcftools --vcf 17893893-17898893.vcf --plink-tped --out 17893893-17898893 

　　

/plink --tfile 17893893-17898893 --recode --out 17893893-17898893

　　

2) 用PLINK確定要研究的位點是否處於連鎖的狀態；生成blocks和blocks.det兩種后綴格式文件；

/plink --file 17893893-17898893 --blocks no-pheno-req --out 17893893-17898893

　　

 

 以上結果說明rs62033246|rs7189274|rs7187782|rs62033247|rs7194607|rs7200159|rs199711091|rs2361776|rs8054992|rs1845376|rs35902958|rs9928657|rs9928743|rs56078865|rs62033249|rs62033250|rs9936044|rs7184460|rs4781927|rs7202082|rs4780669|rs7202439|rs7195939|rs57418698|rs34615631|rs533867711|rs555603620|rs567435894|rs7499772|rs62033251|rs56248612|rs7202488|rs61671028|rs62033253|rs9928974這35個位點是連鎖的

 

3) 提取這35個位點作為接下來的單倍型網絡構建；去掉vcf的頭文件，保存為txt格式，見如下圖，17893893-17898893.txt：

4）准備17893893-17898893_singstring.txt文件，該文件其實就是去掉 0|0，0|1，1|0，1|1的“|”。

 

5）接下來，生成兩條單鏈，用R輸入如下命令：

install.packages("xlsx")
library(xlsx)
kg_ame<-read.table("E:/17893893-17898893.txt")
kg_ame<-data.frame(kg_ame);
library("stringr");
newdata=c()
for (i in 0:34){
  i=i+1
  kk=kg_ame[i,6:2509];
  kk[which(kk=="0|0")]=paste(kg_ame[i,4],kg_ame[i,4]);
  kk[which(kk=="0|1")]=paste(kg_ame[i,4],kg_ame[i,5]);
  kk[which(kk=="1|0")]=paste(kg_ame[i,5],kg_ame[i,4]);
  kk[which(kk=="1|1")]=paste(kg_ame[i,5],kg_ame[i,5]);
  newdata=rbind(newdata,kk)
}
##6：2509指的是2000多個樣本對應的鹼基；
##0：34指的是35個鹼基；這一步是將0|0，0|1，1|0，1|1轉化為鹼基形式，類似vcf轉化為Ped的步驟；
openxlsx::write.xlsx(newdata,file ="E:/17893893-17898893_change.xlsx")
kk=kg_ame[1,7:15];kk
kk[which(kk=="0|0")]=paste(kg_ame[1,4],kg_ame[1,4]);kk
#####分為兩條鏈，轉換為ped格式；
library(xlsx)
kg_ame<-read.table("E:/17893893-17898893_singstring.txt")
kg_ame<-data.frame(kg_ame);
library("stringr");
newdata=c()
for (i in 0:34){
  i=i+1
  kk=kg_ame[i,6:5013];
  kk[which(kk=="0")]=kg_ame[i,4];
  kk[which(kk=="1")]=kg_ame[i,5];
  newdata=rbind(newdata,kk)
}
openxlsx::write.xlsx(newdata,file ="E:/17893893-17898893_changesingstring.xlsx")
se1<-seq(from=1,to=5008,by=2);
se2<-seq(from=2,to=5008,by=2);
stringone<-newdata[,se1];
stringtwo<-newdata[,se2];
openxlsx::write.xlsx(stringone,file ="E:/17893893-17898893_changesingstringone.xlsx")
openxlsx::write.xlsx(stringtwo,file ="E:/17893893-17898893_changesingstringtwo.xlsx")
####以上步驟為轉為兩條單鏈

　　

6) 將上述的兩條單鏈生成fas格式；

 

paste -d "\n"  17893893-17898893_changesingstringone_firstdot.txt 17893893-17898893_changesingstringone.txt > 17893893-17898893_changesingstringone_dot.fas
paste -d "\n"  17893893-17898893_changesingstringtwo_firstdot.txt 17893893-17898893_changesingstringtwo.txt > 17893893-17898893_changesingstringtwo_dot.fas
cat 17893893-17898893_changesingstringone_dot.fas 17893893-17898893_changesingstringtwo_dot.fas >> 17893893-17898893_changesingstring_onetwo_dot.fas #這兩條鏈合並在一起

17893893-17898893_changesingstringone_firstdot.txt 文件格式如下：

17893893-17898893_changesingstringone.txt的文件格式如下：這個文件就是5）生成出來的其中一條單鏈。

 

生成的fas文件如下：

第二條鏈的處理同上。最后一步，就是把這兩條鏈合並在一起

7）用DnaSP6將fas文件生成nex文件

 

8）統計nex文件中不同群體的haplotype的數量

准備17893893-17898893_hap.xlsx文件

准備allsample.txt文件

輸入如下命令：

 

#####統計nex生成的所有haplotype所在的樣本和群體，這里假定生成了100個haplotype
library(dplyr)
library(xlsx)
model_57hanchip<-read.xlsx("E:/17893893-17898893_hap.xlsx",2)
model_57hanchip<-data.frame(model_57hanchip)
kg_ame<-read.table("E:/allsample.txt")
kg_ame<-data.frame(kg_ame);
newdata=c()
for (i in 0:99){
  i=i+1
  HAP1<-kg_ame[,2];
  df1 <- data.frame(id = HAP1, kg_ame[,3],kg_ame[,4])
  HAP2<-model_57hanchip[,i]
  df2 <- data.frame(id = HAP2)
  k=merge(df1, df2, by = "id",all = FALSE)
  j=table(k[,3])
  newdata=rbind(newdata,j)
  print(i)
}
openxlsx::write.xlsx(newdata,file ="E:/17893893-17898893_hapcount.xlsx")

　　

　

9）修改生成的nex文件

DnaSP6生成的nex文件並沒有BEGIN TRAITS的頭文件，因此這一步是需要手動修改的。

這一步，除了圈出來的數字和AFR AMR EAS EUR SAS是需要修改，其他照抄。

 10）用popart構建單體型網絡

導入文件后，選擇適合的算法，即可生成。