从 File_B 中提取与 File_A 具有重叠间隔的名称

从 File_B 中提取与 File_A 具有重叠间隔的名称

两个空格分隔的文件:

文件_A

MT 50000
groupI 7850000
groupI 7950000
groupI 9050000
groupI 21750000
groupII 8750000
groupII 10550000
groupII 16150000
groupII 20850000
groupIII 14750000
groupIII 15250000
groupIII 15450000
groupIII 15550000
groupIII 15650000
groupIV 7850000

第一列是组 ID,第二列是 100,000 单位长的区间的中点组内。例如第一行对应MT组中的区间1-100000,第二行对应区间7800000-7900000,依此类推。

文件_B

MT 2851 3825 Name=mt-nd1
MT 4036 5082 Name=mt-nd2
MT 5465 7015 Name=mt-co1
MT 7173 7863 Name=mt-co2
MT 8097 8780 Name=mt-atp6
groupI 18791 22890 Name=FGF12
groupI 36880 38991 Name=MB21D2
groupI 65279 68049 Name=cldn15lb
groupI 77722 105198 Name=col4a4
groupI 117583 141390 Name=col4a3
groupI 150455 155401 Name=sst1.1
groupI 9050030 9058000 Name=bco2b
groupI 1076088 1085084 Name=SORL1
groupI 1175505 1181937 Name=abcg4b
groupI 1184288 1184688 Name=lyrm9
groupI 1185206 1186192 Name=ift20

File_B的第1列是基因所在的组/染色体名称,第2列和第3列是基因的间隔,其中第2列是开始,第3列是结束。最后,第 4 列是基因名称。我想从File_B的第4列中提取间隔在File_A的100,000间隔内的唯一基因名称。

输出文件

mt-nd1
mt-nd2
mt-co1
mt-co2
mt-atp6
bco2b

我将以下代码用于不同但相似的过程(File_B 有更多列,File_A 的第二列是一个点而不是间隔)。

while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { if (gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) !~ /\s/) print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1); }' <File_B.txt; done < File_A.txt > Output_file.txt

答案1

假设组名称必须相同(从您的描述中不清楚,但数据和预期输出表明如此):

$ sort -k1,1 -k2n,2n <(awk '{print $1, $2-50000, $2+50000, $2}' File_A) File_B |
  awk '
    !gsub(/[^=]*=/, "", $4) {g=$1; s=$2; e=$3; m=$4; next}
    $2 > s && $3 <= e && $1 == g {if(m){print g, m; m=""} print "   "$4}
  '
MT 50000
   mt-nd1
   mt-nd2
   mt-co1
   mt-co2
   mt-atp6
groupI 9050000
   bco2b

没有标题:

$ sort -k1,1 -k2n,2n <(awk '{print $1, $2-50000, $2+50000}' File_A) File_B |
  awk '
    !gsub(/[^=]*=/, "", $4) {g=$1; s=$2; e=$3; next}
    $2 > s && $3 <= e && $1 == g {print $4}
  '
mt-nd1
mt-nd2
mt-co1
mt-co2
mt-atp6
bco2b

相关内容