亲缘关系分析工具:亲属对分类统计+画树+mtY单倍型

Content

(0) why(1) doc(2) func(3) R(3.0) lib+input(3.1) 亲缘关系对分类统计(3.2) 亲缘关系对分类统计(3.3) 画树pedigree(4) pedtools所属系列的书

(0) why

主播做亲缘关系分析的体感：前面从mapping到什么LD等等用了很多高级复杂的算法，但是到统计亲属对时变为了普通的靠勤劳的双手和人畜不分的眼神去数，画树时则是纯粹的PPT对齐边框手稳大赛。主播感觉很疲惫，很恼怒，所以写了脚本。

(1) doc

沿用pedtools的ped。因为这个包内部有一些check环的逻辑，后面还能画图，比较完善，好用。但是ped中祖先个体的父母0非常恼人所以前面加了一段。但是这个ped还是要手动写一段时间的，或者maybe有什么代劳

(2) func

(2.1) 亲缘关系对分类统计：先画deg1有向网络图，然后从common_ancestry找最短路径，得出两两之间的degNum【本文包饺子的醋】。细分另说。(2.2) mtY理论分组，可以和实际结果check。(2.3) 画树pedigree：内置label可以标识属性，e.g.线条虚实=采样,填充颜色=mt,边框颜色=Y

(3) R

(3.0) lib+input

# local({r <- getOption("repos")# r["CRAN"] <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"# options(repos=r)}# )# options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")# options("download.file.method"="libcurl")# options("url.method"="libcurl")## ###————1 Bioconductor，可以兼容后两个————#### if(!requireNamespace("BiocManager",quietly = TRUE))#   install.packages("BiocManager")# # BiocManager::install("GEOquery")# ###————2 在CRAN 存储库的r包————#### install.packages("xxx")# ###————3 在GitHub储存库中r包————#### library(devtools)# install_github("authorName/repositoryName")# # To install ggplot2 from github:# devtools::install_github("tidyverse/ggplot2")library(tidyverse)library(pedtools)library(igraph)library(GWASTools)library(RColorBrewer)library(data.table)library(ggbreak)# "D:/projects/proAAA"pre_DIR <- "D:/pro/proAAA"wkdir <- paste0(pre_DIR,"/pedDF")setwd(wkdir)getwd()usesite <- 'proAAA'###### Part0: input #########int_pedDF <- read.table(paste0('./','pedDF_',usesite,'.csv'),sep=',',header = T,row.names = NULL)# 本次未采样个体的命名特点是以"N"开头，所以step5用这个notsampled_prefixnotsampled_prefix <- "N"unsampled_ids <- int_pedDF$id[grepl(paste0("^", notsampled_prefix), int_pedDF$id)]

(3.1) 亲缘关系对分类统计

