mjmlab / pyinseq Goto Github PK

filterSortCounts() appears to be sorting nucleotides as strings (i.e., 150000 before 180). This is not a problem when mapping to the gene. This is a noteworthy bug when looking at output files so I'm marking it here. However, I think the whole architecture of the pipeline needs to be redesigned and therefore this should get examined at that point.

Consider travis-sphinx

Hi,

Is it possible to run the pyinseq pipeline without sample list? My TnSeq data is already demultiplexed after the sequencing run so I cannot provide additional barcode sequences.

Bests,
Jan

Need overall strategy for logging results into a text file and into the terminal.

Issue #20 changed to focus on appearance. This issue is for content that needs (better) documentation.

Currently all 4 bp barcodes are being carried through the pipeline. Only barcodes in the sample list should be analyzed. Check: assign.py

Typical use case is that a new set of Illumina reads arrives. They immediately need to be demultiplexed. In some cases samples from a single run are analyzed; in other cases a combination of old and new samples are compared for analysis. The current setup requires manual steps to first demultiplex the new data and to then select which samples for the combined analysis. Goal of this issue is to further modularize the pipeline to fully separate demultiplexing, mapping, and analysis. Although these are somewhat separate currently, the code to run them is still too intermingled. For example, they use one strain list that is not passed from function to function but instead just runs in a long script. Goal here is to have these as completely separate functions in pyinseq.py.

Here is a scaffold I have started playing with on a separate branch called modular (will push soon).

