pysam/Pileup.py

   1 '''Tools for working with files in the samtools pileup -c format.'''
   2 import collections
   3 import pysam
   4
   5 PileupSubstitution = collections.namedtuple( "PileupSubstitution",
   6                                     " ".join( (\
   7             "chromosome",
   8             "pos",
   9             "reference_base",
  10             "genotype",
  11             "consensus_quality",
  12             "snp_quality",
  13             "mapping_quality",
  14             "coverage",
  15             "read_bases",
  16             "base_qualities" ) ) )
  17
  18 PileupIndel = collections.namedtuple( "PileupIndel",
  19                                       " ".join( (\
  20             "chromosome",
  21             "pos",
  22             "reference_base",
  23             "genotype",
  24             "consensus_quality",
  25             "snp_quality",
  26             "mapping_quality",
  27             "coverage",
  28             "first_allele",
  29             "second_allele",
  30             "reads_first",
  31             "reads_second",
  32             "reads_diff" ) ) )
  33
  34 def iterate( infile ):
  35     '''iterate over ``samtools pileup -c`` formatted file.
  36
  37     *infile* can be any iterator over a lines.
  38
  39     The function yields named tuples of the type :class:`pysam.Pileup.PileupSubstitution`
  40     or :class:`pysam.Pileup.PileupIndel`.
  41
  42     .. note::
  43        The parser converts to 0-based coordinates
  44     '''
  45
  46     conv_subst = (str,lambda x: int(x)-1,str,str,int,int,int,int,str,str)
  47     conv_indel = (str,lambda x: int(x)-1,str,str,int,int,int,int,str,str,int,int,int)
  48
  49     for line in infile:
  50         d = line[:-1].split()
  51         if d[2] == "*":
  52             try:
  53                 yield PileupIndel( *[x(y) for x,y in zip(conv_indel,d) ] )
  54             except TypeError:
  55                 raise pysam.SamtoolsError( "parsing error in line: `%s`" % line)
  56         else:
  57             try:
  58                 yield PileupSubstitution( *[x(y) for x,y in zip(conv_subst,d) ] )
  59             except TypeError:
  60                 raise pysam.SamtoolsError( "parsing error in line: `%s`" % line)
  61
  62 ENCODE_GENOTYPE = {
  63     'A': 'A', 'C': 'C', 'G': 'G', 'T': 'T',
  64     'AA': 'A', 'CC': 'C', 'GG': 'G', 'TT': 'T', 'UU': 'U',
  65     'AG': 'r', 'GA': 'R',
  66     'CT': 'y', 'TC': 'Y',
  67     'AC': 'm', 'CA': 'M',
  68     'GT': 'k', 'TG': 'K',
  69     'CG': 's', 'GC': 'S',
  70     'AT': 'w', 'TA': 'W',
  71     }
  72
  73 DECODE_GENOTYPE = {
  74     'A': 'AA',
  75     'C': 'CC',
  76     'G': 'GG',
  77     'T': 'TT',
  78     'r': 'AG', 'R': 'AG',
  79     'y': 'CT', 'Y': 'CT',
  80     'm': 'AC', 'M': 'AC',
  81     'k': 'GT', 'K': 'GT',
  82     's': 'CG', 'S': 'CG',
  83     'w': 'AT', 'W': 'AT',
  84     }
  85
  86 ##------------------------------------------------------------
  87 def encodeGenotype( code ):
  88     '''encode genotypes like GG, GA into a one-letter code.
  89     The returned code is lower case if code[0] < code[1], otherwise
  90     it is uppercase.
  91     '''
  92     return ENCODE_GENOTYPE[ code.upper() ]
  93
  94 def decodeGenotype( code ):
  95     '''decode single letter genotypes like m, M into two letters.
  96     This is the reverse operation to :meth:`encodeGenotype`.
  97     '''
  98     return DECODE_GENOTYPE[ code ]
  99
 100 def translateIndelGenotypeFromVCF( vcf_genotypes, ref ):
 101     '''translate indel from vcf to pileup format.'''
 102
 103     # indels
 104     def getPrefix( s1, s2 ):
 105         '''get common prefix of strings s1 and s2.'''
 106         n = min( len( s1), len( s2 ) )
 107         for x in range( n ):
 108             if s1[x] != s2[x]: return s1[:x]
 109         return s1[:n]
 110
 111     def getSuffix( s1, s2 ):
 112         '''get common sufix of strings s1 and s2.'''
 113         n = min( len( s1), len( s2 ) )
 114         if s1[-1] != s2[-1]: return ""
 115         for x in range( -2, -n - 1, -1 ):
 116             if s1[x] != s2[x]: return s1[x+1:]
 117         return s1[-n:]
 118
 119     def getGenotype( variant, ref ):
 120
 121         if variant == ref: return "*", 0
 122
 123         if len(ref) > len(variant):
 124             # is a deletion
 125             if ref.startswith(variant):
 126                 return "-%s" % ref[len(variant):], len(variant) - 1
 127             elif ref.endswith( variant ):
 128                 return "-%s" % ref[:-len(variant)], -1
 129             else:
 130                 prefix = getPrefix( ref, variant )
 131                 suffix = getSuffix( ref, variant )
 132                 shared = len(prefix) + len(suffix) - len(variant)
 133                 # print "-", prefix, suffix, ref, variant, shared, len(prefix), len(suffix), len(ref)
 134                 if shared < 0:
 135                     raise ValueError()
 136                 return "-%s" % ref[len(prefix):-(len(suffix)-shared)], len(prefix) - 1
 137
 138         elif len(ref) < len(variant):
 139             # is an insertion
 140             if variant.startswith(ref):
 141                 return "+%s" % variant[len(ref):], len(ref) - 1
 142             elif variant.endswith(ref):
 143                 return "+%s" % variant[:len(ref)], 0
 144             else:
 145                 prefix = getPrefix( ref, variant )
 146                 suffix = getSuffix( ref, variant )
 147                 shared = len(prefix) + len(suffix) - len(ref)
 148                 if shared < 0:
 149                     raise ValueError()
 150
 151                 return "+%s" % variant[len(prefix):-(len(suffix)-shared)], len(prefix)
 152         else:
 153             assert 0, "snp?"
 154
 155     # in pileup, the position refers to the base
 156     # after the coordinate, hence subtract 1
 157             #pos -= 1
 158
 159     genotypes, offsets = [], []
 160     is_error = True
 161
 162     for variant in vcf_genotypes:
 163         try:
 164             g, offset = getGenotype( variant, ref )
 165         except ValueError:
 166             break
 167
 168         genotypes.append( g )
 169         if g != "*":  offsets.append( offset )
 170
 171     else:
 172         is_error = False
 173
 174     if is_error:
 175         raise ValueError()
 176
 177     assert len(set(offsets )) == 1, "multiple offsets for indel"
 178     offset = offsets[0]
 179
 180     genotypes = "/".join( genotypes )
 181     return genotypes, offset
 182
 183 def vcf2pileup( vcf, sample ):
 184     '''convert vcf record to pileup record.'''
 185
 186     chromosome = vcf.contig
 187     pos = vcf.pos
 188     reference = vcf.ref
 189     allelles = [reference] + vcf.alt
 190
 191     data = vcf[sample]
 192
 193     # get genotype
 194     genotypes = data["GT"]
 195     if len(genotypes) > 1:
 196         raise ValueError( "only single genotype per position, %s" % (str(vcf)))
 197
 198     genotypes = genotypes[0]
 199
 200     # not a variant
 201     if genotypes[0] == ".": return None
 202
 203     genotypes = [ allelles[int(x)] for x in genotypes if x != "/" ]
 204
 205     # snp_quality is "genotype quality"
 206     snp_quality = consensus_quality = data.get( "GQ", [0])[0]
 207     mapping_quality = vcf.info.get( "MQ", [0])[0]
 208     coverage = data.get( "DP", 0)
 209
 210     if len(reference) > 1 or max([len(x) for x in vcf.alt] ) > 1:
 211         # indel
 212         genotype, offset = translateIndelGenotypeFromVCF( genotypes, reference )
 213
 214         return PileupIndel( chromosome,
 215                             pos + offset,
 216                             "*",
 217                             genotype,
 218                             consensus_quality,
 219                             snp_quality,
 220                             mapping_quality,
 221                             coverage,
 222                             genotype,
 223                             "<" * len(genotype),
 224                             0,
 225                             0,
 226                             0 )
 227
 228     else:
 229
 230         genotype = encodeGenotype( "".join(genotypes) )
 231
 232
 233         read_bases = ""
 234         base_qualities = ""
 235
 236         return PileupSubstitution( chromosome, pos, reference,
 237                                    genotype,
 238                                    consensus_quality,
 239                                    snp_quality,
 240                                    mapping_quality,
 241                                    coverage, read_bases, base_qualities )
 242
 243
 244 def iterate_from_vcf( infile, sample ):
 245     '''iterate over a vcf-formatted file.
 246
 247     *infile* can be any iterator over a lines.
 248
 249     The function yields named tuples of the type :class:`pysam.Pileup.PileupSubstitution`
 250     or :class:`pysam.Pileup.PileupIndel`.
 251
 252     Positions without a snp will be skipped.
 253
 254     This method is wasteful and written to support same
 255     legacy code that expects samtools pileup output.
 256
 257     Better use the vcf parser directly.
 258
 259     '''
 260
 261
 262     vcf = pysam.VCF()
 263     vcf.connect( infile )
 264
 265     if sample not in vcf.getsamples():
 266         raise KeyErorr( "sample %s not vcf file" )
 267
 268     for row in vcf.fetch():
 269         result = vcf2pileup( row, sample )
 270         if result: yield result
 271