1
$\begingroup$

I was looking at this script from https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297 and wanted to use it to determine outliers for tentative Horizontal gene transfer.

The data set is provided in the example, under the section 20.2.2 Plot of sequence GC%.

I want to yield a identifier for the x axis instead of just an arbitrary number of genes or cds' so I can actually determine which sequences have the outlying values, as I can not with the current graph.

image taken from biopython

I am not great at python or understanding list comprehensions, so if anyone could augment the code provided in the example to do so, it would be appreciated.

I provided the code below also.

#from : https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297, data set under 
#section 20.2.2 Plot of sequence GC%, ls_orchid.fasta

from Bio import SeqIO
from Bio.SeqUtils import GC

gc_values = sorted(GC(rec.seq) for rec in SeqIO.parse("ls_orchid.fasta", "fasta"))
import pylab
pylab.plot(gc_values)
pylab.title("%i orchid sequences\nGC%% %0.1f to %0.1f" \
            % (len(gc_values),min(gc_values),max(gc_values)))
pylab.xlabel("Genes")
pylab.ylabel("GC%")
pylab.show()
Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ If you goal is to guess putative HGT, is there a reason why not to use one of the pipelines built for it? $\endgroup$
    – Kamil S Jaron
    Apr 22, 2018 at 22:01
  • $\begingroup$ I was not aware of the pipelines. Or i guess i didnt feel comfortable using them. @KamilSJaron did you have any suggestions? $\endgroup$
    – Ro Siv
    Apr 23, 2018 at 11:43

2 Answers 2

Reset to default
2
$\begingroup$

The tricky bit here is that you need a list of sequence records sorted in the same order as the CG values. You can accomplish this by making tuples of the seq records and GC values, then sorting the tuples specifying the GC values as the sorting key:

gc_w_records = sorted(((GC(rec.seq), rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")), key=lambda tup: tup[0])
gc_values = [x[0] for x in gc_w_records]
records = [x[1] for x in gc_w_records]

I don't know much about matplotlib and it's going to be challenging to display 100 IDs, but you'll probably be interested in the extremes of the GC distribution, which you can obtain directly from the sorted list of records:

records[0:2]
# [SeqRecord(seq=Seq('CGTCACGAGGTCTCCGGATGTGACCCTGCGGAAGGATCATTGTTGAGATCACAT...CAT', SingleLetterAlphabet()), id='gi|2765587|emb|Z78462.1|PSZ78462', name='gi|2765587|emb|Z78462.1|PSZ78462', description='gi|2765587|emb|Z78462.1|PSZ78462 P.sukhakulii 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...AGG', SingleLetterAlphabet()), id='gi|2765568|emb|Z78443.1|PLZ78443', name='gi|2765568|emb|Z78443.1|PLZ78443', description='gi|2765568|emb|Z78443.1|PLZ78443 P.lawrenceanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

or:

records[-3:-1]
# [SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet()), id='gi|2765656|emb|Z78531.1|CFZ78531', name='gi|2765656|emb|Z78531.1|CFZ78531', description='gi|2765656|emb|Z78531.1|CFZ78531 C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGCTCATTGTTGAGACCGCAA...AAG', SingleLetterAlphabet()), id='gi|2765625|emb|Z78500.1|PWZ78500', name='gi|2765625|emb|Z78500.1|PWZ78500', description='gi|2765625|emb|Z78500.1|PWZ78500 P.warszewiczianum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]
Improve this answer
$\endgroup$
3
  • $\begingroup$ +1 but would be marginally faster to use operator.itemgetter(0) as the key for sorted, or just omit the key and the sorting will occur by first the GC then the record $\endgroup$
    – Chris_Rands
    Apr 20, 2018 at 9:20
  • $\begingroup$ I did try that, but you can't sort sequence records: "NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest." $\endgroup$
    – heathobrien
    Apr 20, 2018 at 9:31
  • $\begingroup$ Ah yes good point, forgot about that, itemgetter(0) would be the way to go then, although by not essential $\endgroup$
    – Chris_Rands
    Apr 20, 2018 at 9:58
2
$\begingroup$

There are various attempts to guess putative horizontal gene transfer using sequence composition (kmer-frequencies or GC content), you do not have to implement a new method by your own.

I co-authored one of them called SigHunt. It's kind of working, but the code is really messy and there are now bit more precise methods out there. The advantage is SigHunt is that it provides access to all the data calculated on the way (the tetranucleotide frequencies), so if you want to play with it, it's fun.

The method I would probably go for now if I would have searched for HGT is MTGIpick. The most of the tools detects HGT within a genome by a sliding window, but MTGIpick is refining boundaries where the change of kmer frequencies is rapid. I never tried it personally but the approach is sound to me.

Alternatively there are other paradigms for detection of HGT. Phylogeny based methods (there are many, but for instance) or by mapping reads to species boundaries. Phylogeny based methods usually provide stronger evident, but have a disadvantage to be database dependent (i.e. you can not build a tree if you don't have sufficient coverage of homologs across tree of life).

I would also like to mention that right now I do not follow the HGT literature closely, there might be some big new I just missed.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.