For a basic analysis in which all samples are from the same run, first the sample list and paths to the relevant files will be written to a config file (by the function tentatively called pipeline_config()). This will be copied into the beginning of a log for the run as mentioned previously (https://github.com/mandel01/pyinseq/issues/7).

The samples will then be demultiplexed and the results of that will be written to separate files that accompany each demultiplexed read (and added to the log).

The mapping will be conducted on each demultiplexed file against the genome file specified in the initial configuration file. The results of the mapping will be added as separate files with the mapping outputs (and added to the log).

Finally, analyses will be conducted as specified in the config file.

This framework should form the basis for specifying samples for analysis and for combining samples from different Illumina runs in future workflows.

def pipeline_config():
    """organize list of samples; write dictionary-like file (yaml?) that other functions will read."""
    pass

def pipeline_demultiplex():
    """demultiplex samples"""
    pass

def pipeline_mapping():
    """inseq mapping"""
    pass

def pipeline_analysis():
    """analysis"""
    pass

Would it be worth automating a way to add like a '_1 or _2' (or even the date?) to the end of the experiment directory? Just warn the user that it exists and that the new one was created.

(pyinseq) naseersangwan@Naseers-MBP Bt_5951 % pyinseq -i GTAC-134320-RG_1.HL2C3DSX5_CCGCGGTTAT-ACAGCGCTAG_L002_R1.fastq -s bla.txt -g PROKKA_09032022.gbk -e BLA
2023-03-01 20:02 - INFO - pyinseq - Process command line arguments
2023-03-01 20:02 - INFO - utils - Writing configuration file from provided arguments
2023-03-01 20:02 - INFO - utils - Directory BLA already exists. Snakemake will create files that are missing.
2023-03-01 20:02 - INFO - pyinseq - Workflow selected: pyinseq
2023-03-01 20:02 - INFO - pyinseq - ***** BEGIN PYINSEQ WORKFLOW *****
(SNAKEMAKE INFO)
(SNAKEMAKE INFO)
2023-03-01 20:02 - INFO - demultiplex - Preparing to demultiplex reads
Error in job demultiplex while creating output files results/BLA/raw_data/Sample_1.fastq, results/BLA/Sample_1_trimmed.fastq.
RuleException:
Warning in line 84 of /Users/naseersangwan/Desktop/Conda/envs/pyinseq/lib/python3.6/site-packages/pyinseq/workflows/PyinseqWorkflow/Snakefile:
Unknown argument(s): {'colour': 'green'}
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq/lib/python3.6/site-packages/pyinseq/workflows/PyinseqWorkflow/Snakefile", line 84, in __rule_demultiplex
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq/lib/python3.6/site-packages/pyinseq/demultiplex.py", line 62, in demultiplex_fastq
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq/lib/python3.6/site-packages/pyinseq/logger.py", line 78, in init
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq/lib/python3.6/site-packages/tqdm-4.7.2-py3.6.egg/tqdm/_tqdm.py", line 427, in init
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq/lib/python3.6/concurrent/futures/thread.py", line 56, in run
Exiting because a job execution failed. Look above for error message
Exiting because a job execution failed. Look above for error message
2023-03-01 20:02 - INFO - pyinseq - ***** PYINSEQ WORKFLOW FINISHED *****

Add a couple new options (--checkpoint or --restartdir?) to be able to go into the experiment folder and grab the samples.yml file for paths to demultiplexed data

There are (at least!) two overall issues needed for optimization. First is more efficient file reading and writing. Second is a better overall plan for not going over the same data multiple times (likely involving numpy and/or pandas).

Will return to both of these once the full pipeline is available here and can be run with test data and timed.

Many errors that Pyinseq outputs are specific to the computational workflow. To aid in debugging and development, a separate python class should be written that will inherit from Exception.

The -d parameter filters by the location of the insertion within the gene.

Here is the current relevant code in process_mapping.py

if strand == "+":
    three_primeness = (nucleotide - start) / (end - start)
if strand == "-":
    three_primeness = (end - nucleotide) / (end - start)

First problem: denominator should be the length of the gene (i.e., 1 more than the current calculation). Suggested solution:

if strand == "+":
    three_primeness = (nucleotide - start) / length
else:
    three_primeness = (end - nucleotide) / length

Second problem: See example below for a transposon insertion at the 4th bp of a 9 bp gene:

+ strand gene:

5'-ATGTACTAG
      ^

^ first nt in the TA site

Current math = ((3 + start) - start) / 9 = 3/9. However, this should be 4/9, so need to add 1.

What about - strand genes? Note this is the same sequence just flipped to the other strand.

- strand gene:
5'-CTAGTACAT    (gene starts over here)
       ^

Current math = ((end - (end - 4) / 9) = 4/9. If we look at the sample gene and count from the 5' end (i.e., the reverse complement of the ATG on the right side), then it should be 5/9. However, for the sake of biology, we want insertions in the same position to be called the same way whether the gene is on the + or - strand. Therefore, I think that 4/9 is the right call here.

Proposal is for the insertion nucleotide to continue to represent the first position of a TA site (i.e., lowest genome coordinate) regardless of what direction the gene(s) are.

5'-TA
3'-AT
   ^
Always call the "first" nucleotide regardless of gene annotations.

However, when we are processing the site data to collect gene-level data, we will use the nt in the TA site that is closest to the 5' end of the gene.

+ strand gene:
5'-ATGTACTAG
      ^
      *

- strand gene:
5'-CTAGTACAT
       ^*

^ = site-level insertion point
* = gene-level three-primeness call

Therefore... proposed solution:

if strand == "+":
    three_primeness = (1 + nucleotide - start) / length
else:
    three_primeness = (end - nucleotide) / length

Running:
$ pip install requirements.txt
Collecting requirements.txt
Could not find a version that satisfies the requirement requirements.txt (from versions: )
No matching distribution found for requirements.txt

Results in the above error. sudo does not fix the issue either.

While analyzing OSU data, we noticed that this Tn-Seq library contains reads that begin with TA. We believe this presents an issue when pyinseq map the bowtie alignments to the sites sincepyinseq was designed for reads that end with TA sites. Here is an example of how two test reads map to the genome:

1. 5 - AGTTGAACCTCGTCGTA - 3
2. 5 - TCTCTTCAATGAAGTTA - 3

(AGTTGAACCTCGTCG) [TA]   ACTTCATTGAAGAGA
 TCAACTTGGAGCAGC  [AT]  (TGAAGTAACTTCTCT)

• Bowtie reports the insertion start location (0-based) for the reads and pyinseq calculates the coordinate of the T in TA site with:

(+) insertion_NT = insertion_NT + read_length - 1
(-) insertion_NT = insertion_NT + 1

This algorithm holds if reverse alignment is located to the right of TA site or if forward alignment is on the left.

• So for these site mappings:

Read 1 (+) : insertion_NT = 13407 + 17 - 1 = 13423
Read 2 (-): insertion_NT = 13422 + 1 = 13423

• For both reads, TA site index would be 13423, which pyinseq uses to count left and right reads that map around this site.

However, when the reads start TA site, the code used to find this coordinate does not work.

• Example using reverse complement of same reads:

1. 5 - TACGACGAGGTTCAACT - 3
2. 5 - TAACTTCATTGAAGAGA - 3

 AGTTGAACCTCGTCG  [TA]  (ACTTCATTGAAGAGA)
(TCAACTTGGAGCAGC) [AT]   TGAAGTAACTTCTCT

• Read 1 is now aligned to the reverse strand, but its insertion coordinate (which points to the beginning of reading) does not flank the targeted TA site

Read 1 (-) : insertion_NT = 13407 + 1 = 13408
Read 2 (+): insertion_NT = 13422 + 17 - 1 = 13438

• Sites file from a stable pyinseq run on reverse-complement fastq files confirms this because it identifies these two locations as independent sites, but in reality, they are not. In addition, sites have either only left or right flanking reads.

$ cat E001_01_reversed_sites.txt
CP000020	13408	0	11	11	43650.793650793654
CP000020	13438	1	0	1	3968.253968253968

Github lists pyinseq as an HTML project (in searches and in the "bar" at the top of the main page) due to the documentation in the packaged bowtie files. According to this link we should "consider moving the files into one of the paths for vendored code, or use the manual overrides feature to ignore them."

How about changing the packages/ directory to be named /dependencies ? Is that enough?

Would also need to tell bowtie where to find itself in config.py.

It would be ideal to have a user tutorial/development documentation setup so that we can quickly sort out who is working on what, progress, etc. Github also renders IPython notebooks now, so we could use that. Otherwise I suggest using readthedocs.org.

I like how the khmer project is set up.

See: http://khmer.readthedocs.org/en/v1.4.1/

Thoughts?

[This isn't as much as issue as some notes that are relevant for this repo. Will put into another form later!]

First, looking at coverage to guide where tests are needed:

$ coverage run scripts/pyinseq.py -i data/example/example01.fastq -s data/example/example01.txt -g data/example/ES114v2.gb -e example01

...

$ coverage report
Name                        Stmts   Miss  Cover
-----------------------------------------------
scripts/config.py              14      3    79%
scripts/demultiplex.py         98     26    73%
scripts/gbkconvert.py         106     10    91%
scripts/mapReads.py            29      2    93%
scripts/processMapping.py     101      4    96%
scripts/pyinseq.py            129      9    93%
scripts/utils.py               20      5    75%
-----------------------------------------------
TOTAL                         497     59    88%

e.g. File formats (gbk and how it can have multiple contigs; ftt format).

https://pypi.org/project/pyinseq/

Generate unit tests and attempt to catch potential errors within functions. Download test data from GEO?

The summary_gene_table.txt output is not tabulating all of the relevant results from the *genes.txt files. In particular, the totals for the columns should be close to 1E6 (total of the counts per million, minus intergenic regions and regions excluded by -d). However, in some cases the total was much lower, ~15,000. It appears that instead of tabulating results in the summary gene table, only the last entry for each gene from each *genes.txt is being listed. For example:

SPL01_genes.txt:

contig	nucleotide	left_counts	right_counts	total_counts	cpm	three_primeness	locus_tag
CM001400	2279	11	6	17	3.611702084716934	0.996309963099631	VFSR5_0003
CM001400	4041	2	1	3	0.6373591914206354	0.4573643410852713	VFSR5_0004
CM001400	4919	6	2	8	1.6996245104550278	0.8577235772357723	VFSR5_0007
CM001400	4921	100	141	241	51.20118837745771	0.8550135501355014	VFSR5_0007
CM001400	5390	8	4	12	2.5494367656825414	0.21951219512195122	VFSR5_0007
CM001400	5833	28	5	33	7.010951105626989	0.5758928571428571	VFSR5_0008
CM001400	5872	6	5	11	2.336983701875663	0.5178571428571429	VFSR5_0008
CM001400	5907	11	5	16	3.3992490209100557	0.46577380952380953	VFSR5_0008
CM001400	6071	2	3	5	1.0622653190343923	0.22172619047619047	VFSR5_0008
CM001400	6176	6	1	7	1.4871714466481492	0.06547619047619048	VFSR5_0008

summary_gene_table:

Note that the genes with no entries or only one entry are correct. However, VFSR5_0007 should sum the total from the 3 entries (1.70, 51.20, 2.55); instead only 2.55 is reported in the table. Similar for VFSR5_0008.

On my 4th Travis CI build and with my commits last week that should have been squashed first, I have made a mess of the master branch. I'm starting a dev branch for unfinished work and will work on having the master branch be a more linear record from here out.

Only the first barcode was read in.

Specifying the command line argument -d 0.9 should map all hits but then only count those for a particular gene if they fall in the 5'-most 0.9 (i.e., 90%) of the gene. This is not working right now.

First should look at how others do this, but wanted to get these thoughts down as I look through others' Python code.

• Python PEP8 style -- not fanatical, but makes the code more readable and easier for others to contribute and debug.
• Single quote for strings except for docstrings.
• 4 spaces (no tabs) - Python 3 requires consistency here.
• Add documentation for a feature as you commit the feature.

In the first few runs have not seen a gradual distribution of normalized counts (CPM), instead have samples binning near 10, 100, or 1000. Is this a float() error perhaps? Better tests (Issue https://github.com/mandel01/pyinseq/issues/3) should help address this.

Reading fastq files into pyinseq is one of the most slow and memory-heavy processes. However, using yield keyword allows us to create a generator that provides a line by line loop without storing it into virtual memory. This is also related to using Tqdm progress bar since we need a quick way to count lines for setting progression timeline.

The example01_bowtieOutput.txt file contains data like this.

$ head example01/temp/example01_bowtieOutput.txt 
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   +   CP000020    11136   CGACACCGACGACGGTA   FFDDHHHHHCGHHHIJD   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   -   CP000020    11151   TACAGGTCTTCACCACA   DJIHHHGCHHHHHDDFF   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   -   CP000020    13422   TAACTTCATTGAAGAGA   DJIHHHGCHHHHHDDFF   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   +   CP000020    13407   AGTTGAACCTCGTCGTA   FFDDHHHHHCGHHHIJD   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   +   CP000020    13572   CAATTTATCTCAACTTA   FFDDHHHHHCGHHHIJD   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   -   CP000020    13587   TAGAATTTAGATATATG   DJIHHHGCHHHHHDDFF   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   +   CP000020    13572   CAATTTATCTCAACTTA   FFDDHHHHHCGHHHIJD   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   -   CP000020    13587   TAGAATTTAGATATATG   DJIHHHGCHHHHHDDFF   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   +   CP000020    13572   CAATTTATCTCAACTTA   FFDDHHHHHCGHHHIJD   0   
DGL9ZZQ1:720:C6YD0ACXX:2:1101:1246:2185 1:N:0://example01:GAAG:17_Y_I   -   CP000020    13587   TAGAATTTAGATATATG   DJIHHHGCHHHHHDDFF   0

The subsequent steps count up the insertions and do some basic math on them (contig CP000020, nucleotide 11136, orientation +). Goal of this issue is to find a much better way to do this simple task! Currently all of the barcoded samples sit together, but I am working to split that so that each barcode is handled separately (https://github.com/mandel01/pyinseq/issues/13).

Here the goal is to simplify what is currently a daisy-chained set of function calls in processMapping.py. The functions go over the data way many times. Need to streamline this.

A simple first step might be to rework how the table of genome features is handled. Right now that table is rebuilt on every feature lookup (multifasta2dict(fna)). Perhaps this should be built up front and then retained in memory. This is a dictionary but conceptually is a table of 5000 rows and less than 15 columns. I will do some reading to see if I should consider pandas here; could name the columns in the dataframe and then refer to them as needed.

Test passes even if the comparison directory is removed

I suggest either a CC-BY 4.0 license or another completely open license (e.g. BSD) for increasing the user base, more eyes on the code to make it better (once it goes public) etc.

See: http://ivory.idyll.org/blog/2015-on-licensing-in-bioinformatics.html

Currently we have some very basic tests that test whether the output of the example run worked. We need a much more robust testing setup.

In general:

1. Currently the tests only test the output. I think that tests should test the exact functions!! Getting this done right would allow us to set up TravisCI for the repo.
2. There are many ways to screw up the code and have it pass these simple tests. Below I'll start brainstorming some of these ways to inform future tests.

Also, we should consider moving from unittest to pytest as many new projects seem to be using it and the asserts seem quite easy to write and read.

Currently we use Sphinx, but let's port them to MkDocs since we typically write in Markdown. Here are some links...
https://www.mkdocs.org
https://mkdocs-alabaster.ale.sh

Rework beginning of pipeline. User provides the experiment name and the reads file. First demultiplex into individual samples to be processed separately. Using the sample file that connects sample names to barcodes, demultiplex the reads file into:

directory (user provided the experiment name) = (pyinseq)/samples/experiment

• If folder already exists then abort and notify user.
• sample names should not contain spaces or special characters. If they do notify user and have them resubmit (is this already in the demultiplex.py barcodes_prep() code?)
• some reads will not match any barcode - put these with the file prefix _other
• instead of my slow code for going through .fastq files, considering using other code. Consider adapting from:
  • screed
  • khmer fastq-to-fasta

File directory structure will look something like this (with multiple experiments). The user could then proceed to analyze samples from just this experiment or across multiple experiments. Major advantage is that downstream steps can grab & process one file at a time.

Question - store as gzipped or not? I expect this will be often run on a laptop, so it likely makes sense to compromise on speed to save on storage. May want to revisit this or give an option here.

Hi I am trying to run, the below but getting the error, please help

% pyinseq -i GTAC-134320-RG_1.HL2C3DSX5_CCGCGGTTAT-ACAGCGCTAG_L002_R1.fastq -s bla.txt -g PROKKA_09032022.gbk -e BLA
2023-03-01 18:05 - INFO - pyinseq - Process command line arguments
2023-03-01 18:05 - INFO - utils - Writing configuration file from provided arguments
2023-03-01 18:05 - INFO - utils - Directory BLA already exists. Snakemake will create files that are missing.
2023-03-01 18:05 - INFO - pyinseq - Workflow selected: pyinseq
2023-03-01 18:05 - INFO - pyinseq - ***** BEGIN PYINSEQ WORKFLOW *****
ModuleNotFoundError in file /Users/naseersangwan/Desktop/Conda/envs/pyinseq-py37/lib/python3.7/site-packages/pyinseq/workflows/PyinseqWorkflow/Snakefile, line 18:
No module named 'importlib.metadata'
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq-py37/lib/python3.7/site-packages/pyinseq/workflows/PyinseqWorkflow/Snakefile", line 18, in
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq-py37/lib/python3.7/site-packages/pyinseq/demultiplex.py", line 12, in
File "/Users/naseersangwan/Desktop/Conda/envs/pyinseq-py37/lib/python3.7/site-packages/screed/init.py", line 39, in

The --test option of pyinseq, which uses pytest in the rootdir, throws this exception when called:

$ pyinseq --test
2021-08-09 14:55 - INFO - pyinseq - Process command line arguments
=============================================================================================== test session starts ===============================================================================================
platform darwin -- Python 3.6.13, pytest-6.2.4, py-1.10.0, pluggy-0.13.1
rootdir: /Users/emanuelburgos/opt/miniconda3/envs/test_pyinseq_py36/lib/python3.6/site-packages
plugins: cov-2.12.1
collected 0 items / 1 error                                                                                                                                                                                       

===================================================================================================== ERRORS ======================================================================================================
__________________________________________________________________________________________ ERROR collecting test session __________________________________________________________________________________________
../importlib/__init__.py:126: in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
<frozen importlib._bootstrap>:994: in _gcd_import
    ???
<frozen importlib._bootstrap>:971: in _find_and_load
    ???
<frozen importlib._bootstrap>:941: in _find_and_load_unlocked
    ???
<frozen importlib._bootstrap>:219: in _call_with_frames_removed
    ???
<frozen importlib._bootstrap>:994: in _gcd_import
    ???
<frozen importlib._bootstrap>:971: in _find_and_load
    ???
<frozen importlib._bootstrap>:955: in _find_and_load_unlocked
    ???
<frozen importlib._bootstrap>:665: in _load_unlocked
    ???
<frozen importlib._bootstrap_external>:678: in exec_module
    ???
<frozen importlib._bootstrap>:219: in _call_with_frames_removed
    ???
matplotlib/tests/__init__.py:7: in <module>
    'The baseline image directory does not exist. '
E   OSError: The baseline image directory does not exist. This is most likely because the test data is not installed. You may need to install matplotlib from source to get the test data.
============================================================================================= short test summary info =============================================================================================
ERROR  - OSError: The baseline image directory does not exist. This is most likely because the test data is not installed. You may need to install matplotlib from source to get the test data.

I believe this happens because calling pytest.main() does not correctly append the test directories to the PYTHONPATH. For now, you can clone the source code (or cd into site-packages/pyinseq in your python architecture) to verify the software. However, this will not cause issues since TravisCI testing frameworks are implemented in the pyinseq github

Pyinseq 0.3.0 on Pypi has a bug in its project description that does not allow the logo to render properly. To fix this, the README should have a link that points to the Github repository's image which allows PyPi to find it. This fix should be included in pyinseq 0.3.1 including the removal of the disclaimer since it should be available on bioconda. Otherwise, this new version of the program will work just fine.