mjhjmtu

Building on: Extract names from File_B having overlapping intervals with File_A

I want to extract gene names (usually column 10, what's after "Name=") by matching the first column of file_b to file_a, and extracting the gene names if the second column of file_b lies within the gene interval, delineated by columns 4 and 5 of file_b. The first columns must match, such that I only get one gene per row (file_b) but I could in theory have multiple adjacent rows (column_b) match the same gene (e.g. if the second row in file_b was MT 4065)

There are a few problems with the code I have now.

(1)the way it is set up below, the last row of file_b is lost from the output, although this changes if this row (groupVII 17978350) is brought up the list. I would like it to work the way it is set up.

(2)The names are truncated if they have special characters (eg colons and hyphens). I would like to have the entire name after the equal sign.

(3) I would like to match the entry/row of file_b to the gene hits in the output, such that the first two columns is the entry and the third column is the gene hit.

file_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

file_b.tsv

MT 4050
groupIII 7332350
groupIV 5347350
groupVI 11230375
groupVII 17978350

code:

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

output.tsv

mt
si
zgc
bnip4

Desired_output

MT 4050 mt-nd2
groupIII 7332350 si:dkeyp-68b7.10
groupIV 5347350 zgc:153018
groupVI 11230375 bnip4
groupVII 17978350 si:ch211-284e13.4

# 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

Then run

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

Please make sure that the fields in your .tsv files are separated by tabs before processing them. My output:

MT 4050 mt-nd2
groupIII 7332350 si:dkeyp-68b7.10
groupIV 5347350 zgc:153018
groupVI 11230375 bnip4
groupVII 17978350 si:ch211-284e13.4

For the ID MT, the gene hit should be mt-nd2 not mt-nd1.

I still recommend using Python for data processing.

Your expected shown output doesn't look to me consistent(2 lines-->1st and 3rd), if that is a typo then could you please try following.

awk 'FNR==NRa[$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