<p>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:</p>

<pre><code>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] 
</code></pre>

<p>I don't know much about <code>matplotlib</code> 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:</p>

<pre><code>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=[])]
</code></pre>

<p>or:</p>

<pre><code>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=[])]
</code></pre>

CC BY-SA 3.0
2018-04-20T09:20:49.813Z
+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
Chris_Rands
CC BY-SA 3.0
2018-04-20T09:31:59.723Z
I did try that, but you can't sort sequence records: "NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest."
heathobrien
CC BY-SA 3.0
2018-04-20T09:58:26.620Z
Ah yes good point, forgot about that, `itemgetter(0)` would be the way to go then, although by not essential
Chris_Rands
CC BY-SA 3.0
2018-04-20T08:05:41.267Z
2018-04-20T08:05:41.267Z
heathobrien
<p>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. </p>

<p>I co-authored one of them called <a href="https://academic.oup.com/bioinformatics/article/30/8/1081/254636" rel="nofollow noreferrer">SigHunt</a>. 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.</p>

<p>The method I would probably go for now if I would have searched for HGT is <a href="https://academic.oup.com/bib/advance-article-abstract/doi/10.1093/bib/bbw118/2726041?redirectedFrom=fulltext" rel="nofollow noreferrer">MTGIpick</a>. 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.</p>

<p>Alternatively there are other paradigms for detection of HGT. Phylogeny based methods (there are many, but <a href="https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-11-324" rel="nofollow noreferrer">for instance</a>) or by <a href="https://academic.oup.com/bioinformatics/article/32/17/i595/2450753" rel="nofollow noreferrer">mapping reads to species boundaries</a>. 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).</p>

<p>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.</p>

CC BY-SA 3.0
2018-04-23T23:44:32.183Z
2018-04-30T20:16:44.190Z
2018-04-30T20:16:44.190Z
Kamil S Jaron
<p>I was looking at this script from <a href="https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297" rel="nofollow noreferrer">https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297</a> and wanted to use it to determine outliers for tentative Horizontal gene transfer.</p>
<p>The data set is provided in the example, under the section <strong>20.2.2 Plot of sequence GC%</strong>.</p>
<p>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.</p>
<p><a href="https://i.sstatic.net/TjoD4.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/TjoD4.png" alt="image taken from biopython" /></a></p>
<p>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.</p>
<p>I provided the code below also.</p>
<pre><code>#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(&quot;ls_orchid.fasta&quot;, &quot;fasta&quot;))
import pylab
pylab.plot(gc_values)
pylab.title(&quot;%i orchid sequences\nGC%% %0.1f to %0.1f&quot; \
            % (len(gc_values),min(gc_values),max(gc_values)))
pylab.xlabel(&quot;Genes&quot;)
pylab.ylabel(&quot;GC%&quot;)
pylab.show()
</code></pre>

CC BY-SA 3.0
2018-04-22T22:01:09.417Z
If you goal is to guess putative HGT, is there a reason why not to use one of the pipelines built for it?
Kamil S Jaron
CC BY-SA 3.0
2018-04-23T11:43:45.773Z
I was not aware of the pipelines. Or i guess  i didnt feel comfortable using them. @KamilSJaron      did you have any suggestions?
Ro Siv
<gene><genome><python><identifiers>
Obtaining identifier from plot of sequence GC%
# Obtaining identifier from plot of sequence GC%

(Asked by: username_0 on 2018-04-20)

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.

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()

**Tags:** gene, genome, python, identifiers
> **Comments:**
> **username_1**: If you goal is to guess putative HGT, is there a reason why not to use one of the pipelines built for it?  
> **username_0**: I was not aware of the pipelines. Or i guess  i didnt feel comfortable using them. @KamilSJaron      did you have any suggestions?  

---

## 2 answers

### Answer by username_2 on 2018-04-20. Score: 2 | Accepted answer

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=[])]

> **Comments:**
> **username_3**: +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  
> **username_2**: I did try that, but you can't sort sequence records: "NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest."  
> **username_3**: Ah yes good point, forgot about that, `itemgetter(0)` would be the way to go then, although by not essential  

---

### Answer by username_1 on 2018-04-23. Score: 2

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](https://academic.oup.com/bioinformatics/article/30/8/1081/254636). 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](https://academic.oup.com/bib/advance-article-abstract/doi/10.1093/bib/bbw118/2726041?redirectedFrom=fulltext). 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](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-11-324)) or by [mapping reads to species boundaries](https://academic.oup.com/bioinformatics/article/32/17/i595/2450753). 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.


---

