Adding indel support to buildPhylipLineage

[Eric Mukherjee]

It'd be super helpful to add indel support to the buildPhylipLineage function (maybe using the SEQUENCE_INPUT column?) Our data has a ton of indels (both real and erroneous) and we're really interested in seeing how indels increase diversity in our lineages of interest; would be happy to talk more about this by email or whatever.

[Jason Vander Heiden]

Hey @emukherj,

I added a mask_char option to makeChangeoClone and a dist_mat option to buildPhylipLineage, which will, hopefully, let you count indels as mismatches.

For example:

db <- subset(ExampleDb, CLONE == 3138)
clone <- makeChangeoClone(db, 
    text_fields=c("SAMPLE", "ISOTYPE"), 
    num_fields="DUPCOUNT", mask_char="-")

dnapars_exec <- "~/apps/phylip-3.69/dnapars"
graph <- buildPhylipLineage(clone, dnapars_exec, 
    dist_mat=getDNAMatrix(gap=-1), rm_temp=TRUE)

This:

1. Replaces all dot/dash character with dashes (so keeping alignment gaps and converting IMGT gaps to dashes).
2. Treats indels (runs of different length dashes) of any length as a single mismatch (gap=-1). If you instead set gap=1, then this will count each individual gap character as a mismatch (if it's not a gap in both sequences).

I didn't test extensively, but let me know if this accomplishes what you wanted.

[Jason Vander Heiden]

Oh, I forgot to mention. I recently revived the AlignRecords script in changeo (the version on bitbucket), which will allow you to perform a multiple alignment between the germline column and SEQUENCE_VDJ or SEQUENCE_INPUT columns, which should let you preserve indels for lineage building:

http://changeo.readthedocs.io/en/latest/tools/AlignRecords.html

[Eric Mukherjee]

Hey Jason, thanks for all your help. I tried AlignRecords just now but it didn't work; can I send you my data and the error message to see if its something you can fix?

Thanks!

[Jason Vander Heiden]

Sure, please email me the data, error message and command used and I'll see if I can figure out what's wrong.

[Jason Vander Heiden]

No news is good news?