算法：hmp文件转化与表型匹配

引言

分析过程中，如果已经得到了hmp文件，下一步是将表型数据与hmp中的基因型数据一一对应，保证两者的样品ID信息一致，还需要对数据的格式进行规范化处理，用于后续的GWAS分析。

本文提供一种算法，能够实现对hmp文件和表型数据的关联筛选与校正。

主要步骤与设计思路

• 读取hmp文件和表型数据
• 替换hmp文件中的染色体编号格式
• 两表关联后迭代提取匹配的观测值
• 基因型和表型文件整理

项目运行环境

• centos7 linux
• R4.2.3

具体操作步骤

加载R包与数据

library(tidyverse)

chr_ref <- read.table("01_scripts/chr_num2str.txt",header = T)
df <- read_table(paste0("04_hmp/gene_",job,".hmp.txt"),show_col_types = F)
trait <- read_table(paste0("05_trait/","trait.txt"),show_col_types = F)

读取三个数据文件，其中第一个是染色体ID个不同格式对应信息，第二个是基因型hmp.txt文件，第三个是表型数据文件。

染色体格式转换

• chr_id_translate 函数

chr_id_translate <- function(data,type){
  # 输入俩参，一为原始数据，二为类型
  if (type == "1_to_chr1A"){
    # 数字转字符型
    old_id <- as.character(data)
    for (k in 1:nrow(chr_ref)){
      if (as.character(chr_ref$chr_num[k]) == old_id){
        return(chr_ref$chr_str[k])
      }
    }
  }else{
    if (type == "chr1A_to_1"){
      # 字符转数字型
      old_id <- as.character(data)
      for (k in 1:nrow(chr_ref)){
        if (as.character(chr_ref$chr_str[k]) == old_id){
          return(chr_ref$chr_num[k])
        }
      }
    }else{
      if (type == "1_to_1A"){
        old_id <- as.character(data)
        for (k in 1:nrow(chr_ref)){
          if (as.character(chr_ref$chr_num[k]) == old_id){
            new <- paste0(chr_ref$atom7[k],chr_ref$atom3[k],sep="")
            return(new)
          }
        }
      }else{
        print("Please input again! type inaviably")
      }
    }
  }
}

该函数提供了一种对染色体格式的快速转换方法，可以对数字型、字符型、全称之间进行快速转换，第一个参数是原始的编号，第二个参数选择转换方式，返回值是一个新的染色体编码值。

• 批量替换

for (i in 1:nrow(df)){
  df$chrom[i] <- chr_id_translate(
  df$chrom[i],type = "1_to_1A")
}

通过迭代将所有的数值型染色体编号换成数字加字母型。

基因型和表型匹配筛选

• 数据转换与处理

df2 <- rbind(colnames(df),df)
df_gene <- t(df2)
df_add_gene <- matrix(ncol = ncol(df_gene))
df_add_gene <- df_add_gene[-1,]
df_add_trait <- matrix(ncol = ncol(trait))
df_add_trait <- df_add_trait[-1,]
df_gene <- as.data.frame(df_gene)

对原始数据进行转置，目的是为了让基因型中样品ID按行排布，方便后续筛选，定义一个新的数据框用于储存迭代输出信息。

• 迭代提取匹配观测值

for (i in 1:nrow(df_gene)){
  id_gene <- df_gene$V1[i]
  for (k in 1:nrow(trait)){
    id_trait <- trait$ID[k]
    if (id_gene == id_trait){
      my_gene <- df_gene[i,]
      my_trait <- trait[k,]
      df_add_gene <- rbind(df_add_gene,my_gene)
      df_add_trait <- rbind(df_add_trait,my_trait)
    }else{
      next
    }
  }
}

通过上述方法可以找出两个表格中完全匹配的样品，生成的df_add_gene是所有匹配到的基因型文件，df_add_trait是所有对应的表型文件。后续可以直接拿来做GAPIT分析。

结果输出与保存

out_gene <- rbind(df_gene[1:11,],df_add_gene)
out_genet <- t(out_gene)
gene_final <- as.data.frame(out_genet)
write.table(gene_final,paste0("./06_out_gene/",job,".gene.hmp.txt"),
            quote = F,sep = "\\t",col.names = F,row.names = F)
trait_final <- as.data.frame(df_add_trait)

write.table(trait_final,paste0("./07_out_trait/",job,".trait.txt"),
            quote = F,sep = "\\t",col.names = T,row.names = F)
print(paste0(job," hmp and trait formate finished!"))

重新合并头文件并转置，恢复原有结构，然后分别将两个结果保存到对应文件夹中。