How to display differences from Variant Calling File?

[bioinfornatics]

Dear,

I would like to use your awesome library with pysam in order to display the difference between two sequences. I see they are GeneBlock for this purpose. So how to create a such instance from a list of variations ?

Thanks

best regards

[bioinfornatics]

Using pysam.VariantFile I am am able to convert variation to a collection of Bio.SeqFeature which are translated (BiopythonTranslator) to a collection of dna_features_viewer.GraphicRecord.GraphicRecord The result is quite good, I stuck currently to display a coord label for each tick of those features on the ruler. Did you know how to do this ?

[veghp]

Hi, thanks for the question & answer, that looks great! I personally didn't know, so if you don't mind sending the script, then we can include it here and as one of the Examples. To answer the original question, if you can produce the variants as Biopython records, then you can use Geneblocks or CUBA compare-two-sequences.

[Zulko]

+1 for @veghp 's answer above, please share the script if you can! Btw @veghp did you guess what famous sequence this is?

display a coord label for each tick of those features on the ruler. Did you know how to do this ?

One way is to use Matplotlib's ax.set_xticks([list of indices]), but this may create collisions if your indices are too close from one another. Another way is to add the index to the feature label, for instance "535\nReplaced G->T". This way Dna Features Viewer should detect close labels and display them under one another. Cosmetic pro tip: you can use ax2.invert_yaxis() on the bottom ax so that the Replaced G->T annotations will appear below the ax rather than above.

[veghp]

No I didn't look at it closely, but I see now that this is Covid-19! Nice.

[bioinfornatics]

Thanks. The script will be released freely soon :-)