Add step to calculate self-similarity of MT epitope

[susannasiebert]

Original issue: https://github.com/griffithlab/pVAC-Seq/issues/222

Take the peptide sequence and do some sort of similarity search against the ref proteome (BLAST) to see if the peptide naturally exists in the body. If so, it wouldn't be recognized as "foreign" antigen.

1. Implementation: Blast against the reference proteome
2. Implementation: Blast against the personal proteome

[malachig]

Some relevant concepts here: https://bmccancer.biomedcentral.com/articles/10.1186/s12885-018-4325-6

[susannasiebert]

@malachig @jhundal @chrisamiller I have implemented the framework that would use blastp to align the epitopes to the reference proteome and returns the number of alignments found. A couple of follow up things to discuss:

• blastp appears to return alignments for each isoform. Do we want to count different isoform of the same gene? If so, do we want to return identifiers for the matches found to make it easier to identify which alignment were made to the same position in the reference?
• Do we want to hard-filter epitopes that have matches to the reference proteome?
• Do we want to output the match count in the condensed file?
• Do we want to include the match count in the score/rank?

[malachig]

Good questions.

Initial thoughts:

• Probably we just want to count one hit per gene. Otherwise the count returned will largely be a reflection of how many alternative isoforms involve a particular exon rather that how many places in the genome generate a match.

• It would be good to have filtering on these matches be an option, perhaps with the default value to remove candidates that have >=1 match in the wild type proteome.

• If we are generally filtering these out, then including in the condensed file is probably not needed, since the value will always be 0 there. On the other hand, if it is an option to allow them, then in that scenario we would want to see which candidates had matches and which did not. For now, I think we should add it to the condensed file. This will help us confirm that it is working as expected because at least for a while the clinical teams will continue to perform their own manual analysis for these matches. If there are discrepancies between what they do and our implementation here, this will help shine light on them.

• I think we should gain some experience with the output first and decide on whether to include it in the ranking calculation later.

[susannasiebert]

Decision from today's meeting:

• Output a combined fasta with the headers matching the index into the output directory
• Output the index index for each neoeptiope's variant to the all_epitopes and filtered reports
• For each neoepitope, blast the mutant peptide sequence of the variant from the above fasta to the ref_sequence
• If there are any alignments that match at least 8 amino acids, count it as a match
  • If the alignment is a gapped alignment, compare the substring matches
  • Attempt to set a gap penalty that avoids gapped alignments
• Output whether matches were found Y/N in both the filtered and condensed result files
• Don't hard filter
• Stretch goal: implement on optional filter.