text: user2701: 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()
  
  
  user581: 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=[])]
  
  
  user57: 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.
  
original_text: 2701: <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.stack.imgur.com/TjoD4.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/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>
   581: <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>
   57: <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>
  