Incubator for useful bioinformatics code, primarily in Python and R

Even putting numpy >=1.6.1 in setup.py's before biopython, the following message pops up:

Numerical Python (NumPy) is not installed.

This package is required for many Biopython features.  Please install
it before you install Biopython. You can install Biopython anyway, but
anything dependent on NumPy will not work. If you do this, and later
install NumPy, you should then re-install Biopython.

You can find NumPy at http://numpy.scipy.org

Do you want to continue this installation? (y/N):

Apparently install_requires packages are not installed in order, so no dependency order can be defined that way... are you aware of any "pre_install_requires" or similar in setuptools ? Couldn't find it after quickly checking docs :-/

i'm not sure if it's a problem in the latest version of barcode_sort_trim.py

After updating to the pipeline with FastQC, I've got an extra base 'A' in the 3' of read 1.

Brad, It's not really an issue. But I want to know, from your experience, how much time you would save from switching to FastQC from SolexaQA?

Thanks, Paul

Hi,

We have a different setting here where the drive with the bcl files is mounted to the analysis machine and we would do everything there. Do you recommend we keep the messaging system in the pipeline? Just want to get some advices.

Thanks, Paul

Hi Brad,

it seems that it will keep finding CreateSequenceDictionary in /usr/share/java/picard even though I have specify another path in my config file? I have tried doing the setup again after I modified the config files, but still it didn't look up the path I specified.

and I didn't seem to have specified the path of hg19.fa for GATK?

Thanks, Paul

Submitting to GenBank requires converting a GFF file to a Sequin TBL file, which is then converted ASN.1 using tbl2asn. I have searched, and I have not found a good (or any, really) converter from GFF to Sequin TBL. Would you be interested in adding such a tool? Here's the hacky script that I cobbled together for this purpose: gff3-to-tbl. It's not general purpose, but could be a useful starting point.

Hi Brad,

more of a question than an issue. I noticed you've added code (bcbio.pipeline.sample.merge_sample) to merge samples across lanes. I've been using save_diskspace=true in order to remove sam files, but this I noticed also removes the demultiplexed files, right? I just want to make sure because it affects our data delivery routines, as outlined below.

In our setup, we have situations when we run several projects on one lane, which we distinguish with an extra "description" tag in run_info, so in principle each barcode could have a description with a different project name. We then partition fastq files in a lane based on the description tag when delivering data to customers.

On a similar note, when I do analyses for customers, I've been doing it on a project-by-project basis (it makes more sense to me), and therefore written helper scripts (project_*, see EDIT: https://github.com/percyfal/bcbb/tree/develop/nextgen/scripts) for this purpose. project_analysis_pipeline.sh is almost a copy of automated_initial_analysis.py, but starts off with demultiplexed files. Have you had this functionality in mind (or is it even already there)?

Cheers,

Per

Hi Brad,

We are seeing some issues with unexpectedly many reads ending up in the 'unmatched' category after demultiplexing. After digging around a little, we think that this may be related to the trailing 'A' that the Illumina machines add after the barcode.

More specifically, we allow one mismatch and no indels for the demuxing. It seems that the reads that are unexpectedly classified as unmatched have one mismatch in the actual 6-nucleotide barcode and are, in addition, having the trailing 'A' nucleotide miscalled.

Reading the code, it does indeed seem that for Illumina reads, the last 7 nucleotides, including the trailing 'A', of each read are matched when demultiplexing. Can you confirm that this is the case?

Our preference is to match just the 6-mer index sequence, excluding the last nucleotide in the read and it would be nice to have this done by default for Illumina reads, or at least be able to influence this behavior with a configuration option. What do you think?

Thanks /Pontus

I tried following http://biopython.org/wiki/GFF_Parsing to parse UCSC-generated GTF file.

After executing

pprint.pprint(examiner.parent_child_map(handle))

the output was

{}

Similarly,

examiner.available_limits(handle)

produced

3: {'gff_id': {}, 'gff_source': {}, 'gff_source_type': {}, 'gff_type': {}}

Trying to parse that same file with

