如何重命名Beagle文件中的chromosome_position列并将其与索引fai匹配

Question

我有一些文本文件(以制表符分隔(，它们有不同的列。我需要重命名我的chromosome_position列(请参阅下面的MM.beagle.gz文件(，因为我使用的程序不允许在染色体名称中使用多个下划线(由于NC_044592.1_3795不能作为名称使用，导致解析问题(。

我的索引基因组如下：

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

bealgle文件如下所示：

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

在R中，我可以得到染色体和位置：

beag = read.table("MM.beagle.gz", header = TRUE)
chr=gsub("_\d+$", "", beag$marker)
pos=gsub("^[A-Z]*_[0-9]*.[0-9]_", "", beag$marker)

但我无法在适当的位置重命名beagle文件。我想从1:nrow(my.fna.fai)重命名.fai文件中的所有contig，并将其与beagle文件匹配。

所以最后.fai应该看起来像：

head my.fna.fai

期望输出：

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

比格犬档案：

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

会给出：

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

其中22_3795是重叠群22和用_分隔的位置3795的级联。

解决方案将优先在bash中，因为R不实用，因为我的最终压缩beagle文件的文件大小很大(>210GB(

有人提议用这个改变.fai：

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

我现在不能弄清楚的是确保.fai和.beagle文件彼此一致。例如，如果.beagle文件的第一列(标记(被打乱，则应该可以将其与.fai文件匹配，并重命名.beagle文件中的染色体名称。例如，如果.fai中的NC_1234.1重命名为142，则.beagle中的所有NC_1234.1_XXX都应变为142_XXX，其中XXX是数字。

以下是解决方案的尝试：

awk 'BEGIN{OFS="t"}{print $1, NR}' my.fna.fai > my.fna.fai.nr
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
sort file2.sep.txt > file2.1.s.txt
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
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

Answer 1

为了确保新的FASTA索引和修改后的Beagle文件一致，我们可以使用FASTA索引以及关联数组来存储染色体名称及其行号。这样我们就可以解析Beagle文件，并使用染色体名称从数组中检索行号。这里有一种使用awk的方法：

rename_chroms.awk:的内容

BEGIN {
FS=OFS="t"
}
FNR==NR {
arr[$1]=NR
next
}
FNR==1 {
print
next
}
{
n = split($1, a, "_")
chrom = substr($1, 0, length($1) - (length(a[n]) + 1))
pos = a[n]
print arr[chrom] "_" pos, $2, $3
}

运行使用：

awk -f rename_chroms.awk my.fna.fai <(zcat MM.beagle.gz)

结果：

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