Capsid alignment and variants

[schorlton-bugseq]

Hi, thanks for the great tool! I am interested in capsid variants causing lenacapavir resistance, as HIVDB now includes those in its XML.

I am happy to contribute this addition via PR but don't even know where to start after browsing the code. I tried adding it to

sierra-local/sierralocal/nucaminohook.py

Line 388 in 2896b83

min_overlap = {'PR': 40, 'RT': 60, 'IN': 30}

but no success on reporting. Any pointers where modifications are needed for capsid variants/LEN resistance would be appreciated.

Thanks!

[ArtPoon]

Hi @schorlton-bugseq, thanks for the offer to help! If HIVdb has updated its XML to include resistance mutations in gag then I would start with checking how those parts of the XML are being parsed correctly. Next, there are some hard-coded components that are specific to pol that would need to be updated, e.g.:

sierra-local/sierralocal/nucaminohook.py

Lines 69 to 80 in 2896b83

	self.pi_dict = self.prevalence_parser('PIPrevalences.tsv')
	self.rti_dict = self.prevalence_parser('RTIPrevalences.tsv')
	self.ini_dict = self.prevalence_parser('INIPrevalences.tsv')
	# initialize gene map
	self.pol_start = 2085
	self.pol_nuc_map = {
	'PR': (2253, 2549),
	'RT': (2550, 4229), # incorrectly includes RNAse, emulating sierrapy
	'IN': (4230, 5096)
	}
	self.gene_map = self.create_gene_map()

It's been a long time since I've looked at this code - @WilliamZekaiWang might be more help.

[WilliamZekaiWang]

moving comments from #109 to here so we have a concise log ab CA alignments

adding "CA" to the self.pol_nuc_map puts it into here to reindex it in respects to pol

sierra-local/sierralocal/nucaminohook.py

Lines 360 to 371 in 2896b83

	def create_gene_map(self):
	"""
	Returns a dictionary with the AMINO ACID position bounds
	for each gene in Pol, based on the HXB2 reference annotations.
	@return pol_aa_map: dict
	"""
	# start and end nucleotide coordinates in HXB2 pol
	convert = lambda x: int((x - self.pol_start) / 3)
	pol_aa_map = {}
	for key, val in self.pol_nuc_map.items():
	pol_aa_map[key] = (convert(val[0]), convert(val[1]))
	return pol_aa_map

The result of this indexes CA in negative indexes, which then gets filtered out here:

sierra-local/sierralocal/nucaminohook.py

Lines 373 to 413 in 2896b83

	def get_genes(self, pol_aligned_sites, pol_first_aa, pol_last_aa):
	"""
	Determines the first POL gene that is present in
	the query sequence, by virtue of gene breakpoints
	TODO: sierra uses different minimum numbers of sites per gene
	(40, 60 and 30 for PR, RT and IN)
	@param pol_aligned_sites: list, sublist holds alignment program
	output aligned POL sites
	@param pol_first_aa: int, location of first amino acid in pol
	@param pol_last_aa: int, location of last amino acid in pol
	@return: list, list of length 1
	[gene, first aa position in pol, last aa position in pol
	first na position in pol, last na position in pol]
	"""
	# good here
	min_overlap = {'PR': 40, 'RT': 60, 'IN': 30}
	genes = []
	for gene, bounds in self.gene_map.items():
	aa_start, aa_end = bounds
	geneLength = aa_end - aa_start + 1
	try:
	overlap = min(aa_end, pol_last_aa) - max(aa_start, pol_first_aa)
	if overlap < min_overlap[gene]:
	# discard alignment of this gene, too short
	continue
	aligned_sites = filter(lambda x: x['PosAA'] >= aa_start and x['PosAA'] <= aa_end,
	pol_aligned_sites)
	aligned_sites = list(aligned_sites)
	first_aa = max(pol_first_aa - aa_start, 1)
	last_aa = min(pol_last_aa - aa_start + 1, geneLength)
	first_na = aligned_sites[0]['PosNA']
	last_na = aligned_sites[-1]['PosNA'] - 1 + aligned_sites[-1]['LengthNA']
	genes.append((gene, first_aa, last_aa, first_na, last_na))
	except:
	# if post-align returns nothing, need to raise message
	# "There were no PR, RT, and IN genes found, refuse to process."
	pass
	return genes