How to search the human genome for sequences that differ from a given sequence by a set number of mismatches?

Question

[azhai]

Part of a program I'm writing involves taking in an input query sequence that is between 18-21 bp and finding all, or as many as possible, sequences that differ from that sequence by 1, 2, 3, and 4 bp. Basically any alignments with 1, 2, 3, or 4 mismatches. I've been trying to use blastn to do this so far, but the problem is that I'm only getting 1 hit, the actual place where the query sequence aligns. It seems that BLAST doesn't let you specify how many mismatches you want to allow in the sequence. Is there maybe a way to set the score cutoff for what is considered a hit to be lower? I'm also open to suggestions on other ways to do this not using BLAST, since I've experimented a lot with it and had little success. Thanks in advance for any advice people give!

Answer 1

[Mensur Dlakic]

BLAST operates with E-value thresholds. If there are multiple hits for your sequences that satisfy the E-value threshold, with or without mismatches, BLAST would identify them. I think the default threshold is E=10. Depending on the E-value of your perfect match, you may want to increase the threshold to 100 or so. It is possible that there are no hits with 1-4 mismatches that BLAST can identify with its default parameters, so you may have to play with word size as well. If I remember correctly, smaller word sizes will make the search slower but more sensitive to short matches.

Answer 2

[harold.smith.tarheel]

See this answer.

Answer 3

[Pierre Lindenbaum]

Simple and stupide solution in C. Not really tested.

	#include <stdio.h>
	#include <string.h>
	#include <stdlib.h>
	#include <ctype.h>
	#define MAX_SEQ_LENGTH 100
	static char compl(char c) {
	switch(c) {
	case 'A': return 'T';
	case 'T': return 'A';
	case 'G': return 'C';
	case 'C': return 'G';
	default:return 'N';
	}
	}
	int main(int argc,char** argv)
	{
	size_t i=0;
	size_t position=0;
	size_t qlen;
	char* query;
	size_t maxdiff=0;
	size_t win_size=0;
	char seqName[MAX_SEQ_LENGTH];
	int c;
	if(argc!=3) {
	fprintf(stderr,"usage %s query maxdiff < in.fa\n",argv[0]);
	return EXIT_FAILURE;
	}
	query=strdup(argv[1]);
	qlen = strlen(query);
	maxdiff = (size_t)atoi(argv[2]);
	if(qlen==0 ||maxdiff >= qlen) {
	fprintf(stderr,"check parameters");
	return EXIT_FAILURE;
	}
	for(i=0;i<qlen;i++) query[i]=toupper(query[i]);
	char *window = (char*)malloc((qlen+1)*sizeof(char));
	memset((void*)window,0,(qlen+1)*sizeof(char));
	seqName[0]=0;
	while((c=fgetc(stdin))!=EOF)
	{
	if(c=='>')
	{
	i=0UL;
	while((c=fgetc(stdin))!=EOF && c!='\n')
	{
	seqName[i++]=c;
	}
	seqName[i]=0;
	position=0;
	win_size=0;
	continue;
	}
	if(isspace(c)) continue;
	window[win_size++]=toupper(c);
	position++;
	if(win_size==qlen) {
	size_t count=0UL;
	char strand='+';
	for(i=0;i< qlen && count<= maxdiff;i++) {
	if(query[i]!=window[i]) count++;
	}
	if(count > maxdiff) {
	strand = '-';
	count=0UL;
	for(i=0;i< qlen && count<= maxdiff;i++) {
	if(query[i]!=compl(window[(qlen-1)-i])) count++;
	}
	}
	if(count <= maxdiff && i==qlen) {
	fprintf(stdout,"%s\t%d\t%c\t%s\t%s\t%d\n",seqName,position+1,strand,query,window,count);
	}
	memmove((void*)window,&window[1],sizeof(char)*(qlen-1));
	win_size--;
	}
	}
	return EXIT_SUCCESS;
	}

view raw biostar9528296.c hosted with ❤ by GitHub

compilation:

gcc -Wall -o biostar9528296 biostar9528296.c

usage:

./biostar9528296 ATATCATCGTACTAGCGATGTTTAGGGGA 2 < in.fa