# Picking illumina index sequence combinations

June 24, 2013

Illumina uses a short read with a separate sequencing primer to read index sequences contained in adapter sequences. If there is not enough heterogeneity in the index sequences the read quality suffers. This script generates combinations of indices that maximize heterogeneity given a set of allowed indices and the number of samples to be combined in a lane.

The script shown at the bottom ranks all possible combinations of n index sequences from the given list of index sequences and displays the best m combinations. This is done by (1) determining what fraction of all pairs of indices in a set has 4 or more differences and (2) what fraction of positions in the index has close to optimal sequence diversity (measured as being within 90% of maximal sequence entropy). Repository on bitbucket.

Given the following inputs

cat > setB.adapters <<EOF
1|ATCACG
3|TTAGGC
8|ACTTGA
9|GATCAG
10|TAGCTT
11|GGCTAC
20|GTGGCC
21|GTTTCG
22|CGTACG
23|GAGTGG
25|ACTGAT
27|ATTCCT
EOF


#./goodadapter_combinations.py n m
#    n: number of samples (indices) per lane
#    m: number of highest ranking combinations to return


and returns the following results to stdout:

***************************  ADAPTER SET  ****************************
1       ATCACG
3       TTAGGC
8       ACTTGA
9       GATCAG
10      TAGCTT
11      GGCTAC
20      GTGGCC
21      GTTTCG
22      CGTACG
23      GAGTGG
25      ACTGAT
27      ATTCCT
max sequence entropy at each cycle: 1.91829583405
top  sequence sets sorted by
(1) the fraction of positions with entropy > 90% of max. entropy
(2) the fraction of pairs with 4 or more differences
**********************************************************************
---------------
100% of positions had close to max. entropy
100% of pairs had more than 4 differences
---------------
3       TTAGGC
8       ACTTGA
10      TAGCTT
11      GGCTAC
22      CGTACG
27      ATTCCT
---------------
100% of positions had close to max. entropy
100% of pairs had more than 4 differences
---------------
3       TTAGGC
8       ACTTGA
10      TAGCTT
11      GGCTAC
22      CGTACG
23      GAGTGG
---------------
100% of positions had close to max. entropy
100% of pairs had more than 4 differences
---------------
3       TTAGGC
8       ACTTGA
9       GATCAG
10      TAGCTT
11      GGCTAC
22      CGTACG
---------------
100% of positions had close to max. entropy
93% of pairs had more than 4 differences
---------------
3       TTAGGC
8       ACTTGA
10      TAGCTT
11      GGCTAC
22      CGTACG
25      ACTGAT


Below is the original version of the python script. Current versions are available from the bitbucket repo.

#! /usr/bin/env python
"""
Arguments:
n:  size of the desired adapter sets (i.e. how many samples to multiplex)
m:  how many adapter combinations to show
Ranks adapter combinations by (1) how many positions/cycles
have close to maximal sequnce entropy (90%) (2) how
many adapter pairs have 4 or more differences.
"""

import sys
import itertools
import array
import math

def entropy_ratio(seq_list, max_h):
hf = []
for s in itertools.izip(*seq_list):
hf.append(entropy(s) / max_h)
return hf

def entropy(s):
nf = float(len(s))
h = 0
for nt in ["A", "C", "G", "T"] :
f = s.count(nt) / nf
if f > 0:
h += - f * math.log(f, 2)
return h

def dist(a, b):
return sum(1 for a, b in itertools.izip(a, b) if a != b)

def pairwise_diff(seq_list):
d = []
for pair in itertools.combinations(seq_list, 2):
d.append(dist(*pair))
return d

def max_entropy(n):
max_h = 0
for s in itertools.combinations_with_replacement("ACGT", n):
h = entropy(s)
if h > max_h:
max_h = h
return max_h

combo_list = []
for combo in itertools.combinations(kit, n):
seq_list = [b for a, b in combo]
hf = entropy_ratio(seq_list, max_h)
f_hf = sum(1 for x in hf if x > 0.9) / float(adapter_length)
d = pairwise_diff(seq_list)
f_d = sum(1 for x in d if x >= 4) / float(len(d))
combo_list.append((f_hf, f_d, combo))
combo_list.sort(key = lambda x: (x[0], x[1]))
combo_list.reverse()
return combo_list

if __name__ == "__main__":
if (len(sys.argv) != 3):
print >>sys.stderr, __doc__
sys.exit(1)
n = int(sys.argv[1])
m = int(sys.argv[2])
adapters = [x.strip().split("|") for x in sys.stdin]
print >>sys.stderr, "adapters of differing lenghts found in input"
sys.exit(1)
max_h = max_entropy(n)
print "\n".join("{}\t{}".format(a, b) for a, b in adapters)
print "max sequence entropy at each cycle: {}".format(max_h)
print "top {} sequences sorted by ".format(m)
print "  (1) the fraction of positions with entropy > 90% of max. entropy"
print "  (2) the fraction of pairs with 4 or more differences"
print "*" * 70