from BCBio import GFF
for rec in GFF.parse(handle):
    print rec

produced

ID: chr1
Name: <unknown name>
Description: <unknown description>
Number of features: 2
UnknownSeq(14409, alphabet = Alphabet(), character = '?')

Here are the first 10 lines from the GTF in question

chr1 hg19_knownGene exon 11874 12227 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3"; chr1 hg19_knownGene exon 12613 12721 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3"; chr1 hg19_knownGene exon 13221 14409 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3"; chr1 hg19_knownGene start_codon 12190 12192 0.000000 + . gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; chr1 hg19_knownGene CDS 12190 12227 0.000000 + 0 gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; chr1 hg19_knownGene exon 11874 12227 0.000000 + . gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; chr1 hg19_knownGene CDS 12595 12721 0.000000 + 1 gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; chr1 hg19_knownGene exon 12595 12721 0.000000 + . gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; chr1 hg19_knownGene CDS 13403 13636 0.000000 + 0 gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; chr1 hg19_knownGene stop_codon 13637 13639 0.000000 + . gene_id "uc010nxq.1"; transcript_id "uc010nxq.1";

Hey Brad, we are trying to use the distributed version of the pipeline.

We have a couple of test sets that we use to quickly see if the pipeline is working. One that takes the normal pipeline about 3 hours to finish, and another much smaller that takes about 7 minutes (this is with 8 cores).

When running the small test set on the messaging variant all files get generated as they should, and the program exits properly. Note that this small set consists of fastq files which are only 12 lines each, and I'm guessing much of the analysis gets skipped due to a lack of data.

When we run the messaging version of the pipeline for the larger set, the programs work for a while (time varies, but say between 45 minutes and 1 hour 30 minutes), but then one of the jobs crashes with a socket.timeout error, (this specific job I believe is some master that coordinates what the other jobs should be doing.

I'll include the output of that job here:

[2012-02-25 02:55:26,856] Found YAML samplesheet, using /proj/a2010002/nobackup/illumina/pipeline_test/archive/000101_SN001_001_AABCD99XX/run_info.yaml instead of Galaxy API
Traceback (most recent call last):
  File "/bubo/home/h10/vale/.virtualenvs/devel/bin/automated_initial_analysis.py", line 7, in <module>
    execfile(__file__)
  File "/bubo/home/h10/vale/bcbb/nextgen/scripts/automated_initial_analysis.py", line 117, in <module>
    main(*args, **kwargs)
  File "/bubo/home/h10/vale/bcbb/nextgen/scripts/automated_initial_analysis.py", line 48, in main
    run_main(config, config_file, fc_dir, work_dir, run_info_yaml)
  File "/bubo/home/h10/vale/bcbb/nextgen/scripts/automated_initial_analysis.py", line 65, in run_main
    lane_items = run_parallel("process_lane", lanes)
  File "/bubo/home/h10/vale/bcbb/nextgen/bcbio/distributed/messaging.py", line 28, in run_parallel
    return runner_fn(fn_name, items)
  File "/bubo/home/h10/vale/bcbb/nextgen/bcbio/distributed/messaging.py", line 67, in _run
    while not result.ready():
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/result.py", line 306, in ready
    return all(result.ready() for result in self.results)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/result.py", line 306, in <genexpr>
    return all(result.ready() for result in self.results)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/result.py", line 108, in ready
    return self.status in self.backend.READY_STATES
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/result.py", line 196, in status
    return self.state
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/result.py", line 191, in state
    return self.backend.get_status(self.task_id)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/backends/base.py", line 237, in get_status
    return self.get_task_meta(task_id)["status"]
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/backends/amqp.py", line 128, in get_task_meta
    return self.poll(task_id)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/celery/backends/amqp.py", line 153, in poll
    with self.app.pool.acquire_channel(block=True) as (_, channel):
  File "/sw/comp/python/2.7.1_kalkyl/lib/python2.7/contextlib.py", line 17, in __enter__
    return self.gen.next()
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/kombu/connection.py", line 789, in acquire_channel
    yield connection, connection.default_channel
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/kombu/connection.py", line 593, in default_channel
    self.connection
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/kombu/connection.py", line 586, in connection
    self._connection = self._establish_connection()
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/kombu/connection.py", line 546, in _establish_connection
    conn = self.transport.establish_connection()
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/kombu/transport/amqplib.py", line 252, in establish_connection
    connect_timeout=conninfo.connect_timeout)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/kombu/transport/amqplib.py", line 62, in __init__
    super(Connection, self).__init__(*args, **kwargs)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/amqplib/client_0_8/connection.py", line 129, in __init__
    self.transport = create_transport(host, connect_timeout, ssl)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/amqplib/client_0_8/transport.py", line 281, in create_transport
    return TCPTransport(host, connect_timeout)
  File "/bubo/home/h10/vale/.virtualenvs/devel/lib/python2.7/site-packages/amqplib/client_0_8/transport.py", line 85, in __init__
    raise socket.error, msg
