Adding Electronic Markers to FPC Maps

Supplemental information to provide detailed information on the data sets and experiments used to compare different strategies of electronic marker addition, using the search tools MegaBLAST (Zhang et. al., 2002), BLAST (Altschul et. al., 1997) and BLAT (Kent, 2001).

To add markers to a map, one uses an alignment program to return hits of the marker sequences against sequence associated to the clones in the map (this can be draft, finished, or bac-end sequence (BES)). Two primary decisions must be made: which program to use, and how to rank the resulting hits. These are discussed in sections I and II below, for a target of draft or finished sequence. BES targets are covered in section III.

In order to compare different programs and hit-ranking approaches, we have used rice markers provided by the Japanese Rice Genome Research Project(JGRP). The set consisted of 3444 ESTs from 2191 cDNAs, and 338 STSs with 263 unique ids. The average EST sequence length was 421, while the STS average was 325.

Our target database of clone sequences was the full set of draft sequence for rice, as of Feb. 2003, consisting of 11126 sequenced contigs from 3497 clones. These were downloaded from Genbank with the query

    Oryza [ORGN] AND (30000[SLEN]:250000 [SLEN]) AND ((htg [KYWD] OR BAC [ALL] OR 
    chromosome [TITL] OR PAC [ALL]) NOT (marker [TITL] OR EST [TITL] OR mRNA [TITL]
    OR RAPD [TITL] OR GSS [KYWD] OR telomere [TITL] OR protein[TITL]))

Our BES database consisted of 127457 rice BESs of clones in the rice FPC map (Chen et al. 2001; Mao et al. 2000).

We scored hits as true-positive or false-positive, depending on whether the hit was to the correct chromosome, where the chromosome of the marker was known from the JGRP map, while that for the clone was known from the Genbank record (for draft sequence) or the FPC map (for BES). This is an admittedly simplistic measure, but if bad hits are randomly distributed, it should only mislabel a bad hit as good 1/12 of the time (since there are 12 chromosomes). To check whether bad hits might be concentrated on the "correct" chromosome, 10 hits labeled as true-positive were examined manually, using the detailed map location of the given marker and neighboring markers, and in all cases the hit was confirmed.

Another consideration is that the markers are not entirely specific; i.e., they may hybridize in the lab to chromosomes other than the one mapped; if these extra hits are also found electronically, they will be scored as false-positives according to our "wrong chromosome" test. Because of this, the error rate does not go to zero at high stringencies, but goes to a "base rate" which is determined by the specificity of the markers.Only error in excess of this base rate can be attributed to the search strategy. For the EST markers in our test set this rate is about 9%, while for the STS markers it is much lower.

I. Scoring Strategy: %Match vs. E-value

In a traditional BLAST search, hits are ranked according to the BLAST E-value, which measures the probability of that hit occurring by chance. This is a sensible measure to use for similarity searches, where one wants to find any target sequences that probably have the same function or evolutionary origin as the target.

When adding electronic markers, however, one isn't looking for similarity, but for identity. The marker sequence is supposed to be exactly present in at least one location in the genome. A perfect match of, say, half of the marker sequence, is not the kind of hit one is looking for, even though it could have a very low E-value.

For this reason, it seems sensible to score electronic marker hits according to the percentage of the total marker length which is matched; this score will be denoted "%Match". (Note that this is not the same as the "percent identity" from BLAST reports, because that score only gives the fraction of bases that match within the alignment region, which may not cover the whole marker.)

The %Match strategy is clearly seen to be superior to E-value in tests on the STS subset of JRGP markers (table I below). On this dataset, using MegaBLAST, a screen on %Match ≥ 95 finds 462 hits, with a false-positive rate of 2%;by contrast, taking E-value ≥ 170 gives only 179 hits, with 3% false-positives.[ %Match was implemented by parsing the output file and taking %Match = (match length/query length) *(percent identity) .]

It is understandable that E-value should work poorly in this case, because the sequences have a variety of lengths, and a given E-value is only well-suited for a small range of lengths. E-170 corresponds on this dataset to a match of about 300bp, and since the average STS marker has length only 325, there are many which are shorter than 300. Hits of these markers will never be found at E-170; however, if the E-value cutoff is lowered to accept the shorter markers, it becomes too low for the longer ones, and false hits are introduced. For example, marker L704 has length 268, and scores a 100% hit on clone AC018929, but this hit has E-value only E-149. To accept this hit, the E-value cutoff would have to be lowered at least to E-149, at which point the cumulative false-positive rate is 6%, twice as high as at E-170.

Another problem with E-value screens, which may not be sufficiently appreciated, is that the value returned for a given alignment varies with the size of the databases being compared. The larger the databases from which the two sequences are drawn, the larger will be the returned E-value. This means that hits which are accepted at a given E-value cutoff, could be rejected if the database sizes are increased, even if the two sequences are unchanged. So for example if one runs an E-value screen of some markers against batches of sequence that come in periodically over time, and then runs the same screen later against the full set of sequences, the set of hits will be different. By contrast, the %Match strategy is independent of database size, so this potential confusion is absent.

From the data in table I, these STS marker sequences are seen to be very accurate, so that one can expect to find typically at least 95% of the marker sequence in the target genome; this is the ideal situation for the %Match strategy. If the marker sequences are less accurate, or contain bad reads at either end, then the advantage of %Match over E-value will be reduced or eliminated. This effect can be seen within the EST marker set, which also has the additional complication of spliced introns, as discussed in the next section.

It is worth mentioning that not all of the markers get a hit at %Match ≥ 95, or indeed at any cutoff. At 95%, hits are found for 220 unique markers, or 83% of the total 263 markers in the set. Lowering the cutoff to 90% brings the marker total to 235, while going all the way down to 70% brings in only 5 more for a total of 240. So it seems that about 10% of the markers simply do not hit, which is not unreasonable given the gaps in clone coverage and in draft sequence, and possible errors in marker sequencing.

Finally we note that although MegaBLAST data was used in this comparison, BLAT gives the same result; furthermore, in the context of BSS, the %Match scoring is supported only through BLAT (the reason for this, as discussed below, is that %Match requires BLAT in order to work well on EST data). MegaBLAST data was used here merely to emphasize that it is the scoring strategy, not the program, which is being tested.

II. Unsplicing Introns: BLAST vs. BLAT

Marker sequences which have undergone splicing present a problem for the %Match strategy, because although the full marker sequence may be present in the target, it is in disconnected exons. BLAST (and MegaBLAST) returns each exon hit as a separate high-scoring pair, but from the report it is not straightforward to put those pieces together to produce the total %Match.

However, there is a search tool which does this automatically - BLAT. Given this capability, and the superiority of the %Match strategy seen above on STS data, it is reasonable to believe that BLAT + %Match will be superior on the EST data. In fact we will see that it is, but the margin of difference is much smaller than on the STS data.