如何重命名Beagle文件染色体位置列并匹配fai索引

问题背景 #

我有多个制表符分隔的文本文件，所用程序不允许染色体名称包含多个下划线（比如NC_044592.1_3795会导致解析错误），因此需要重命名Beagle文件（MM.beagle.gz）中的marker列（原染色体-位置列）。

现有文件内容 #

基因组索引文件 my.fna.fai #

执行命令：

head my.fna.fai

输出内容：

NC_044571.1     115307910       88      80      81
NC_044572.1     151975198       116749435       80      81
NC_044573.1     113180500       270624411       80      81
NC_044574.1     71869398        385219756       80      81

Beagle文件 MM.beagle.gz #

执行命令：

zcat MM.beagle.gz | head | cut -f 1-3

输出内容：

marker  allele1 allele2
NC_044571.1_3795        G       T
NC_044573.1_3796        G       T
NC_044572.1_3801        T       C
NC_044574.1_3802        G       A

需求说明 #

• 将.fai文件中的所有contig名称替换为从1开始的连续数字编号（即1:nrow(my.fna.fai)），修改后格式如下：

  1     115307910       88      80      81
  2     151975198       116749435       80      81
  3     113180500       270624411       80      81
  4     71869398        385219756       80      81
  

• 同步修改Beagle文件的marker列：将每个标记中的原始contig名称替换为对应数字编号，格式为[数字编号]_[位置]，修改后前3列如下：

  marker  allele1 allele2
  1_3795        G       T
  3_3796        G       T
  2_3801        T       C
  4_3802        G       A
  

• 即使Beagle文件的marker列顺序打乱，也要保证contig名称与.fai文件的映射完全匹配（比如.fai中NC_1234.1对应编号142，那么Beagle中所有NC_1234.1_XXX都要改为142_XXX）。
• 由于Beagle文件超过210GB，优先采用Bash解决方案，R处理大文件不适用。

可行解决方案 #

1. 生成contig到数字编号的映射表 #

awk 'BEGIN{OFS="\t"}{print $1, NR}' my.fna.fai > my.fna.fai.nr

2. 拆分Beagle文件的marker列并保留原始行顺序 #

awk -F'\t' -v OFS='\t' '{split($1,a,"_"); print $0,a[1]"_"a[2],a[3]}' MM.beagle.txt | awk 'NR!=1 {print}' | awk 'BEGIN{OFS="\t"}{print $0, NR}' > file2.sep.txt

3. 排序文件以用于关联 #

sort file2.sep.txt > file2.1.s.txt

4. 关联映射表并生成修改后的Beagle内容 #

join -1 4 -2 1 file2.1.s.txt my.fna.fai.nr | sort -k6 -n | awk 'BEGIN{OFS="\t"}{$1=$2=$6=""; print $7"_"$5,$0}' | awk 'BEGIN{OFS="\t"}{$4=$5=""; print $0}' > file4.txt

5. 拼接表头和修改后的内容 #

echo $(awk 'NR==1 {print}' MM.beagle.txt); cat file4.txt

最终输出结果 #

marker allele1 allele2
1_3795  G       T
3_3796  G       T
2_3801  T       C
4_3802  G       A

修改.fai文件的单独命令 #

awk 'BEGIN{OFS="\t"}{print NR, $2, $3, $4, $5}' my.fna.fai

内容的提问来源于stack exchange，提问作者M. Beausoleil