###### Part1: pedDF #########deg1caseDF_2_pedDF <- function(tmp0_df){# 给非0的DEFmom/dad添加上父母信息，用于画pedigree# 已经是0的就说明到此为止# 如果不早点赋0，最后画出来太乱了没法看  use_momids <- setdiff(unlist(unique(tmp0_df$momid)),'0')  use_dadids <- setdiff(unlist(unique(tmp0_df$dadid)),'0')  add0_momids <- setdiff(use_momids,tmp0_df$id)  add0_dadids <- setdiff(use_dadids,tmp0_df$id)if(length(add0_momids)>0){    add0_mom_df <- data.frame(      id=add0_momids,      momid=rep(0,length(add0_momids)),      dadid=rep(0,length(add0_momids)),      sex=rep('F',length(add0_momids))    )    tmp0_df <- tmp0_df %>% rbind(add0_mom_df)  }if(length(add0_dadids)>0){    add0_dad_df <- data.frame(      id=add0_dadids,      momid=rep(0,length(add0_dadids)),      dadid=rep(0,length(add0_dadids)),      sex=rep('M',length(add0_dadids))    )    tmp0_df <- tmp0_df %>% rbind(add0_dad_df)  }  tmp1_df <- tmp0_df %>%    mutate(Sex=ifelse(sex=='F', 2, 1)) %>%    select(-sex) %>%    rename(sex=Sex)  tmp1_ped <- ped(id=tmp1_df$id,fid=tmp1_df$dadid,mid=tmp1_df$momid,sex=tmp1_df$sex)return(tmp1_ped)}pedDF <- deg1caseDF_2_pedDF(tmp0_df=int_pedDF)GWAS_pedDF <- as.data.frame(pedDF)###### __1.1 pcsib_pairdtfm #########pc_pairdtfm <- data.frame(pair1 = character(),                     pair2 = character(),                     relation = character(),                     sex1 = integer(),                     sex2 = integer(),                     stringsAsFactors = FALSE)for (i in1:nrow(GWAS_pedDF)) {# 提取当前行的值  id <- GWAS_pedDF$id[i]  fid <- GWAS_pedDF$fid[i]  mid <- GWAS_pedDF$mid[i]  sex <- GWAS_pedDF$sex[i]  new_row1 <- data.frame(pair1 = id, pair2 = fid, relation = 'pc_21', sex1 = sex, sex2 = 1)  new_row2 <- data.frame(pair1 = id, pair2 = mid, relation = 'pc_21', sex1 = sex, sex2 = 2)  pc_pairdtfm <- rbind(pc_pairdtfm, new_row1, new_row2)}sib_pairdtfm <- data.frame(pair1 = integer(), pair2 = integer(), relation = character(), sex1 = integer(), sex2 = integer())for (i in1:(nrow(GWAS_pedDF)-1)) {# 获取当前行的信息  current_row <- GWAS_pedDF[i, ]  id1 <- current_row$id  sex1 <- current_row$sex  fid1 <- current_row$fid  mid1 <- current_row$midif (!is.na(fid1) && !is.na(mid1) && (fid1 != 0 || mid1 != 0)) {# 遍历当前行之后的所有行for (j in (i + 1):nrow(GWAS_pedDF)) {      next_row <- GWAS_pedDF[j, ]      fid2 <- next_row$fid      mid2 <- next_row$midif (!is.na(fid2) && !is.na(mid2) && fid1 == fid2 && mid1 == mid2) {        id2 <- next_row$id        sex2 <- next_row$sex        sib_pairdtfm <- rbind(sib_pairdtfm, data.frame(pair1 = id1, pair2 = id2, relation = "sib", sex1 = sex1, sex2 = sex2))      }    }  }}pcsib_pairdtfm <- pc_pairdtfm %>%  filter((pair2!=0) &(fid!=0) & (mid!=0)) %>%  rbind(sib_pairdtfm)###### __1.2 kinship-net ####### 亲缘关系在有向图中的表现形式：两者有共同祖先same_upper，即在有向图中至少存在一个点，对这两点都存在强连接（符合有向边本身方向性的连接）# (1)一级关系，pc和sib，这部分实际上只能手动整理后作为输入，表现为两点之间存在1条边直接相连；# (2)二级及以上关系，即利用pc和sib得到有向网络图后，有same_upper的两者连接所需边数即为deg数###### ____1.2.1 构建net  ########## ！！！！！！Notice: pcsib_pairnet不能用0，而且net的fromto在前面pair12列存在时失效# 父母->孩子是单向边# 兄弟姐妹之间是双向边nodes <- union(pcsib_pairdtfm$pair1,pcsib_pairdtfm$pair2)one_way_edges <- pcsib_pairdtfm %>%  filter(relation == "pc_21") %>%  mutate(from = pair2,to = pair1) %>%  select(from,to,relation,sex1,sex2)pcsib_pairnet <- graph_from_data_frame(one_way_edges, directed = TRUE)# 创建双向边sib_edges <- pcsib_pairdtfm %>%  filter(relation == "sib") %>%  mutate(from = pair2,to = pair1) %>%  select(from,to,relation,sex1,sex2)# 正向边：sib_edges_forward <- c(rbind(sib_edges$from, sib_edges$to))# 反向边：sib_edges_backward <- c(rbind(sib_edges$to, sib_edges$from))pcsib_pairnet <- pcsib_pairnet %>%  add_edges(sib_edges_forward) %>%  add_edges(sib_edges_backward)###### ____1.2.2(略慢) 计算亲属对，得到antientDF+same_upper #########antientDF <- data.frame(upper = character(), lower = character(), re_filter = character(), stringsAsFactors = FALSE)vertices <- V(pcsib_pairnet)$name# 生成所有可能的节点对组合，排除自身到自身的组合pairs <- expand.grid(upper = vertices, lower = vertices)pairs <- pairs[pairs$upper != pairs$lower, ]results <- vector("list", nrow(pairs))# 遍历所有节点对for (idx in seq_len(nrow(pairs))) {  i <- pairs$upper[idx]  j <- pairs$lower[idx]# 计算从 i 到 j 的有向路径长度  path <- shortest_paths(pcsib_pairnet, from = i, to = j, mode = "out")$vpath[[1]]if (length(path) >= 2) {  # 路径长度至少为 2    results[[idx]] <- data.frame(upper = i, lower = j, re_filter = "need", stringsAsFactors = FALSE)  }}# 将结果列表合并为一个数据框antientDF <- do.call(rbind, results)here_uppers <- unique(antientDF$upper)# 使用data.table的rbindlist高效合并same_upper_list <- list()for (upper in here_uppers) {  lowers <- antientDF$lower[antientDF$upper == upper]if (length(lowers) > 1) {# 生成所有组合并创建数据框    pair_mat <- combn(lowers, 2)    same_upper_list[[upper]] <- data.table(      lower_pair1 = pair_mat[1, ],      lower_pair2 = pair_mat[2, ],      re_filter = "need"    )  }}same_upper <- rbindlist(same_upper_list)same_upper <- unique(same_upper)###### ____1.2.3 直接赋值deg1，得到pcsibdeg_relaDF #########all_pairs <- expand.grid(V(pcsib_pairnet)$name, V(pcsib_pairnet)$name)all_pairs <- all_pairs[all_pairs$Var1 != all_pairs$Var2, ]# 从 pcsib_pairdtfm 中查找已知的关系pcsib_relaDF <- data.frame(  pair1 = as.character(all_pairs$Var1),  pair2 = as.character(all_pairs$Var2),  stringsAsFactors = FALSE)pcsib_relaDF <- pcsib_relaDF %>%  left_join(pcsib_pairdtfm,by=c('pair1'='pair1','pair2' = 'pair2')) %>%  rename(relatedness_A = relation) %>%  left_join(pcsib_pairdtfm,by=c('pair1'='pair2','pair2' = 'pair1')) %>%  rename(relatedness_B = relation) %>%  mutate(relatedness = ifelse(is.na(relatedness_A), relatedness_B, relatedness_A)) %>%  select(pair1, pair2, relatedness)  # 保留匹配到的关系# 如果未找到关系，则根据图的结构判断# pcsibdeg_relaDF这步实际上只是把n个人与其他(n-1)个人的连接边数全部列举出来，但是存在pair1/2重复且并未过滤same_upperpcsibdeg_relaDF <- pcsib_relaDFsetDT(pcsibdeg_relaDF)# 只处理NA行pcsibdeg_relaDF[is.na(relatedness),                relatedness := {                  path <- shortest_paths(pcsib_pairnet, from = pair1, to = pair2, mode = "all")$vpath[[1]]if(length(path) > 0) paste0("deg", length(path)-1) else"un"                },                by = .I]# need这块过滤得到的结果，才是真正存在kinshipneed_pcsibdeg_relaDF <- pcsibdeg_relaDF %>%  left_join(antientDF,by=c('pair1'='upper','pair2' = 'lower')) %>%  rename(antient_filter_A=re_filter) %>%  left_join(antientDF,by=c('pair2'='upper','pair1' = 'lower')) %>%  rename(antient_filter_B=re_filter) %>%  mutate(antient_filter = ifelse(is.na(antient_filter_A),antient_filter_B, antient_filter_A)) %>%  select(pair1, pair2, relatedness,antient_filter)%>%  left_join(same_upper,by=c('pair1'='lower_pair1','pair2' = 'lower_pair2')) %>%  rename(sameupper_filter_A=re_filter) %>%  left_join(same_upper,by=c('pair1'='lower_pair2','pair2' = 'lower_pair1')) %>%  rename(sameupper_filter_B=re_filter) %>%  mutate(sameupper_filter = ifelse(is.na(sameupper_filter_A),sameupper_filter_B, sameupper_filter_A)) %>%  filter((sameupper_filter=='need')|(antient_filter=='need')) %>%  select(pair1, pair2, relatedness) %>%  filter(pair1<pair2)###### ____1.2.4(可跳过) 借助GWAS细分deg2/3 ########## GWAS包对于deg2/3内部有更详细的区分函数，所以添加了相应列(直接应用这个函数的结果)，如有必要可以用于后续分析# https://www.rdocumentation.org/packages/GWASTools/versions/1.18.0/topics/pedigreePairwiseRelatednessGWAS_tmp <- data.frame(  family=rep('f1',nrow(GWAS_pedDF)),  individ=GWAS_pedDF$id,  father=GWAS_pedDF$fid,  mother=GWAS_pedDF$mid,  sex12=GWAS_pedDF$sex) %>%  mutate(sex = case_when(    sex12 == 1 ~ 'M',TRUE ~ 'F'  )) %>%  select(-sex12)GWAS_rst <- GWASTools::pedigreePairwiseRelatedness(GWAS_tmp)GWAS_PAIRrst <- GWAS_rst$relativeprs###### ____1.2.5 去掉未采样个体，得到sampled_relaDF #########sampled_relaDF <- need_pcsibdeg_relaDF %>%  filter(!(pair1 %in% unsampled_ids)) %>%  filter(!(pair2 %in% unsampled_ids)) %>%  left_join(GWAS_PAIRrst[,c('Individ1','Individ2','relation')],            by=c('pair1'='Individ1','pair2' = 'Individ2')) %>%  rename(GWAS_relation_A=relation)  %>%  left_join(GWAS_PAIRrst[,c('Individ1','Individ2','relation')],            by=c('pair1'='Individ2','pair2' = 'Individ1')) %>%  rename(GWAS_relation_B=relation) %>%  mutate(GWAS_relation = ifelse(is.na(GWAS_relation_A), GWAS_relation_B, GWAS_relation_A)) %>%  select(-c(GWAS_relation_A,GWAS_relation_B)) %>%  mutate(deg = case_when(    relatedness %in% c("pc_21", "sib") ~ 'deg1',TRUE ~ relatedness  )) %>%  mutate(re12_deg3more = case_when(    relatedness == "pc_21" ~ "pc",    relatedness == "sib" ~ "sib",    GWAS_relation == "GpGc" ~ "gr",    GWAS_relation == "HS" ~ "hsib",    GWAS_relation == "Av" ~ "avu",TRUE ~ relatedness  )) %>%  mutate(re1_deg2more = case_when(    relatedness == "pc_21" ~ "pc",TRUE ~ relatedness  )) %>%  rename(re123_deg4more=GWAS_relation) %>%  select(pair1, pair2, deg,re1_deg2more,re12_deg3more,re123_deg4more)write.csv(sampled_relaDF,paste0('./',usesite,'_relaDF_deg_re12_re123',".csv"), row.names = F, fileEncoding = "UTF-8")###### ____1.2.6 添加un关系 #######all_ids <- labels(pedDF)sampled_ids <- setdiff(all_ids, unsampled_ids)# R语言利用tidyverse，已有数据框sampled_relaDF，其中pair1和pair2两列内容都是sampled_ids的子集，但是现在想赋值新的数据框allpairs_relaDF，# 列名和sampled_relaDF相同，其中pair1和pair2列为sampled_ids元素，要求pair1元素所在顺序比pair2的小，# 如果这一对组合在sampled_relaDF没有相应的行就添加这一行并且给除pair1/2的其他几列赋值NA，如果有的话就直接用sampled_relaDF对应的行all_combinations <- expand.grid(pair1 = sampled_ids,                               pair2 = sampled_ids,                               stringsAsFactors = FALSE) %>%  filter(pair1 < pair2) %>%  arrange(pair1, pair2)# 获取原始数据框的其他列名（除了pair1和pair2）other_cols <- setdiff(names(sampled_relaDF), c("pair1", "pair2"))# 合并原始数据并填充"un"allpairs_relaDF <- all_combinations %>%  left_join(sampled_relaDF, by = c("pair1", "pair2")) %>%  mutate(across(all_of(other_cols), ~ ifelse(is.na(.), "un", .))) %>%  select(names(sampled_relaDF))write.csv(allpairs_relaDF,paste0('./',usesite,'_addun_relaDF_deg_re12_re123',".csv"), row.names = F, fileEncoding = "UTF-8")####### __1.3 bar_plot统计deg #######most_deg <- length(unique(allpairs_relaDF$deg))-1allpairs_relaDF$deg <- factor(allpairs_relaDF$deg,                                 levels=c(paste0('deg',1:most_deg),'un')                                   )allpairs_relaDF$re1_deg2more <- factor(allpairs_relaDF$re1_deg2more,                                 levels=c('pc','sib',paste0('deg',2:most_deg),'un')                                   )allpairs_relaDF$re12_deg3more <- factor(allpairs_relaDF$re12_deg3more,                                 levels=c('pc','sib','gr','avu',paste0('deg',3:most_deg),'un')                                   )####### ____1.3.1 col: deg #######cut_use <- as.data.frame(table(allpairs_relaDF$deg))cut_low <- cut_use %>%  filter(Var1 != 'un') %>%  summarise(max_freq = max(Freq)) %>%  pull(max_freq)+20cut_high <- cut_use$Freq[cut_use$Var1=='un']-50colors <- colorRampPalette(brewer.pal(7, "Set2"))(length(unique(allpairs_relaDF$deg)))p <- ggplot(allpairs_relaDF, aes(x = deg)) +  geom_bar(fill = colors) +  labs(title = "Kinship Statistics (deg)", x = "kinship", y = "number of pairs") +  theme_minimal()p <- p+scale_y_break(c(cut_low,cut_high),#截断位置及范围                space = 0.3,#间距大小                scales = 0.6)#上下显示比例，大于1上面比例大，小于1下面比例大p <- p + geom_text(stat = "count", aes(label = ..count..), vjust = -0.5)max_value <- cut_high + 150p <- p + scale_y_continuous(limits = c(0, max_value))# 显示图形ggsave(filename = "Kinship_Stat_deg.png", plot = p, width = 8, height = 6, dpi = 300)####### ____1.3.2 col: re1_deg2more #######cut_use <- as.data.frame(table(allpairs_relaDF$re1_deg2more))cut_low <- cut_use %>%  filter(Var1 != 'un') %>%  summarise(max_freq = max(Freq)) %>%  pull(max_freq)+20cut_high <- cut_use$Freq[cut_use$Var1=='un']-50colors <- colorRampPalette(brewer.pal(7, "Set2"))(length(unique(allpairs_relaDF$re1_deg2more)))p <- ggplot(allpairs_relaDF, aes(x = re1_deg2more)) +  geom_bar(fill = colors) +  labs(title = "Kinship Statistics (re1_deg2more)", x = "kinship", y = "number of pairs") +  theme_minimal()p <- p+scale_y_break(c(cut_low,cut_high),#截断位置及范围                space = 0.3,#间距大小                scales = 0.6)#上下显示比例，大于1上面比例大，小于1下面比例大p <- p + geom_text(stat = "count", aes(label = ..count..), vjust = -0.5)max_value <- cut_high + 150p <- p + scale_y_continuous(limits = c(0, max_value))# 显示图形ggsave(filename = "Kinship_Stat_re1_deg2more.png", plot = p, width = 8.5, height = 6, dpi = 300)####### ____1.3.3 col: re12_deg3more #######cut_use <- as.data.frame(table(allpairs_relaDF$re12_deg3more))cut_low <- cut_use %>%  filter(Var1 != 'un') %>%  summarise(max_freq = max(Freq)) %>%  pull(max_freq)+20cut_high <- cut_use$Freq[cut_use$Var1=='un']-50colors <- colorRampPalette(brewer.pal(7, "Set2"))(length(unique(allpairs_relaDF$re12_deg3more)))p <- ggplot(allpairs_relaDF, aes(x = re12_deg3more)) +  geom_bar(fill = colors) +  labs(title = "Kinship Statistics (re12_deg3more)", x = "kinship", y = "number of pairs") +  theme_minimal()p <- p+scale_y_break(c(cut_low,cut_high),#截断位置及范围                space = 0.3,#间距大小                scales = 0.6)#上下显示比例，大于1上面比例大，小于1下面比例大p <- p + geom_text(stat = "count", aes(label = ..count..), vjust = -0.5)max_value <- cut_high + 150p <- p + scale_y_continuous(limits = c(0, max_value))# 显示图形ggsave(filename = "Kinship_Stat_re12_deg3more.png", plot = p, width = 9, height = 6, dpi = 300)

