一、准备工作
需要安装devtools、vcfR、binmapr
install.packages("devtools")
install.packages("vcfR")
#加载devtools时提示需要usethis包
install.packages("usethis")
library(usethis)
library(devtools)
library(vcfR)
devtools::install_github("xuzhougeng/binmapr")
- 安装devtools时出现问题，详见笔记：
R-Rstudio-server无法安装devtools解决方法 - 简书 (jianshu.com) - 安装binmapr时出现问题，详见笔记：
R-devtools无法安装binmapr包 - 简书 (jianshu.com)
二、读取数据
#设置工作路径
setwd("~/workspace/BSA/practice/")
#vcfR读取VCF文件
vcf <- read.vcfR("4.variants_filter/snps.vcf")
#从VCF对象中提取两个关键信息，AD(Allele Depth)和GT(Genotype)
gt <- extract.gt(vcf)
ad <- extract.gt(vcf, "AD")
gt:基因型矩阵，过滤后仅剩"0/0", "0/1","1/1"
ad:等位基因的count数.
> head(gt)
SRR6327815 SRR6327816 SRR6327817 SRR6327818
chr01_1151 "0/0" "1/1" "0/0" "0/1"
chr01_6918 "0/0" "1/1" "0/0" "0/1"
chr01_17263 "0/0" "1/1" "0/0" "0/1"
chr01_21546 "1/1" "1/1" "0/0" "0/1"
chr01_24732 "1/1" "0/0" "0/0" "0/0"
chr01_33667 "1/1" "1/1" "0/0" "0/1"
> head(ad)
SRR6327815 SRR6327816 SRR6327817 SRR6327818
chr01_1151 "14,0" "0,18" "10,0" "11,3"
chr01_6918 "31,0" "0,34" "32,0" "21,5"
chr01_17263 "46,0" "0,37" "20,0" "21,9"
chr01_21546 "0,26" "0,20" "17,0" "16,3"
chr01_24732 "0,30" "36,0" "22,0" "31,0"
chr01_33667 "0,28" "0,34" "28,0" "29,12"
SRR6327815(KY131), SRR6327816(DN422),
SRR6327817(T-pool), SRR6327818(S-pool)
需要测试是选择KY131是"0/0", "DN422" 是 "1/1"的点
还是选择KY131是 "1/1"，"DN422" 是"0/0"的点
#xxx15对应KY131，xxx16对应DN422
#选取gt中KY0/0,DN1/1
mask <- which(gt[,"SRR6327815"] == "0/0" & gt[,"SRR6327816"] == "1/1")
length(mask)#结果显示[1] 140788
#选取gt中KY1/1,DN0/0
verse_mask <- which(gt[,"SRR6327815"] == "1/1" & gt[,"SRR6327816"] == "0/0")
length(mask)#结果显示[1] 221712
#先测试0/0,1/1的点
ad_flt <- ad[mask,c("SRR6327817", "SRR6327818")]
colnames(ad_flt) <- c("T_Pool", "S_Pool")
接下来根据AD来计算SNP-index简单理解SNP-index就是与亲本基因型不同的计数，即
SNP-index=AD（ALT）/AD(REF)+AD(ALT)
AD:allele depth,指样本中每一种等位基因的reads覆盖度。vcf文件中，一般对于二倍体或多倍体，会用逗号进行分隔，前者是REF，后者是ALT
> freq <- calcFreqFromAd(ad_flt, min.depth = 10, max.depth = 100)
#弹出来了这个报错，探索了半天，感觉可能是因为ad_flt是matrix，不支持$选择特定列？
$ operator is invalid for atomic vectors
#输入函数，看看是什么
> calcFreqFromAd
function (x, min.depth = 10, max.depth = 200)
{
Ref_COUNT <- x$refmat #逗号前的数值
Alt_COUNT <- x$altmat #逗号后的数值
AD_Depth <- Ref_COUNT + Alt_COUNT #Allele Depth
NA_pos <- which(AD_Depth < min.depth | AD_Depth > max.depth) #超过最小和最大深度的值
Ref_COUNT[NA_pos] <- NA
Alt_COUNT[NA_pos] <- NA #将Ref和Altcount中这些超出阈值的值设置为NA
AD_freq <- round(Alt_COUNT/(Ref_COUNT + Alt_COUNT), 2) #计算SNP-index
return(AD_freq)
}
<bytecode: 0x557d38229660>
<environment: namespace:binmapr>
然后我发现，vcfR包里是有相关函数(masplit)可以用来提取字符的
ref_count <- masplit(ad_flt, record = 1, sort = 0)
alt_count <- masplit(ad_flt, record = 2, sort = 0)
ad_depth <- ref_count + alt_count
NA_pos <- which(ad_depth < 10 | ad_depth >100)
ref_count[NA_pos] <- NA
alt_count[NA_pos] <- NA
freq <- round(alt_count/(ref_count + alt_count), 2)
#参照函数的步骤，成功得到了freq
freq2 <- freq[Matrix::rowSums(is.na(freq)) == 0, ]
接下来进行画图
par(mfrow = c(3,4))
for (i in paste0("chr", formatC(1:12, width = 2, flag=0)) ){
#i对应chr01、chr02……chr12
#根据染色体编号，将freq2分为12个矩阵
freq_flt <- freq2[grepl(i,row.names(freq2)), ]
#提取位置信息：substring从第七位（即除去染色体编号+“_”）开始的字符串
pos <- as.numeric(substring(row.names(freq_flt), 7))
#开始绘图，以pos为横坐标，freq_flt第二列和第一列差值为纵坐标，
#y轴界限为-1~1
plot(pos, freq_flt[,2] - freq_flt[,1], ylim = c(-1,1),
#pch：指定描绘点时的符号，cex：符号大小
pch = 20, cex = 0.2,
#x轴标题和y轴标题
xlab = i,
#expression，注明数学符号△
ylab = expression(paste(Delta, " " ,"SNP index")))
}
KY1/1DN0/0中
以上为两种情况下所做的△SNP_index图
文献中指出，候选基因被定位在六号染色体的20-25M处，从图中可知KY0/0，DN1/1这种情况与文献中情况类似，因此在最初的实验设计中应该就是以KY0/0,DN1/1为标准
0：REF，1：ALT
文献中除了用到了SNP-index，还用到了ED（Euclidean distance）方法进行定位
where Aaa, Caa, Gaa and Taa represent the frequency of bases A, C, G and T in the T-pool, respectively. Aab, Cab, Gab and Tab represent the frequency of bases A, C, G and T in the S-pool, respectively.
Aaa、Caa、Gaa、Taa：T-pool的碱基
Aab、Cab、Gab、Tab：S-pool的碱基
The depth of each base in different pools and the ED value of each SNP loci were calculated.
计算不同池中的depth及每一个SNP位点的ED值
To eliminate the background noise, the ED value was powered, and the ED^5 was used as the associated value.
为消除噪音，返回值为ED^5
mask <- which(gt[,"SRR6327815"] == "0/0" & gt[,"SRR6327816"] == "1/1")
ad_flt <- ad[mask,c("SRR6327817", "SRR6327818")]
#apply函数避免for循环，margin = 1：对行进行操作
ED_list <- apply(ad_flt, 1, function(x){
#strsplit()，按逗号分割字符串，返回列表
#unlist()，将列表转化为向量，as.numeric()将数据类型改为numeric
count <- as.numeric(unlist(strsplit(x, ",",fixed = TRUE,useBytes = TRUE)))
#depth1：T-pool；depth2：S-pool
depth1 <- count[1] + count[2]
depth2 <- count[3] + count[4]
#sqrt()：开方
#count[3]/depth2:S-pool
ED <- sqrt((count[3] / depth2 - count[1] / depth1)^2 +
(count[4] / depth2- count[2] /depth1)^2)
return(ED^5)
})
par(mfrow = c(3,4))
for (i in paste0("chr", formatC(1:12, width = 2, flag=0)) ){
ED_flt <- ED_list[grepl(i,names(ED_list))]
pos <- as.numeric(substring(names(ED_flt), 7))
plot(pos, ED_flt,
pch = 20, cex = 0.2,
xlab = i,
ylab = "ED")
}
与文献对比
