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!"))
重新合并头文件并转置，恢复原有结构，然后分别将两个结果保存到对应文件夹中。
END
© 素材来源于网络，侵权请联系后台删除
笔记合集，点击直达 