(3.2) 亲缘关系对分类统计

####### Part2: mtYHap理论分组 ############## __2.1 取出mtYHap非落单者的子集 #######use_pc_pairdtfm <- pc_pairdtfm %>%   filter(pair2!=0)mt_pc_pairdtfm <- filter(use_pc_pairdtfm,sex2==2)Y_pc_pairdtfm <- filter(use_pc_pairdtfm,sex1==1 & sex2==1)####### __2.2 构建mt/Y_net ######## ！！！Notice:  父母->孩子是单向边nodes <- union(mt_pc_pairdtfm$pair1,mt_pc_pairdtfm$pair2)one_way_edges <- mt_pc_pairdtfm %>%  filter(relation == "pc_21") %>%  mutate(from = pair2,to = pair1) %>%  select(from,to,relation,sex1,sex2)mt_net <- graph_from_data_frame(one_way_edges, directed = TRUE)l <- layout_with_fr(mt_net)V(mt_net)$size <- 3E(mt_net)$width <- 5plot(mt_net, edge.arrow.size = 0.1)nodes <- union(Y_pc_pairdtfm$pair1,Y_pc_pairdtfm$pair2)one_way_edges <- Y_pc_pairdtfm %>%  filter(relation == "pc_21") %>%  mutate(from = pair2,to = pair1) %>%  select(from,to,relation,sex1,sex2)Y_net <- graph_from_data_frame(one_way_edges, directed = TRUE)l <- layout_with_fr(Y_net)V(Y_net)$size <- 3E(Y_net)$width <- 5plot(Y_net, edge.arrow.size = 0.1)####### __2.3 分组 ######## 已知mt_net是一个igraph对象，要求把其中不同连通分量对应的node的name属性按连通分量进行分组，最后得到一个数据框，# 一列叫SampleID即name属性赋值结果，一列叫mt_group，内容是按连通分量分组后从前往后自动编号为mtHap1/mtHap2等等# 获取连通分量编号####### ____2.3.1 mt取分组并给落单人续编 #######membership <- components(mt_net)$membership# 创建数据框mtHap_group <- data.frame(  SampleID = V(mt_net)$name,  mt_group = paste0("mtHap", membership)) %>%  filter(!(SampleID %in% unsampled_ids)) %>%  arrange(mt_group)mtHap_pedDF <- mtHap_group %>%  mutate(mt_group = factor(mt_group)) %>%  mutate(mt_group = paste0("mtHap", as.numeric(mt_group)))# 找到mtHap_pedDF中mt_group的最大值max_mt_group <- max(as.numeric(gsub("mtHap", "", mtHap_pedDF$mt_group)))# 为mt_left_ids创建新的mt_group编号mt_left_ids <- setdiff(sampled_ids,mtHap_pedDF$SampleID)mt_left_pedDF <- data.frame(  SampleID = mt_left_ids,  mt_group = paste0("mtHap", max_mt_group + 1:length(mt_left_ids)))# 将mt_left_pedDF添加到mtHap_pedDF中mtHap_pedDF <- rbind(mtHap_pedDF, mt_left_pedDF)####### ____2.3.2 Y取分组并给落单人续编  #######membership <- components(Y_net)$membershipYHap_group <- data.frame(  SampleID = V(Y_net)$name,  Y_group = paste0("YHap", membership)) %>%  filter(!(SampleID %in% unsampled_ids)) %>%  arrange(Y_group)YHap_pedDF <- YHap_group %>%  mutate(Y_group = factor(Y_group)) %>%  mutate(Y_group = paste0("YHap", as.numeric(Y_group)))# 找到YHap_pedDF中Y_group的最大值max_Y_group <- max(as.numeric(gsub("YHap", "", YHap_pedDF$Y_group)))# 为Y_left_ids创建新的Y_group编号male_ids <- intersect(int_pedDF$id[int_pedDF$sex == 'M'], sampled_ids)Y_left_ids <- setdiff(male_ids,YHap_pedDF$SampleID)Y_left_pedDF <- data.frame(  SampleID = Y_left_ids,  Y_group = paste0("YHap", max_Y_group + 1:length(Y_left_ids)))# 将Y_left_pedDF添加到YHap_pedDF中YHap_pedDF <- rbind(YHap_pedDF, Y_left_pedDF)write.csv(mtHap_pedDF,paste0('./',usesite,'_mtHap_group',".csv"), row.names = F, fileEncoding = "UTF-8")write.csv(YHap_pedDF,paste0('./',usesite,'_YHap_group',".csv"), row.names = F, fileEncoding = "UTF-8")

