Twease/FastaTwease

From Icbwiki

Contents 1 Searching biological sequences with information retrieval approaches 1.1 Corpus preparation 1.2 First sequence query with BM25 1.3 Modify the method 1.4 Next test: TAS1R3 1.5 With rat TAS1R3 1.6 Considering all terms 1.7 TAS1R3 human (again) 1.8 TAS1R3 rat (again) 1.9 human ENSP00000298762 novel protein 1 TM 1.10 human dopamin receptor D1 2 Formal evaluation 3 Method development

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.

Considering all terms

The approach is changed to include all the 3-grams instead of the just the 100 less frequent.

TAS1R3 human (again)

 102	Q0	41031	0	1682.7220721696349	twease
 102	Q0	41032	1	1682.7220721696349	twease
 102	Q0	8748	2	485.30564056151843	twease
 102	Q0	8749	3	464.95326073290676	twease
 102	Q0	3530	4	459.53838213925894	twease
 102	Q0	3529	5	459.4655007399344	twease
 102	Q0	3532	6	458.6147755418392	twease
 102	Q0	3533	7	458.6147755418392	twease
 102	Q0	37027	8	436.85086103161836	twease
 102	Q0	10132	9	434.79656462374174	twease

Index 8748 is ENSG00000008710

TAS1R3 rat (again)

 103	Q0	41031	0	892.2105666838777	twease
 103	Q0	41032	1	892.2105666838777	twease
 103	Q0	8748	2	465.38689433912015	twease
 103	Q0	3530	3	449.603795256827	twease
 103	Q0	8749	4	449.2513348708423	twease
 103	Q0	3532	5	448.6877217481474	twease
 103	Q0	3533	6	448.6877217481474	twease
 103	Q0	3529	7	447.4830146476012	twease
 103	Q0	37027	8	434.8688649287805	twease
 103	Q0	10132	9	432.78236120582835	twease

The results are now much more comparable between the human and rat search. Also, the match to LRP2 is gone from the 100 best matches in each topic.

human ENSP00000298762 novel protein 1 TM

The matches are:

 104	Q0	19367	0	562.7758519492568	twease
 104	Q0	31035	1	108.75911540085959	twease
 104	Q0	31036	2	108.75911540085959	twease
 104	Q0	27700	3	101.82761853863441	twease
 104	Q0	11763	4	101.00175274089645	twease
 104	Q0	11764	5	100.3089591542421	twease
 104	Q0	12482	6	99.93174646468118	twease
 104	Q0	12483	7	99.51344641663302	twease
 104	Q0	12481	8	99.46956241499163	twease
 104	Q0	40275	9	98.69795735774433	twease

human dopamin receptor D1

The matches are:

 108	Q0	34812	0	1154.2163464235011	twease
 108	Q0	35259	1	461.1347114631305	twease
 108	Q0	35260	2	461.1347114631305	twease
 108	Q0	39426	3	258.94310514309365	twease
 108	Q0	45070	4	257.96947802567473	twease
 108	Q0	34016	5	250.88603680618343	twease
 108	Q0	28008	6	247.62729392202303	twease
 108	Q0	46961	7	246.05595526861794	twease
 108	Q0	46960	8	245.46204881341686	twease
 108	Q0	48184	9	243.83515621691683	twease

35259 is ENSG00000197004 a short protein with 'Rhodopsin-like GPCR superfamily' domain.

35260 is ENSP00000306129 is DRD5 human.

39426 is ENSP00000366426 a novel protein with the following domain annotations:

 IPR002233 	Adrenergic receptor 
 IPR000276 	Rhodopsin-like GPCR superfamily 
 IPR000497 	Dopamine 1B receptor

Formal evaluation

Evaluating FastaTwease objectively will require a formal evaluation protocol. Ideas for such an evaluation are listed on a separate page.

Method development

• Can we calculate P-values for the results, instead of a score? some observations.