基于Vcf文件进行基因单倍型分析
2021-04-10 本文已影响0人
DumplingLucky
单倍型，是单倍体基因型的简称，在遗传学上是指在同一染色体上进行共同遗传的多个基因座上等位基因的组合；通俗的说法就是若干个决定同一性状的紧密连锁的基因构成的基因型。
单倍型分析举例：
不同样本中TtBtr-A基因的单倍型分布
TaHY5-like基因的单倍型全球分布
我们以目的基因X来进行实战分析
(1) 提取基因X的基因区SNP
bcftools filter All.vcf.gz --regions chr:start-end > GeneX.vcf
#chr：基因X所在的染色体；start：基因X起始位点；end：基因X终止位点。可以检索基因组gff文件获得
(2) 重新编码vcf文件的基因型
我们这里主要针对双等位基因vcf文件进行分析，vcf中包含./.，0/0，0/1，1/1四种基因型。
编码对应表（转换后 - 转换前）：
- -1 - ./.
- 0 - 0/0
- 1 - 0/1
- 2 - 1/1
使用vcftools可以实现基因型的转换：
vcftools --vcf test.vcf --012
Java脚本也可以实现替换基因型
public class HaploType {
public static void main(String[] args) throws IOException {
this.changeVcfCode("input.txt","output,txt");
}
public static void changeVcfCode(String input, String output) throws IOException {
BufferedReader rd = new BufferedReader(new FileReader(new File(input)));
BufferedWriter wt = new BufferedWriter(new FileWriter(new File(output)));
String line;
List<String> listTab;
while ((line = rd.readLine()) != null ) {
if (line.startsWith("#")) continue;
listTab = Arrays.asList(line.split("\t"));
StringBuilder bf = new StringBuilder();
for (int i = 9; i < listTab.size(); i++) {
if(listTab.get(i).startsWith("1/1")){
bf.append("2\t");
}else if(listTab.get(i).startsWith("0/0")){
bf.append("0\t");
}else if(listTab.get(i).startsWith("0/1") | listTab.get(i).startsWith("1/0") ){
bf.append("1\t");
}else if(listTab.get(i).startsWith("./.")){
bf.append("-1\t");
}else{
bf.append(listTab.get(i)+"\t");
}
}
wt.write(bf.toString()+"\n");
}
rd.close();
wt.close();
}
}
生成haplo.txt文件，没有表头。
haplo.txt
(3) 基因单倍型图谱
- 输入文件
（1）haplo.txt：行代表snp位点，列代表样本。
（2）Annotation_sample.txt：可选的，用于注释单倍型图谱的样本信息，行代表样本，tab隔开。
#Annotation_sample.txt
Growing_Habit Region
Winter SCA
Spring SCA
Winter EA
- 画图
#加载所需的软件包
library(pheatmap)
require(reshape)
require (rworldmap)
require(rworldxtra)
#加载注释内容，根据自己的注释文件进行修改,145是样本数目
annotation_col <- read.table("Annotation.txt",header=T,stringsAsFactors = T)
rownames(annotation_col) = c(1:145)
labels_col = rep(c(""), 145)
ann_colors = list(
Growing_Habit = c(Facultative = "yellow", Spring="red",Winter="blue"),
Region = c(AM = "#228B22",EA = "#D95F02",SCA="#1E90FF",WA="#E7298A", EU="#000000", AF="#FFD700",EAF="#8B008B")
)
#加载haplotype文件
GeneX <- read.table("haplo.txt", header=F, stringsAsFactors = F)
colnames(GeneX) <- c(1:145)
#画图
out <- pheatmap(GeneX, show_colnames=FALSE, legend_breaks = -1:2, legend_labels = c("./.", "0/0", "0/1", "1/1"), cutree_cols=3, cluster_row = FALSE, annotation_col = annotation_col, annotation_colors = ann_colors, main="")
out
基因单倍型分布.png
(4) 基因单倍型地理分布图谱
- 输入文件
（1）haplo.txt：行代表样本，列代表snp位点。
（2）Annotation_sample_location.txt：行代表样本，tab隔开。
#Annotation_sample_location.txt
Lat Lon
36.65 61.15
33.8 62.2
34.36 111.61
- 画图
#加载所需的软件包
library(pheatmap)
require(reshape)
require (rworldmap)
require(rworldxtra)
#加载样本经纬度信息
latlon <- read.table("latlon.txt",header=T,stringsAsFactors = F, fill=TRUE)
rownames(latlon) <- c(1:145)
#对样本经纬度进行聚类
library(cluster)
library(factoextra)
data2 <- latlon[!is.na(latlon$Lat),]
df = scale(data2)
#先求样本之间两两相似性
result <- dist(df, method = "euclidean")
#产生层次结构
result_hc <- hclust(d = result, method = "ward.D2")
#确定样本地理位置聚类的数目，这里分成了40个地理group
data2$type <- cutree(result_hc, k=40)
lat_mean <- tapply(data2[,1],data2$type,mean,na.rm = TRUE)
lon_mean <- tapply(data2[,2],data2$type,mean,na.rm = TRUE)
data2$cluster1 <- NA
data2$cluster2 <- NA
#145是样本数，40是聚类数目
for(i in 1:145) {
for(j in 1:40){
if(data2[i,3] == j ){
data2[i,4] <- as.numeric(lat_mean[j])
data2[i,5] <- as.numeric(lon_mean[j])
}
}
}
new_latlon <- data2[,c(4,5)]
rownames(new_latlon) <- c(1:145)
#根据上述单倍型图谱确定每个地理group的样本以及标注样本对应的单倍型(1，2，3...)
names <- as.numeric(colnames(GeneX[,out$tree_col[["order"]]]))
dataset <- new_latlon[names,]
new_latlon$id <- rownames(new_latlon)
dataset$sub <- NA
colnames(dataset)[1:3] <- c("Lat", "Lon", "id",)
#下面几行代码是确定不同的单倍型cluster的分界线。上面单倍型图谱中分了三种不同的单倍型，从图中确定界限的样本编号(41, 63, 123)，根据下面代码的返回值为样本确定单倍型group(1,2,3)
which(dataset$id==41,)
which(dataset$id==63,)
which(dataset$id==123,)
#确定样本的单倍型类型
dataset$sub[1] <- 1
dataset$sub[2:11] <- 2
dataset$sub[12:145] <- 3
colnames(dataset)[4] <- "Sub.population.3"
#画图
dataset$value <- 1
reshape <- cast(dataset,Lat+Lon~Sub.population.3)
reshape2 <- as.data.frame(reshape)
mapPies(reshape2,xlim=c(-120,140),ylim=c(35,40),nameX="Lon",nameY="Lat",nameZs=c('1','2','3'),symbolSize=1, zColours=c('green','blue','orange'),barOrient='vert',oceanCol="#D1EEEE",landCol="#FFDAB9",main="29")
样本单倍型地理分布.png
参考文献：
【1】https://www.nature.com/articles/s41588-020-00722-w
【2】https://www.nature.com/articles/s41467-020-18738-5
