step2：计算变异位点类型
首先进行初始化，设置参数确定连续变异位点的评判标准，然后读取参考基因组，新建内存空间用于添加候选标记位点。
# 根据GWAS得到的显著性SNP，将其标注到参考基因组上对应位置（300kb内为连续，大于300kb为单标记） ### 数据导入与参数设置 ==================================================================== windows_near <- 300*1000 #默认300kb内为连续 # 连续位点的评判标准是任意两个相邻位点的距离小于windows_near id_list <- list.files("./out/1_GWAS_Result_txt2csv/") Ref <- read_csv("./data/Ref.csv",show_col_types = FALSE) Ref_chr <- Ref all_single <- list() #汇总单标记 all_near <- list() #汇总连续标记
对于每个文件都要进行循环操作，以下列举其中一个文件的计算思路：
file_name <- paste0("./out/1_GWAS_Result_txt2csv/",id) # loc <- substr(file_name,14,18) # job <- substr(file_name,10,26) atom <- str_split(id,"[.]") print(file_name) data <- read_csv(file_name,show_col_types = FALSE) loc <- data$loc[1] job <- paste0(atom[[1]][1],"_",atom[[1]][2],"_",atom[[1]][3]) way <- atom[[1]][3]
以下进行单标记位点的筛选：主要的思路是对每一个位点，计算其前后相邻位点的距离，针对开头和结尾进行特殊考虑，中间的位点识别其是否相邻距离小于阈值。
# 计算基因位置间距 data$longH <- NA data$longQ <- NA data$class <- NA for (i in 2:nrow(data)){ a <- data$POS[i] i <- i+1 b <- data$POS[i] i <- i-2 c <- data$POS[i] i <- i+1 if(i == nrow(data)){ data$class[i] <- "wei" break } if(a-c<0){ data$class[i] <- "shou" next } if(b-a<0){ data$class[i] <- "wei" next } data$longH[i] <- (b-a) data$longQ[i] <- (a-c) } data$class[1] <- "shou" # 对距离进行区分,按照windows_near为区分阈值 for (i in 1:nrow(data)){ if (is.na(data$longH[i]) | is.na(data$longQ[i])){ next } if (data$longH[i]>windows_near & data$longQ[i]>windows_near){ data$class[i] <- "single" } }
找出了连续的标记位点后，需要将其进行整理，转换为数据框格式，第一列为物理位置，第二列为信息，第三列为染色体，第四列为表型区，方法如下：
# 单标记位置信息写入single single <- data.frame(matrix(ncol = 4)) single <- single[-1,] colnames(single) <- c("positon","info","chr","loc") for (i in 1:nrow(data)){ tem_class <- data$class[i] tem_add <- c(data$POS[i],data$ws[i],data$chr[i],data$loc[i]) ifelse(tem_class == 'single',single <- rbind(single,tem_add),"1") ifelse(tem_class == 'shou',single <- rbind(single,tem_add),"2") ifelse(tem_class == 'wei',single <- rbind(single,tem_add),"3") } colnames(single) <- c("positon","info","chr","loc")
接下来进行连续型标记位点的查找筛选，首先初始化一个矩阵用于储存中间变量。
需要注意的是连续型位点有两个位置参数，分别代表这个连续区段的开始和结束，另外需要统计该区段内有多少的实际位点。
near <- data.frame(matrix(ncol = 5)) #初始化矩阵 near <- near[-1,] colnames(near) <- c("p1","p2","info","chr","number")
然后按照染色体进行迭代，每个染色体分开计算，主要步骤如下：
for (x in c(1:7)){ for (y in c("A","B","D")){ chr_id <- paste0(x, y, sep="") # 将x和y连接起来，形成染色体编号chr_id foot <- 0 # 步长，用于迭代计算阅读框，初始化为0 for (i in which(data$chr == chr_id)) { # 遍历data中chr等于chr_id的索引（位点） if (sum(data$chr == chr_id) < 2) { next # 若某个染色体的位点数小于3则跳过当前循环 } if (foot > 0) { foot <- foot - 1 next # 如果foot变量大于0，说明两个位点存在跨越关系，将foot减1，然后跳过当前循环 } n_pos_1 <- data$POS[i] # 获取当前位点的位置信息 for (m in which(data$chr == chr_id)) { n_pos_2 <- data$POS[m] # 获取另一个位点的位置信息 if (n_pos_2 - n_pos_1 < 0) { next # 后一个位点位置小于前一个位点时，跳过当前循环 } else { if (n_pos_2 - n_pos_1 == 0) { next # 两个位点位置相等时，跳过当前循环 } else { if (n_pos_2 - n_pos_1 < windows_near) { foot <- foot + 1 # 任意两个位点距离小于预设窗口大小，增加foot的值，跨越一个位点 if (is.na(data$class[i])) { data$class[i] <- "near" # 如果位点尚不属于特定类别，则将其标记为"near" } } else { if (foot > 10) { n_wn <- paste0(data$SNP[i], "-", data$SNP[m-1], ", [", foot, "], Find in ", job) } else { n_wn <- paste0(data$SNP[i], "-", data$SNP[m-1], ", Find in ", job) } # 构建一个文本信息n_wn，描述了位点之间的关系和发现位置 n_add <- c(n_pos_1, data$POS[m-1], n_wn, data$chr[i], foot) # 将相关信息存储在n_add向量中 if (is.na(data$class[i])) { data$class[i] <- "near" } # 如果位点尚不属于特定类别，则将其标记为"near" break } } } } if (length(data$POS[m-1]) > 0) { if (n_pos_1 != data$POS[m-1]) { near <- rbind(near, n_add) n_add <- c() } } } }} colnames(near) <- c("p1","p2","info","chr","number")
主要是对位点数据进行迭代和处理，根据位点位置以及预设的窗口大小，识别和处理相邻的位点之间的关系，最终将处理后的信息存储在不同变量中。
删除重复数据
由于部分信息统计时出现重复折叠，因此通过下述代码对结果中重复信息进行过滤。
# 删除重复行 near_new <- data.frame(matrix(ncol =5)) near_new <- near_new[-1,] for (i in 1:nrow(near)){ if (!identical(near$p1[i],near$p2[i])){ new_add <- c(near$p1[i],near$p2[i],near$info[i],near$chr[i],near$number[i]) near_new <- rbind(near_new,new_add) } } colnames(near_new) <- c("p1","p2","info","chr","number") near <- near_new
按染色体进行标注
连续型变异区段的标注方法：
for (i in 1:nrow(near)){ if (identical(near$chr[i],"Chr")){ my_a <- which.min(abs(as.numeric(near$p1[i]) - as.numeric(chr_1A$Start.1))) my_b <- which.min(abs(as.numeric(near$p2[i]) - as.numeric(chr_1A$Start.1))) chr_1A$out[my_a:my_b] <- near$info[i] } }
单独型变异位点标注方法：
for (i in 1:nrow(single)){ if (identical(single$chr[i],"1A")){ index <- which.min(abs(as.numeric(single$positon[i]) - as.numeric(chr_1A$X3G.Start.1))) chr_1A$out[index] <- single$info[i] OK <- OK+1 } }
step3：比对参考基因组
第二步中已经对标记位点进行分类整理，接下来将其填入参考基因组对应区间，看看有哪些基因或QTL，主要流程如下：
Ref_file <- "data/Ref.csv" Ref <- read_csv(Ref_file) id_list <- list.files("./out/3_Gene_Maping_Result/") index <- 8 for (id in id_list){ now <- read.csv(paste0("./out/3_Gene_Maping_Result/",id)) atom <- str_split(id,"[.]") job <- atom[[1]][1] tem <- bind_cols(tem,now[,8]) colnames(tem)[index] <- job index <- index + 1 print(id) } for (i in 1:nrow(tem)){ # 判断每个基因有多少个对应的位点 tem$End[i] <- ncol(tem) - sum(tem[i,8:ncol(tem)] == "") - 7 } colnames(tem)[4] <- "CountSNP"
如果想将多个表型或者环境的结果合并到一起，可以通过下面的代码进行操作：
tem$all <- NA # 将多个环境区的数据合并到同一列----------- for (i in 1:nrow(tem)){ var_add <- c() for (m in 8:(ncol(tem)-1)){ if (tem[i,m] == ""){ next }else{ var_add <- c(var_add,tem[i,m]) } } add_info <- str_c(var_add,sep = "",collapse = " ; ") tem$all[i] <- add_info } all_col <- tem[,ncol(tem)] out <- cbind(tem[,1:7],all_col,tem[,8:(ncol(tem)-1)])
由于上述过程涉及到两次迭代循环，计算量比较大，因此对其进行重构优化，实现并行计算，方法如下：
# 多线程重构并行运算----------- library(foreach) library(doParallel) library(stringr) # 获取计算机可用线程数 num_cores <- parallel::detectCores() - 2 cl <- makeCluster(num_cores) registerDoParallel(cl) n <- nrow(tem) result <- foreach(i = 1:n, .combine = rbind) %dopar% { var_add <- c() for (m in 8:(ncol(tem)-1)){ if (tem[i,m] == ""){ next } else { var_add <- c(var_add, tem[i,m]) } } add_info <- stringr::str_c(var_add, sep = "", collapse = " ; ") tem$all[i] <- add_info return(tem[i, ]) } # 关闭并行计算 stopCluster(cl) registerDoSEQ() # 将结果写回原始数据框tem tem <- result
step4：特殊结构变异区段查找
有一些大范围的结果变异，比如CNV等也比较重要，因此通过以下算法实现CNV类型初步判断与整理：
opt_snp_near <- 300*1000 # 相邻片段的最大间距,默认为300kb opt_windows <- 2*1000*1000 # CNV片段的最小宽度，默认阈值为2MB，若连续位点跨度大于该值表示CNV # CNV的评判标准：在连续型位点中，第一个位点和最后一个位点的间距跨度大于opt_windows CNV_find <- as.data.frame(matrix(ncol = 7,nrow = 0)) colnames(CNV_find) <- c("p1","p2","info","chr","number","gap","type") for (page in 1:length(SNP_near_all)){ job <- names(SNP_near_all)[page] print(job) job_name <- str_replace(job,".csv","") if (nrow(SNP_near_all[[page]]) < 3){next} # 过滤小于2行的数据集 df <- SNP_near_all[[page]] df$gap <- as.numeric(df$p2) -as.numeric(df$p1) df$type <- NA for (i in 1:nrow(df)){ # 判断三种特殊情况，过滤具有单个位点的染色体 if (i == 1 & df$chr[i] != df$chr[i+1]){df$type[i] <- "del";next} if (i == 1 & df$chr[i] == df$chr[i+1]){ if (as.numeric(df$p1[i+1]) - as.numeric(df$p2[i]) < opt_snp_near){ df$type[i] <- str_c(job_name,"_CNV_1") next }else{next} } if (i == nrow(df) & df$chr[i] != df$chr[i-1]){df$type[i] <- "del";next} if (df$chr[i] != df$chr[i-1] & df$chr[i] != df$chr[i+1]){df$type[i] <- "del";next} } df <- df[which(is.na(df$type)),] # 删除单个染色体 if (nrow(df) == 0){next} # 若删除后数据为空则跳过 # 进行CNV判断 temp_type <- 1 # CNV片段分类计数 for (i in 1:nrow(df)){ foot <- 0 # 记录前进步数 if (!is.na(df$type[i])){ next } is_add <- F # 记录是否添加CNV # 最后一行特殊处理 if (i == nrow(df)){ if (as.numeric(df$p1[i]) - as.numeric(df$p2[i-1]) < opt_snp_near){ df$type[i] <- df$type[i-1] }else{next} } for (index in (i+1):nrow(df)) { #从下一行开始查找连续区段，直到终止 # if (foot != 0){foot <- foot -1 ;next} if (as.numeric(df$p1[index]) - as.numeric(df$p2[i+foot]) < opt_snp_near){ df$type[i+foot] <- str_c(job_name,"_CNV_",temp_type) df$type[index] <- str_c(job_name,"_CNV_",temp_type) is_add <- T foot <- foot + 1 }else{ break } } if (is_add){ temp_type <- temp_type +1 } } out <- df[which(!is.na(df$type)),] CNV_find <- rbind(CNV_find,out) }
行文至此，流程结束！以上内容中可能还有一些不完善的地方，欢迎加我好友交流，感谢你能看到这里，如果愿意的话点个在看或转发支持一下吧。
END
© 素材来源于网络，侵权请联系后台删除
笔记合集，点击直达 