(3.3) 画树pedigree

####### Part3: plot pedigree with sample-info+mtYHap ######## https://cran.r-project.org/web/packages/pedtools/vignettes/pedtools.html####### __3.1 采样=实线，未采样=虚线 ######## 初始化 lty 向量，所有个体的线型默认为实线 (lty = 1)。将未采样个体的线型设置为虚线 (lty = 2)# lty 的作用通常是通过底层绘图函数或绘图包的特定参数来实现的，而不是直接传递给 plot 函数lty <- rep(1, length(all_ids))lty[all_ids %in% unsampled_ids] <- 2####### __3.2 mtHap=fill颜色 #######mtHap_as_fill <- rep(NA, length(all_ids))mtHap_colors <- colorRampPalette(brewer.pal(12, "Set3"))(length(unique(mtHap_pedDF$mt_group)))names(mtHap_colors) <- unique(mtHap_pedDF$mt_group)for ( mtHap_index in1:length(unique(mtHap_pedDF$mt_group)) ) {  this_mtHap <- unique(mtHap_pedDF$mt_group)[mtHap_index]  these_ids <- mtHap_pedDF$SampleID[mtHap_pedDF$mt_group == this_mtHap]  mtHap_as_fill[all_ids %in% these_ids] <- mtHap_colors[mtHap_index]}####### __3.3 YHap=border颜色 #######YHap_as_border <- rep('black', length(all_ids))YHap_colors <- colorRampPalette(brewer.pal(9, "Set1"))(length(unique(YHap_pedDF$Y_group)))names(YHap_colors) <- unique(YHap_pedDF$Y_group)for ( YHap_index in1:length(unique(YHap_pedDF$Y_group)) ) {  this_YHap <- unique(YHap_pedDF$Y_group)[YHap_index]  these_ids <- YHap_pedDF$SampleID[YHap_pedDF$Y_group == this_YHap]  YHap_as_border[all_ids %in% these_ids] <- YHap_colors[YHap_index]}####### __3.4 pedtools::plot #######png("Pedigree_with_theoretical_mtYHap.png", width = 40, height = 20, units = 'in', res = 150)# 划分绘图区域：左边 20% 用于图例，右边 80% 用于 plotlayout(matrix(c(1, 2), nrow = 1), widths = c(0.15, 0.85))# 绘制图例（左）par(mar = c(1, 0, 1, 0)) #上，左，下，右plot.new()legend("left",  legend = c("Sampled", "Unsampled"),   lty = c(1, 2),       col = "black", title = "Line Types", cex = 1.4, bty = "n")legend("center", legend = names(mtHap_colors), col = mtHap_colors, pch = 21,       pt.bg = mtHap_colors, cex = 1.5, bty = "n")legend("right", legend = names(YHap_colors), col = YHap_colors, pch = 21,       pt.bg = YHap_colors, cex = 1.5, bty = "n")# 绘制系谱图（右）par(mar = c(2, 0, 2, 10)) #上，左，下，右plot(pedDF,     fill = mtHap_as_fill,     col = YHap_as_border,     lwd = 5,     lty = lty)dev.off()

(4) pedtools所属系列的书

《Pedigree Analysis in R》，作者MagnusDehli Vigeland ，单位Department of Medical Genetics, Oslo University Hospital and University of Oslo, Oslo, Norway。可以当做包的doc。