建立在:从 File_B 中提取与 File_A 具有重叠间隔的名称
我想通过将 file_b 的第一列与 file_a 匹配来提取基因名称(通常是第 10 列,“Name=”之后的内容),如果 file_b 的第二列位于由第 4 列和文件_b.5第一列必须匹配,这样我每行只能得到一个基因(file_b),但理论上我可以让多个相邻行(column_b)匹配相同的基因(例如,如果 file_b 中的第二行是MT 4065
)
我现在的代码存在一些问题。
(1)按照下面的设置方式,file_b 的最后一行从输出中丢失,尽管如果该行 ( groupVII 17978350
) 出现在列表中,情况会发生变化。我希望它按照设置的方式工作。
(2) 如果名称中含有特殊字符(例如冒号和连字符),名称将被截断。我想在等号后面加上完整的名字。
(3) 我想将 file_b 的条目/行与输出中的基因命中相匹配,使得前两列是条目,第三列是基因命中。
文件_a.tsv
MT insdc gene 2851 3825 . + . ID=gene:ENSGACG00000020925 Name=mt-nd1 biotype=protein_coding description=NADH dehydrogenase 1%2C mitochondrial [Source:ZFIN%3BAcc:ZDB-GENE-011205-7] gene_id=ENSGACG00000020925 logic_name=mt_genbank_import version=1
MT insdc gene 4036 5082 . + . ID=gene:ENSGACG00000020929 Name=mt-nd2 biotype=protein_coding description=NADH dehydrogenase 2%2C mitochondrial [Source:ZFIN%3BAcc:ZDB-GENE-011205-8] gene_id=ENSGACG00000020929 logic_name=mt_genbank_import version=1
groupIII ensembl gene 7332324 7334769 . - . ID=gene:ENSGACG00000015265 Name=si:dkeyp-68b7.10 biotype=protein_coding description=si:dkeyp-68b7.10 [Source:ZFIN%3BAcc:ZDB-GENE-070912-667] gene_id=ENSGACG00000015265 logic_name=ensembl version=1
groupIV ensembl gene 1368026 1374881 . + . ID=gene:ENSGACG00000016447 Name=hnrnpa0b biotype=protein_coding description=heterogeneous nuclear ribonucleoprotein A0b [Source:ZFIN%3BAcc:ZDB-GENE-030131-6154] gene_id=ENSGACG00000016447 logic_name=ensembl version=1
groupIV ensembl gene 5347339 5349041 . - . ID=gene:ENSGACG00000017010 Name=zgc:153018 biotype=protein_coding description=zgc:153018 [Source:ZFIN%3BAcc:ZDB-GENE-060929-752] gene_id=ENSGACG00000017010 logic_name=ensembl version=1
groupV ensembl gene 120615 125489 . + . ID=gene:ENSGACG00000002103 Name=zdhhc6 biotype=protein_coding description=zinc finger%2C DHHC-type containing 6 [Source:ZFIN%3BAcc:ZDB-GENE-030131-3189] gene_id=ENSGACG00000002103 logic_name=ensembl version=1
groupVI ensembl gene 11230354 11232784 . + . ID=gene:ENSGACG00000009527 Name=bnip4 biotype=protein_coding description=BCL2 interacting protein 4 [Source:ZFIN%3BAcc:ZDB-GENE-051113-212] gene_id=ENSGACG00000009527 logic_name=ensembl version=1
groupVII ensembl gene 2271611 2277214 . + . ID=gene:ENSGACG00000019012 Name=sf3b2 biotype=protein_coding description=splicing factor 3b%2C subunit 2 [Source:ZFIN%3BAcc:ZDB-GENE-070928-1] gene_id=ENSGACG00000019012 logic_name=ensembl version=2
groupVII ensembl gene 15815857 15824549 . + . ID=gene:ENSGACG00000020296 Name=mpp1 biotype=protein_coding description=membrane protein%2C palmitoylated 1 [Source:ZFIN%3BAcc:ZDB-GENE-031113-4] gene_id=ENSGACG00000020296 logic_name=ensembl version=1
groupVII ensembl gene 17978322 17982388 . + . ID=gene:ENSGACG00000020399 Name=si:ch211-284e13.4 biotype=protein_coding description=si:ch211-284e13.4 [Source:ZFIN%3BAcc:ZDB-GENE-060526-161] gene_id=ENSGACG00000020399 logic_name=ensembl version=1
文件_b.tsv
MT 4050
groupIII 7332350
groupIV 5347350
groupVI 11230375
groupVII 17978350
代码:
while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { if (gensub(/.*Name=([A-Za-z0-9]*).*/, "\\1", 1) !~ /\s/) print gensub(/.*Name=([A-Za-z0-9]*).*/, "\\1", 1); }' <file_a.tsv; done < file_b.tsv > output.tsv
输出.tsv
mt
si
zgc
bnip4
期望输出
MT 4050 mt-nd2
groupIII 7332350 si:dkeyp-68b7.10
groupIV 5347350 zgc:153018
groupVI 11230375 bnip4
groupVII 17978350 si:ch211-284e13.4
答案1
# save this as script.awk or whatevernameyouwant.awk
function within_range(val, lower, upper, proximity) {
# you can specify the "proximity" as required
return val > lower - proximity && val < upper + proximity
}
BEGIN {
OFS="\t"
}
$1 == id && within_range(pos, $4, $5, 100) {
name = gensub(/.*Name=([^\t]*).*/, "\\1", 1)
if (name ~ /[^[:space:]]+/)
print id, pos, name
}
然后运行
while read -r id pos
do
awk -v id=$id -v pos=$pos -f script.awk file_a.tsv
done < file_b.tsv > output.tsv
请确保在.tsv
处理文件中的字段之前用制表符分隔它们。我的输出:
MT 4050 mt-nd2
groupIII 7332350 si:dkeyp-68b7.10
groupIV 5347350 zgc:153018
groupVI 11230375 bnip4
groupVII 17978350 si:ch211-284e13.4
对于ID来说MT
,基因命中应该是mt-nd2
不行的mt-nd1
。
我还是推荐使用Python来进行数据处理。
答案2
您预期的显示输出在我看来不一致(2 行 - >第 1 行和第 3 行),如果这是一个拼写错误,那么您可以尝试以下操作吗?
awk 'FNR==NR{a[$1]=$2;next} ($1 in a) && (a[$1]>=$4 && a[$1]<=$5){sub("Name=","",$10);print $1,a[$1],$10}' b.tsv a.tsv > output.tsv