Twease/FastaTwease
From Icbwiki
Contents |
Searching biological sequences with information retrieval approaches
Corpus preparation
A version of Twease that can search a corpus made of biological sequences.
I made a corpus with the ensembl human proteins (Homo_sapiens.NCBI36.41.pep.all.fa) and indexed with Textractor/MG4J with overlapping 3-grams as words.
First sequence query with BM25
Searched this corpus with BM25 scoring and disjunctive queries (i.e., AAP | HIL| ...) and an example query chosen at random among the sequences in the indexed corpus:
>ENSP00000331971 pep:known chromosome:NCBI36:21:30774235:30774507:-1 gene:ENSG00 000184351 transcript:ENST00000331764 MSHYGSYYGGLGYSCGGFGGLGYGYGCGCGSFCRRGYGYGSGFGSYGYGSGFGGYGYGSG
At this stage, I could not build a doc store, so the index of the sequence in the input fasta file is used as document identifier (/PMID) in the trec_eval output.
ENSP00000331971 has index 132-1=131 in the input Fasta file (sequence numbers start at 0 with textractor, but grep numbers lines from 1).
Fabien Campagne@PC120871 /cygdrive/d/dev/ensembl-data grep '>' Homo_sapiens.NCBI36.41.pep.all.fa|grep -n ENSP00000331971 132:>ENSP00000331971 pep:known chromosome:NCBI36:21:30774235:30774507:-1 gene:EN SG00000184351 transcript:ENST00000331764
The best matches in trec_eval format are:
100 Q0 131 0 42.46582283201506 twease 100 Q0 148 1 28.144966050948458 twease 100 Q0 133 2 25.34762250828596 twease 100 Q0 135 3 23.00047774168498 twease 100 Q0 141 4 22.69864217728251 twease 100 Q0 134 5 22.67280744193289 twease 100 Q0 139 6 22.358672780465746 twease 100 Q0 144 7 22.312434886092404 twease 100 Q0 132 8 21.778693703352136 twease 100 Q0 136 9 21.145654858966886 twease
Therefore, the query sequence is found at rank 0 (best). That's a good start. The next best match is the sequence at index 148 in the input file.
$ grep '>' Homo_sapiens.NCBI36.41.pep.all.fa |head -149 |tail -1 >ENSP00000335566 pep:known-ccds chromosome:NCBI36:21:31049328:31049567:-1 gene:ENSG00000187005 transcript:ENST00000335093 CCDS13606.1 Verification: $ grep '>' Homo_sapiens.NCBI36.41.pep.all.fa|grep -n ENSP00000335566 149:>ENSP00000335566 ...
ENSP00000335566 is a novel gene, but marked as belonging to the keratin family (ENSF00000003793 : KERATIN).
The query sequence is Keratin-associated protein 21-1 marked as a member of the keratin family (ENSF00000003793).
Modify the method
Log inspection indicates that Twease uses only 8 terms in the query (the terms with the lowest corpus frequency). I changed the behavior to include the 100 terms with the lowest corpus frequency. This should help find more matches.
Next test: TAS1R3
The sweet taste receptor should be an interesting test. It has low homology to other GPCRs and human proteins.
The following matches are returned:
102 Q0 41031 0 336.1844751538954 twease 102 Q0 41032 1 336.1844751538954 twease 102 Q0 4167 2 53.182348855402594 twease 102 Q0 14048 3 50.04884551615295 twease 102 Q0 2071 4 49.95855069089579 twease 102 Q0 2072 5 49.95855069089579 twease 102 Q0 39985 6 49.741951152235266 twease 102 Q0 5039 7 49.5103894224437 twease 102 Q0 40045 8 48.447158334399354 twease 102 Q0 39904 9 47.41193634057475 twease 102 Q0 39905 10 47.41193634057475 twease
Index 41031 is ENSG00000169962/ENST00000339381/ENSP00000344411, TAS1R3 itself.
Index 4167 is ENSG00000142449/ENST00000270509/ENSP00000270509, FBN3 Fibrilin 3. The link is not obvious.
Index 14048 is ENSG00000184613/ENST00000333837/ENSP00000327988, NELL2, Protein kinase C-binding protein NELL2 precursor (NEL-like protein 2) (Nel-related protein 2)
With rat TAS1R3
Let's check that the rat TAS1R3 finds the human one:
103 Q0 41031 0 147.95588360052815 twease 103 Q0 41032 1 147.95588360052815 twease 103 Q0 21081 2 50.66315874422973 twease 103 Q0 21082 3 50.66315874422973 twease 103 Q0 21083 4 50.257226358209856 twease 103 Q0 40045 5 45.977610216284475 twease 103 Q0 36374 6 44.11439908886646 twease 103 Q0 4962 7 43.908705176987375 twease 103 Q0 32389 8 43.60676044312044 twease 103 Q0 6164 9 43.29170704576298 twease
It does, but most of the other matches do not overlap with the initial human search. An exception is sequence at index 40045 (rank 5 in rat search, rank 8 in human search).
Index 40045 is ENSG00000081479, LRP2.
LRP2 occured also in previous searches. This strongly suggests that something is wrong.
