使用geneHapR做单倍型分析
1、以hmp文件为例
chr1_hmp<- read.delim("data/Haplotype/chr1.hmp.txt", check.names=F, header=T)chr1_hmp<- chr1_hmp[,-c(1,5:11)]# 去掉无用的列chr1_hmp_df<- chr1_hmp %>% tidyr::separate(alleles, into=c("REF","Alt"), sep="/")# 将alleles拆分为REF、ALT列chr1_hmp_df<- data.frame(CHR=chr1_hmp_df$chrom, POS=chr1_hmp_df$pos, REF=chr1_hmp_df$REF, Alt=chr1_hmp_df$Alt, INFO="NA", chr1_hmp_df[,-c(1:4)])# 创建新的数据框hap_data<- chr1_hmp_df %>% dplyr::filter(CHR=="1", POS<=20000&POS>=1)# 取一定范围内的数据# 从表格形式的基因型数据开始单倍型分析hapResult<- table2hap(hap_data, hapPrefix="Hap", hetero_remove=TRUE,# 移除包含杂合位点的样本na_drop=TRUE)# 移除包含基因型缺失的样本hapSummary<- hap_summary(hapResult, hapPrefix="Hap")# 官方可视化热图plotHapTable(hapSummary,# 单倍型结果hapPrefix="Hap",# 单倍型名称前缀angle=45,# 物理位置的角度displayIndelSize=0,# 图中展示最大的Indel大小title="Haplotype", ALLELE.color="white")# 自定义热图hapSummary_df<- hapSummary[c(5:nrow(hapSummary)),]%>% dplyr::select(-c("Accession"))heatmap_data<- hapSummary_df %>%# 创建副本，保留原始频率as_tibble()%>% mutate(row_id=row_number())%>%# 转换为长格式pivot_longer(cols=-c(Hap, freq, row_id),# 排除这些列names_to="position_name",# 原始列名values_to="allele"# 等位基因)%>%# 从列名中提取位置信息mutate(# 提取位置数字（假设列名是pos1, pos2, pos3或类似格式）position=as.numeric(gsub("\\D+","", position_name)),# 提取数字# 如果列名没有数字，则使用列索引position=ifelse(is.na(position), as.numeric(factor(position_name, levels=unique(position_name))), position), hap=Hap)%>%# 排序和选择arrange(desc(freq), hap, position)%>% select(hap, position, allele, freq)%>% ungroup()# 创建热图p3<- ggplot(heatmap_data, aes(x=factor(position), y=reorder(hap, freq), fill=allele))+ geom_tile(color="white", size=0.8, width=0.9, height=0.9)+ geom_text(aes(label=allele), color="white", fontface="bold", size=4)+ scale_fill_brewer(palette="Set1")+ scale_y_discrete(labels=function(x){# 在y轴标签上添加频率freq_values<- heatmap_data %>% distinct(hap, freq)%>% arrange(desc(freq))%>% pull(freq)paste0(x," (", sprintf("%s", freq_values),")")})+ labs(title="Haplotype heatmap",#subtitle = "Visualization of alleles at each polymorphic site",x="Variant Position", y=NULL, fill="Allele")+ theme_minimal(base_size=12)+ theme(axis.text.x=element_text(face="bold", angle=45, color="black"), axis.text.y=element_text(face="bold", color="black"), panel.grid=element_blank(), legend.position="bottom", plot.title=element_text(face="bold", hjust=0.5), plot.subtitle=element_text(hjust=0.5, color="black"))+ guides(fill=guide_legend(nrow=1, byrow=TRUE))print(p3)# LD heatmapolor.rgb<- colorRampPalette(rev(c("white","red")),space="rgb")assign("name_gap",0.2, envir=.GlobalEnv)p<- plot_LDheatmap(hap=hapResult,# 单倍型结果add.map=FALSE,# 是否添加基因模式图gff=NULL,# 注释信息Chr=NULL,# 染色体名称start=NULL,# 基因的起始位置end=NULL, title=NULL, color=olor.rgb(20), newpage=F)LDheatmap::LDheatmap(gdat=p, title=NULL, text=FALSE,# 单元格不显示数字SNP.name=F)