socket.timeout: timed out
[INFO/MainProcess] process shutting down
[DEBUG/MainProcess] running all "atexit" finalizers with priority >= 0
[DEBUG/MainProcess] running the remaining "atexit" finalizers

Have you encountered any issues with socket.timeout? Any ideas what we might be doing wrong?

Overview

After upgrading to Biopython 1.68, GFF.parse() is now failing where it had no issues before.

To Reproduce

In a new virtualenv environment, run:

pip install numpy
pip install biopython
pip install bcbio-gff

wget http://tritrypdb.org/common/downloads/release-27/TcruziCLBrenerEsmeraldo-like/gff/data/TriTrypDB-27_TcruziCLBrenerEsmeraldo-like.gff

Next, launch python and run:

>>> from BCBio import GFF
>>> gff = 'TriTrypDB-27_TcruziCLBrenerEsmeraldo-like.gff'
>>> x=list(GFF.parse(gff))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/keith/.virtualenvs/gff/lib/python3.5/site-packages/BCBio/GFF/GFFParser.py", line 737, in parse
    target_lines):
  File "/home/keith/.virtualenvs/gff/lib/python3.5/site-packages/BCBio/GFF/GFFParser.py", line 327, in parse_in_parts
    cur_dict = self._results_to_features(cur_dict, results)
  File "/home/keith/.virtualenvs/gff/lib/python3.5/site-packages/BCBio/GFF/GFFParser.py", line 367, in _results_to_features
    results.get('child', []))
  File "/home/keith/.virtualenvs/gff/lib/python3.5/site-packages/BCBio/GFF/GFFParser.py", line 428, in _add_parent_child_features
    children)
  File "/home/keith/.virtualenvs/gff/lib/python3.5/site-packages/BCBio/GFF/GFFParser.py", line 471, in _add_children_to_parent
    cur_child, _ = self._add_children_to_parent(cur_child, children)
  File "/home/keith/.virtualenvs/gff/lib/python3.5/site-packages/BCBio/GFF/GFFParser.py", line 477, in _add_children_to_parent
    cur_parent.sub_features.append(cur_child)
AttributeError: 'SeqFeature' object has no attribute 'sub_features'
>>> import Bio
>>> Bio.__version__
'1.68'

The same code worked with Biopython 1.67, so it seems likely to be an issue resulting from changes made in the 1.68 release.

There are small typos in:

• posts/conferences/bosc2018_day1a.md
• posts/seminars/tumor_heterogeneity_carter.md

Fixes:

• Should read suppressors rather than supressors.
• Should read service rather than serivce.

Semi-automated pull request generated by https://github.com/timgates42/meticulous/blob/master/docs/NOTE.md

Hi,

I'd like to add this issue for those who'd like to use the script with an alternative codon table.

Example: If you want to translate the sequence with a codon table 6 you need to change the script like the following, protein_seq = gene_seq.translate(6)

Regards, Zeynep

Hi,

Thank you for this very useful script!!

I was wondering if it's possible to create the protein multi fasta file with the order of ref_file (contigs) instead of the glimmer_file (gff) file?

I've re-ordered my assembly from large contigs to smaller ones. However, when I run this script I got my protein multi fasta file in the order of glimmer.gff instead of the contigs.

Regards, Zeynep

I was just wondering if it would be possible for a new release of the gff package to be made sometime soon. The fix from #126 would be really nice to have in a released version.

Hi!

I annotated a bacterium (Acidipropionibacterium acidipropionici - strain FAM19036) with NCBI PGAP.

I wanted to create SeqIO-objects from the gff file, but it failed:

import pprint
from BCBio import GFF
from BCBio.GFF import GFFExaminer
examiner = GFFExaminer()
with open('data/FAM19036/annot.gff') as in_handle:
    pprint.pprint(examiner.available_limits(in_handle))
print("------------------------------------------------------------")
with open('FAM19036/annot.gff') as in_handle:
    for rec in GFF.parse(in_handle):
        print(rec)

{'gff_id': {('CP040634.1',): 6772},
 'gff_source': {('.',): 3361,
                ('GeneMarkS-2+',): 360,
                ('Local',): 1,
                ('Protein Homology',): 2916,
                ('cmsearch',): 24,
                ('tRNAscan-SE',): 110},
 'gff_source_type': {('.', 'exon'): 8,
                     ('.', 'gene'): 3208,
                     ('.', 'pseudogene'): 137,
                     ('.', 'rRNA'): 8,
                     ('GeneMarkS-2+', 'CDS'): 360,
                     ('Local', 'region'): 1,
                     ('Protein Homology', 'CDS'): 2916,
                     ('cmsearch', 'RNase_P_RNA'): 1,
                     ('cmsearch', 'SRP_RNA'): 1,
                     ('cmsearch', 'exon'): 7,
                     ('cmsearch', 'rRNA'): 4,
                     ('cmsearch', 'riboswitch'): 10,
                     ('cmsearch', 'tmRNA'): 1,
                     ('tRNAscan-SE', 'exon'): 55,
                     ('tRNAscan-SE', 'tRNA'): 55},
 'gff_type': {('CDS',): 3276,
              ('RNase_P_RNA',): 1,
              ('SRP_RNA',): 1,
              ('exon',): 70,
              ('gene',): 3208,
              ('pseudogene',): 137,
              ('rRNA',): 12,
              ('region',): 1,
              ('riboswitch',): 10,
              ('tRNA',): 55,
              ('tmRNA',): 1}}
------------------------------------------------------------

Error
Traceback (most recent call last):
  File "/usr/lib64/python3.7/unittest/case.py", line 59, in testPartExecutor
    yield
  File "/usr/lib64/python3.7/unittest/case.py", line 628, in run
    testMethod()
  File "/project/gene_loci_comparison/test_gene_loci_comparison.py", line 129, in test_recreate_gff_bug
    for rec in GFF.parse(in_handle):
  File "/project/venvs/gene_loci_comparison/lib64/python3.7/site-packages/BCBio/GFF/GFFParser.py", line 746, in parse
    target_lines):
  File "/project/venvs/gene_loci_comparison/lib64/python3.7/site-packages/BCBio/GFF/GFFParser.py", line 327, in parse_in_parts
    cur_dict = self._results_to_features(cur_dict, results)
  File "/project/venvs/gene_loci_comparison/lib64/python3.7/site-packages/BCBio/GFF/GFFParser.py", line 369, in _results_to_features
    base = self._add_directives(base, results.get('directive', []))
  File "/project/venvs/gene_loci_comparison/lib64/python3.7/site-packages/BCBio/GFF/GFFParser.py", line 388, in _add_directives
    val = (val[0], int(val[1]) - 1, int(val[2]))
IndexError: tuple index out of range

To recreate the bug, here is the relevant gff file.

Thanks in advance.

Edit: bcbio-gff version 0.6.6

I'm trying to parse a gff file downloaded from NCBI (GCA_001536265) and when I iterate on the parser it gives me this error Did not find remapped ID location: gene670, [[42143, 44074], [44736, 45087], [45979, 46332], [47064, 47369]], [42143, 47369]

Inspecting the GFF with GFFExaminer gives no error at all.