---
# System prepended metadata

title: Testing HybPiper 2.0

---

# Testing HybPiper 2.0

[TOC]

Testing using data from [Slimp et al. 2021](https://bsapubs.onlinelibrary.wiley.com/doi/full/10.1002/aps3.11419). [Data available on NCBI](https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA667845). Target file using [mega353.fasta from NewTargets](https://github.com/chrisjackson-pellicle/NewTargets). Should be run with `--cov_cutoff 4` and `--bwa` to match the parameters in the paper.

`conda create -n hybpiper -c chrisjackson-pellicle hybpiper`

Testing with unpaired reads and BWA
```bash=
time hybpiper assemble --readfiles ../reads/GUMOseq-0192.R*.paired.fastq --baitfile ../mega353.fasta --bwa --unpaired ../reads/GUMOseq-0192.R1.unpaired.fastq --prefix GUMOseq-0192.bwa.unpaired --cpu 50
```

I forgot to `--run_intronerate`. Can we re-run with the `-no-blast` and `-no-assemble` flags and get them anyway?

```bash=
time hybpiper assemble --readfiles ../reads/GUMOseq-0192.R*.paired.fastq --baitfile ../mega353.fasta --bwa --unpaired ../reads/GUMOseq-0192.R1.unpaired.fastq --prefix GUMOseq-0192.bwa.unpaired --cpu 50 --no-blast --no-distribute --no-assemble --run_intronerate
```

It worked!

Comparison with published data not exact because I forgot we used `-cov_cutoff 4` in the paper.

## Matt Todo before release

- Replace wiki pages with updates
    - Including the new Figure
- Provide links to old wiki pages
- Update test dataset zip file
- Update `run_tests.sh`
- Test on TTU HPCC
- E-mail Elliot Garner about updates
- Write wiki on target file construction


## ToDo list for discussion
- Add a 'no-chimeras' option to `retrieve_sequences` command? Currently paralog_retreiver produces two folders, one with all seqs/paralogs, and another without any putative chimera sequences. However, for folks who don't care about paralogs, should we add the option to skip recovery of *.FNA and *.FAA seqs if they're putative chimeras? **[DONE]**
- Add functionality to `hybpiper paralog_retriever` to get gene names from the target file, rather than manually having to provide a list? Avoid running in a while loop with a list of gene names. **[DONE]**
- For `hybpiper paralog_retriever`, add as output 1) A list of genes that have paralogs form at least one sample; 2) A table of gene_name vs sample_name, with a matrix of paralog/sequence counts. **[DONE]**
- Add a heatmap output to `hybpiper paralog_retriever`, so it's easy to visualise which genes have lots of paralogs? **[DONE]**. **MJ** also some default stats to print out: the number of genes with warnings in X samples (10%?) and the number of samples with X genes (10%). Set to 0 if the report should contain *any* genes/samples with warnings. **[DONE]**
- Fix paralog `*.main` selection to include selecting by depth first, then similarity. Currently selecting via similarity only. **CJJ 24 March 2022:** My mistake - selection by depth then similarity is already in place for chosing the `*.main` paralog sequence. I meant when removing subsumed hits when creating a stitched contig - currently if more than one hit has identical ranges, the hit with the greatest similarity to the targetfile query is chosen. If the similarity is identical, a hit it chosen arbitralily. Reckon this needs a depth filter first as well?
- Remove some of the files currently written? e.g. the chimera test mapping file before filtering, the individual nucleotide and amino acid Exonerate hit sequence that get used to make the stitched contig, the individual SPAdes contigs that get used to make the supercontig, intronerate_query_stripped.fasta **[DONE]**
- Move the `--check_dependencies_only` to the main `hybpiper` command, and check for depedencies for all subcommands (not just `assemble`)? **[DONE]** (and renamed to `--check_dependencies`)
- add files other than read files to test_dataset download - **MJ do**
- installation page: remove manual installation stuff, have singularity/docker instructions
- start tutorial with conda installation instructions, link to install page
- keep old wiki pages linked somehow
- CJJ move changelog to .md file, and only describe highlights on main readme (linking to the .md file)
- Get HP2.0 on bioconda
- Work out when biopython 1.80 is released
 


## For future development
- Map reads against exon-intron supercontig, allowing hardmasking if read depth falls beneath a threshold. In retrieve sequences, allow this threshold to be specified, and mask supercontig with Ns in region below the threshold.
- Incorporate `fasta_merge.py` as a subcommand of `hybpiper`; this will be run after users have created alignments, and should output a partition file for e.g. RAxML and IQTREE.
- Matt suggested that `retreive_sequences` returns the single longest Exonerate hit for any sequence that is a putative chimera, when using the `--no_chimeric_genes` flag. 


## HybPiper Wiki Markdowns
[Home](https://hackmd.io/@mossmatters/Sy_t3Lhk5/edit)
[Installation](https://hackmd.io/@mossmatters/SyNSaUhy5/edit)
[HPC](https://hackmd.io/@mossmatters/S1rD6LhJ5/edit)

## Fixed Issues

### BWA index run twice
In the logs - why is bwa index being run twice? Because of the unpaired file?

***CJJ:*** Default Hybpiper did check for existing BWA or BLAST databases, but I think there was a bug in the code that meant it was looking in the wrong place (unless I introduced that bug when refactoring the code). I've now added in checks for existing BWA/BLAST/DIAMOND databases. So, **fixed**, I think.  



```
[NOTE]:    Making nucleotide bwa index in current directory.
[CMD]:     bwa index mega353.fasta
[CMD]:     time bwa mem -t 72 mega353.fasta /scratch/joh97948/GUMO/reads/GUMOseq-0192.R1.unpaired.fastq  | samtools view -h -b -S - >  GUMOseq-0192.bwa.unpaired_unpaired.bam
[NOTE]:    Making nucleotide bwa index in current directory.
[CMD]:     bwa index mega353.fasta
```

```
[NOTE]:    Generated sequences from 259 genes!
[NOTE]:    WARNING: Potential long paralogs detected for 15 genes!
[NOTE]:    WARNING: Potential paralogs detected via contig depth for 36 genes!
```

### Log baitfile checks to file rather than screen, or gets too busy
Lots of reports of genes not being a multiple of 3. Maybe this could just be output to a file and only report the number of targets?

***CJJ:*** Yep, **fixed**. The report to screen now just tells the user the number of genes affected, and points the user to the log file for corresponding gene names. 


### Supercontig file names and headers

Should we change the supercontig file name back to what it was originally? `[gene]_interonerate_supercontig_with_Ns.fasta` is pretty unwieldy. Also, need to fix the FASTA headers:

```
>Intronerated_supercontig_with_Ns <unknown description>
CCTGTAACTTACTGTGGTGGACGCATACGGGAAATGTGAACTAAATTTGTCTATCCTTTT
TCTAAACAAAAGAATTTAGAATCATGATCTATCTTTTTAGAGCACCAGTTCTTTGATATC
CATGCCAGATGCCTAATTGGTCTTGATCTGCTTTTCAGGTTTATTATGCCACTTGATTTG
GGAGCGAAAGGTAGTTGCCACATTGGCGGAAACATTTCAACTAATGCTGGTGGCCTGCGT
TTCATACGCTATGGTTCACTTCATGGAAGTGTACTTGGTAGGTACAAGTATAGAGCTTAT
AAATTTAATTAACAAAAAATTGCATTGATTTGAACTTACACGGTATTCCCTTTTAATTAT
GCTTGCTCCGATACAGGTCTTGAAGCTGTCCTGGCTGATGGGACTGTCCTTGATATGCTT
ACTACCTTACGGAAAGATAATACTGGGNNNNNNNNNNGTATGATCTGAAGCACTTATTTA
TTGGTATGTACTGAGCGTGCTATGTTTACCATGACTACTGTGCAAGCAATATGTTTTTTA
CAATATCAGACAGTAAATACACAACTTATTGTTCAGGTAGTGAAGGCTCCCTAGGAGTGG
TCACCAAAATTTCAGTACTGACACCTGCAAAACTACCTTCAACTAATGTTGCATTTCTTT
CCTGCAATGACTACACAAGCTGCCAGGTTCATTTTCCTTGAAGCATTCTCTTATTTATTT
CTGCCAAAGAAAAGCAGGCATTTGTTATGTGTGTATCTAGTCTTCAAGACAGTGATTACT
GCTTACGGGTTGTGTTTCCCAGAAATTACTGCTTGCAGCTAGGAGAAGCTTGGGTGAGAT
CCTGTCTGCATTTGAATTCATGGATAATCACTGCATTGATCTGGTATGTACATAAGAACC
CCACCTTGTAATGAAGTATGCTNNNNNNNNNNCTGATTTTTCTTAGGGAGAATTAGTAAC
TATTTTACTTTTACCTGAAAGTTGACACTGTCATATTGACACGTACTTGAGAGAAAATTG
ACAACAGAGTTGTTTTAGATCAATTATCATGCCAGCTTGCATTTACTAAGATTAGGCTTA
TTGTTCATCTGCGCATCTCATGCTTAAAATAAGTCTAGGTTCTCATCTTGCATTTTGTAT
TTGGCTCTTTTGTAGAACAAAACTGGAAGCTTTTTTGTTGCGTTCCATGGAAGATGGTCT
GGTAGCCGATGGAGTTATAGCGCAGGATATCAGCCAAGCATCAAATTTTTGGCGGATCCG
TGAGGTTGTGGTCCTTTGATTGATACTATGATCGTCCCTTATACCTTAATGATGACAGAC
AAATGTTACATCTGTGCTGCAATTTGTCTACTATCAATATCTAACTGTTCATGTGGTTAG
GGA

```

***CJJ:*** Regarding supercontig names as  written by Intronerate, I see that in [Intronerate.py](https://github.com/mossmatters/HybPiper/blob/develop/intronerate.py) the [sequence_names](https://github.com/mossmatters/HybPiper/blob/6c6ec6bf22b7342b6b43deb325a68cd223d565cb/intronerate.py#L39) were of the form "sample-gene"; the addition of "-gene" (in comparison to the *.FAA and *.FNA sequences) means that the supercontig sequences can't be concatenated. Is there a reason for this, or should we remove the "-gene" part?

***CJJ/MJ:*** remove `-gene` when extracting from all samples, but add function to retrieve sequences to extract only from a single sample, in which case retain `-gene` (for all supercontig, FNA, FAA etc.)

Idea for the future: add a new function for `retrieve_sequences` to retrieve all genes for one sample (instead of all samples for one gene).  ***CJJ:*** **fixed**

### Exonerate/Intronerate Issues

#### NoneType error

Running in a loop got this error - I think from GUMOseq-0118: 
```
[NOTE]:    Finished running Exonerate for gene 5469, 64/190<Future at 0x7f4faf196250 state=finished raised TypeError>: error returned: 'NoneType' object is not iterable
exception calling callback for <Future at 0x7f4faf196250 state=finished raised TypeError>
concurrent.futures.process._RemoteTraceback:
"""
Traceback (most recent call last):
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/process.py", line 243, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/site-packages/hybpiper/assemble.py", line 882, in exonerate
    exonerate_result = exonerate_hits.parse_exonerate_and_get_supercontig(
  File "/home/joh97948/.conda/envs/hybpiper/bin/exonerate_hits.py", line 393, in parse_exonerate_and_get_supercontig
    exonerate_result = Exonerate(searchio_object=exonerate_hits_from_alignment,
  File "/home/joh97948/.conda/envs/hybpiper/bin/exonerate_hits.py", line 479, in __init__
    self.hits_subsumed_hits_removed_dict = self._remove_subsumed_hits()
  File "/home/joh97948/.conda/envs/hybpiper/bin/exonerate_hits.py", line 757, in _remove_subsumed_hits
    hit_comparisons = itertools.permutations(self.hits_filtered_by_pct_similarity_dict, 2)
TypeError: 'NoneType' object is not iterable
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 329, in _invoke_callbacks
    callback(self)
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/site-packages/hybpiper/assemble.py", line 794, in done_callback
    return future_returned.result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 438, in result
    return self.__get_result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 390, in __get_result
    raise self._exception
TypeError: 'NoneType' object is not iterable
```

If gene 5469 had completed, it might mean the next gene in `exonerate_genelist.txt` is the culprit, gene 6148. The exonerate log for this gene has:

```
2022-02-10 09:54:45,592 - assemble.py - 6148 - exonerate - DEBUG - prefix is: 6148/GUMOseq-0118.bwa.unpaired
2022-02-10 09:54:45,593 - assemble.py - 6148 - exonerate - DEBUG - perform_supercontig_chimera_test is: True
2022-02-10 09:54:45,593 - assemble.py - 6148 - exonerate - DEBUG - path_to_interleaved_fasta is: 6148/6148_interleaved.fasta
2022-02-10 09:54:45,594 - assemble.py - 6148 - exonerate - DEBUG - spades_assembly_dict is: {'NODE_1_length_586_cov_3.921569': SeqRecord(seq=Seq('CAGCAGCATACCTTGATTTCCAAGACAAGACAAAGGAAATCAAGAATCGGGCTA...ACA'), id='NODE_1_length_586_cov_3.921569', name='NODE_1_length_586_cov_3.921569', description='NODE_1_length_586_cov_3.921569', dbxrefs=[]), 'NODE_2_length_514_cov_3.741602': SeqRecord(seq=Seq('AGTGACAAGGAAATCAAGAAATGACTATTTTGTTTAGGTTGCAGGTATGGCGTC...TGC'), id='NODE_2_length_514_cov_3.741602', name='NODE_2_length_514_cov_3.741602', description='NODE_2_length_514_cov_3.741602', dbxrefs=[])}
2022-02-10 09:54:45,594 - assemble.py - 6148 - exonerate - DEBUG - best_protein_ref_dict is: {'TORX-6148': SeqRecord(seq=Seq('LQQYTDPKEGFTLLIPSSWIKVEKAGATVLFEDSEKKSNNIGVVVNPTRISSLG...LHS'), id='TORX-6148', name='TORX-6148', description='TORX-6148', dbxrefs=[])}
2022-02-10 09:54:45,594 - exonerate_hits.py - 6148 - initial_exonerate - DEBUG - Exonerate command is: exonerate -m protein2genome --showalignment yes --showvulgar no -V 0 --refine full --ryo ">%ti,%qi,%qab,%qae,%pi,(%tS),%tab,%tae\n%tcs\n" 6148/6148_baits.fasta 6148/6148_contigs.fasta > 6148/GUMOseq-0118.bwa.unpaired/exonerate_results.fasta
2022-02-10 09:54:45,732 - exonerate_hits.py - 6148 - initial_exonerate - DEBUG - Exonerate ran with --refine
2022-02-10 09:54:45,739 - exonerate_hits.py - 6148 - parse_exonerate_and_get_supercontig - DEBUG - nosupercontigs is: False
2022-02-10 09:54:45,739 - exonerate_hits.py - 6148 - _parse_searchio_object - DEBUG - Number of HSPs before filtering via percent_similarity: 2
2022-02-10 09:54:45,740 - exonerate_hits.py - 6148 - _parse_searchio_object - DEBUG - Number of HSPs after filtering via percent_similarity: 2
2022-02-10 09:54:45,742 - assemble.py - 6148 - exonerate - DEBUG - exonerate_result.hits_subsumed_hits_removed_overlaps_trimmed_dict for gene 6148 is: defaultdict(<class 'dict'>, {'NODE_2_length_514_cov_3.741602,TORX-6148,57,152,74.74,(-1),57,512': {'query_range': (57, 152), 'query_range_all': [(57, 73), (73, 119), (119, 152)], 'hit_range': (57, 512), 'hit_range_all': [(464, 512), (158, 290), (57, 156)], 'hit_inter_ranges': [(290, 464), (156, 158)], 'hit_strand': -1, 'hit_sequence': SeqRecord(seq=Seq('ACTCCTCAGTTTGTAGCAGATAAACTTATCCAGGCTGAAAGGCGTAAGGAAAGT...TCC'), id='NODE_2_length_514_cov_3.741602,TORX-6148,57,152,74.74,(-1),57,512', name='NODE_2_length_514_cov_3.741602,TORX-6148,57,152,74.74,(-1),57,512', description='Single hit after filtering: N/A', dbxrefs=[]), 'hit_spades_contig_depth': 3.741602, 'hit_similarity': 74.74}})
2022-02-10 09:54:45,742 - assemble.py - 6148 - exonerate - DEBUG - Sequence for gene 6148 is derived from a single Exonerate hit with no introns - intronerate will not be run for this gene

```

Maybe this causes a `NoneType` error when there are no introns to be recovered? Or maybe this is the wrong file, since `parallel` is controlling `exonerate` across genes.

***CJJ:*** This occurred when all Exonerate hits were below the similarity threshold. **Fixed** now.

###  Generate supcontig sequences even if only a single Exponerate hit with no internal introns

Unrelated to the error: I thought that we'd decided that the supercontig should still be generated even if there are no internal introns, because flanking coding sequence is still useful in the supercontig?

***CJJ:*** ah, I recall this conversation in the context of when the Exonerate run within Intronerate fails to identify an intron, meaning that there's no `introns.fasta` file. In these cases supercontigs are still returned. In the above case the gene sequence corresponds to a single Exonerate hit; if there are regions in the SPAdes contig flanking this hit, they may or may not be exonic depending on whether the protein query was trucated at the N or C terminus relative to the sequence recovered in the SPAdes contigs. Due to this latter issue we'd decided to only report introns that were flanked by exons on both sides, and so I think I was applying the same logic here. But I'll think about this further.***14Feb2022*** Ah yeah - the supercontig sequences from genes that _aren't_ derived from a single hit will also contain this flanking sequence, so it makes sense to include it for single-hit sequences with no internal introns as well. The flanking regions will be un-annotated in the gff file produced by Exonerate, so we won't know if they're intronic or exonic; probably something to flag in the Wiki. **Fixed**, I think. 

### Value Errors

A different kind of error for GUMOseq-0098:

```
"""
Traceback (most recent call last):
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/process.py", line 243, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/site-packages/hybpiper/assemble.py", line 906, in exonerate
    exonerate_hits.intronerate(exonerate_result, spades_assembly_dict, logger=logger,
  File "/home/joh97948/.conda/envs/hybpiper/bin/exonerate_hits.py", line 141, in intronerate
    raise ValueError(f'A hit for contig {spades_name_only} has already been processed!')
ValueError: A hit for contig NODE_2_length_451_cov_13.141975 has already been processed!
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 329, in _invoke_callbacks
    callback(self)
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/site-packages/hybpiper/assemble.py", line 794, in done_callback
    return future_returned.result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 438, in result
    return self.__get_result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 390, in __get_result
    raise self._exception
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/site-packages/hybpiper/assemble.py", line 1040, in exonerate_multiprocessing
    logger.info(f'result is {future.result()}')
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 438, in result
    return self.__get_result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 390, in __get_result
    raise self._exception
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/site-packages/hybpiper/assemble.py", line 1029, in exonerate_multiprocessing
    gene_name, prot_length = future.result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 438, in result
    return self.__get_result()
  File "/home/joh97948/.conda/envs/hybpiper/lib/python3.9/concurrent/futures/_base.py", line 390, in __get_result
    raise self._exception
ValueError: A hit for contig NODE_2_length_451_cov_13.141975 has already been processed!

```

***CJJ:*** This was an error that I inserted in Intronerate, to pick up on cases where a single SPAdes contig returned multiple Exonerate hits from the intial Exonerate run. When creating the supercontig target for the run of Exonerate within Intronerate, I didn't want to add the same SPAdes contig twice. In the Slimp et. al. sample SRR12785795, for e.g. gene 6494 this seems to occur due to two slightly overlapping Exonerate hits (see below). I have now added code that, for a given gene/sample:

- Checks whether any SPAdes contig has more than one Exonerate hit, and, if so, checks that:
    - Hits occur consecutively with respect to the protein query sequence. That is, we don't get cases of "contig_1, contig_2, contig_1". If the latter does occur, currently a debug message is logged - do something more drastic? ***CJJ/MJ:*** Print to screen, and don't recover supercontig/intron files in these cases. **done!**
     - Hits occur on the same DNA strand in the SPAdes contig. That is, we don't get e.g. one hit from the positive strand and another from the negative strand. If the latter does occur, currently a debug message is logged - do something more drastic? ***CJJ/MJ:*** Print to screen, and don't recover supercontig/intron files in these cases. **done!**
     
- When iterating over ordered SPAdes contigs to add them to the Intronerate supercontig reference, check whether a given contig is in a list of 'contigs already added' (see below). If so, skip this hit (i.e. don't add the same contig twice). If the contig is not in the list, check whether it has multiple Exonerate hits. If so, recover details for the hit with the largest 3' coordinate (taking in to account revcomp hits). If any 3' trimming was carried out on this hit, trim the SPAdes contig accordingly before adding it to the Intronerate supercontig list, otherwise add the whole contig. Add the name of the contig to a list of 'contigs already added'. Continue iterating over the ordered SPAdes contigs.


Slimp et. al. sample SRR12785795, gene 6494, Exonerate hits snippet:

```
C4 Alignment:
------------
         Query: AQGE_grafted_with_BSEY-6494 grafted sequence longest
        Target: NODE_1_length_695_cov_5.734155 [revcomp]
         Model: protein2genome:local
     Raw score: 228
   Query range: 398 -> 452
  Target range: 482 -> 320

 399 : LeuSerHisAlaGlnIleTrpAlaLeuGlnAspGlyThrSerValValValGlyGlyLysThrA : 420
       |||||||||  !:!!||||||||||||||||||||||||:!!||||||||||||!.!|||||||
       LeuSerHisGlnLysIleTrpAlaLeuGlnAspGlyThrAlaValValValGlyAlaLysThrA
 482 : CTCTCTCATCAAAAGATATGGGCCTTACAAGATGGAACGGCAGTGGTTGTTGGAGCAAAGACTA : 419

 421 : snArgAlaLysAlaGluPheProAspAspMetAsnArgLeuLeuAspGlnLeuProGluLeuVa : 441
       ||||||||||||||||||||||||||!!:||||||||||||||||||||||||||||||:!!||
       snArgAlaLysAlaGluPheProAspGluMetAsnArgLeuLeuAspGlnLeuProGluIleVa
 418 : ATCGAGCGAAAGCCGAATTTCCGGACGAGATGAATCGCTTGCTCGATCAGCTTCCAGAAATAGT : 356

 442 : lGluSerSerProAlaLys***PheGlySerIle : 452
       |||||||||| !!:!!|||||||||.!!|||:!!
       lGluSerSerSerSerLys***PheSerSerLeu
 355 : TGAATCGTCTTCTTCCAAGTAATTCAGTAGCTTA : 321

>NODE_1_length_695_cov_5.734155,AQGE_grafted_with_BSEY-6494,398,452,81.48,(-),482,320
CTCTCTCATCAAAAGATATGGGCCTTACAAGATGGAACGGCAGTGGTTGTTGGAGCAAAGACTAATCGAG
CGAAAGCCGAATTTCCGGACGAGATGAATCGCTTGCTCGATCAGCTTCCAGAAATAGTTGAATCGTCTTC
TTCCAAGTAATTCAGTAGCTTA


C4 Alignment:
------------
         Query: AQGE_grafted_with_BSEY-6494 grafted sequence longest
        Target: NODE_1_length_695_cov_5.734155 [revcomp]
         Model: protein2genome:local
     Raw score: 160
   Query range: 390 -> 434
  Target range: 132 -> 0

 391 : MetLeuThrThrCysIleSerTyrLeuSerHisAlaGlnIleTrpAlaLeuGlnAspGlyThrS : 412
       :!::!!..!  !!  :!!! !!:!! !:!!|||  !:!!||||||||||||||||||||||||:
       LeuValValAspLeuLeuPhePheTyrAlaHisGlnLysIleTrpAlaLeuGlnAspGlyThrA
 132 : TTAGTGGTTGATTTACTTTTCTTCTATGCTCATCAAAAGATATGGGCCTTACAAGATGGAACGG :  69

 413 : erValValValGlyGlyLysThrAsnArgAlaLysAlaGluPheProAspAspMetAsnArgLe : 433
       !!||||||||||||!.!|||||||||||||||||||||||||||||||||!!:|||||||||||
       laValValValGlyAlaLysThrAsnArgAlaLysAlaGluPheProAspGluMetAsnArgLe
  68 : CAGTTGTTGTTGGAGCAAAGACTAATCGAGCGAAAGCCGAATTTCCGGACGAGATGAATCGCTT :   6

 434 : uLeu : 434
       ||||
       uLeu
   5 : GCTC :   1

>NODE_1_length_695_cov_5.734155,AQGE_grafted_with_BSEY-6494,390,434,65.91,(-),132,0
TTAGTGGTTGATTTACTTTTCTTCTATGCTCATCAAAAGATATGGGCCTTACAAGATGGAACGGCAGTTG
TTGTTGGAGCAAAGACTAATCGAGCGAAAGCCGAATTTCCGGACGAGATGAATCGCTTGCTC
```

### Stats

Should we just merge `get_seq_lengths` and into stats? ***CJJ:*** **done!**

`hybpiper get_seq_lengths ../mega353.fasta test_namelist.txt dna > test_seq_lengths.txt`

***CJJ:*** Make sure to allow user name for seq_lengths.tsv file. **done!**

***CJJ:*** I just noticed an issue with `get_seq_lengths.py`. In the script [docstring](https://github.com/mossmatters/HybPiper/blob/6c6ec6bf22b7342b6b43deb325a68cd223d565cb/get_seq_lengths.py#L11) it states that the parameter `aa` or `dna` refers to "whether the bait file is DNA or AA". However, in the develop branch you can also provide the parameter `supercontig`, so the parameter no longer corresponds to the nature of the baitfile.This could lead to user confusion, and if the user has a protein baitfile but uses the parameter `dna`, the ``****WARNING! Sequence length for {} is more than 50% longer than {} reference!`` ***MJ*** Don't output this warning if supercontigs are requested. ***CJJ:** **this was already in place**

***CJJ:** I could add code to automatically multiply the baitfile sequence length values by 3 in the latter case, and change the help for the 'sequence_type' parameter to something like "Sequence type (dna, aa or supercontig) to recover lengths for"? ***CJJ***: Standardise to always output values in nucleotides, even if a protein baitfile is provided. Multiply the baitfile sequence length values by 3 if it's protein. Get user to specify whether baitfile is protein or dna. Always calculate gene lengths recovered from samples via the nucleotide *.FNA sequences. ***CJJ:*** **done!**

### New Tests 2022-02-24

`while read name; do time hybpiper assemble --readfiles ../reads/$name.R*.paired.fastq --baitfile ../mega353.fasta --bwa --unpaired ../reads/$name.R1.unpaired.fastq --prefix $name.bwa.unpaired --cov_cutoff 4 --cpu 24 --run_intronerate; done < test_seqlist.txt`

### Issues 2022-04-10

Release candidate doesn't appear to work with `mega353.fasta` target file from here: https://github.com/chrisjackson-pellicle/NewTargets/raw/master/mega353.fasta.zip

```
hybpiper assemble -t mega353.fasta -r ABRONIAseq050.R1.trim.fastq.gz ABRONIAseq050.R2.trim.fastq.gz --bwa --prefix ABRONIAseq050 --cpu 4
```
```
[NOTE]:    Checking for external dependencies:

blastx               found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/blastx
exonerate            found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/exonerate
parallel             found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/parallel
makeblastdb          found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/makeblastdb
spades.py            found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/spades.py
bwa                  found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/bwa
samtools             found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/samtools
bbmap.sh             found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/bbmap.sh
bbmerge.sh           found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/bbmerge.sh
diamond              found at /opt/apps/anaconda/anaconda3/envs/hybpiper/bin/diamond

[NOTE]:    Everything looks good!
[NOTE]:    Checking target file FASTA header formatting...
[NOTE]:    The target file FASTA header formatting looks good!
[NOTE]:    The target file contains at least one sequence for 353 unique genes.
[NOTE]:    The target file appears to contain protein sequences.
[ERROR]:   Your target file appears to contain protein sequences. You need a nucleotide target file for BWA!
```

Adding the `--targetfile_ambiguity_codes N` does not appear to help. I suspect this is related to the decision to assume protein files or to translate by default.

Removing `--bwa` doesn't help because `assemble` still thinks it's a protein file and will not translate `mega353.fasta`

I had to specify the entire list of ambiguity codes: `--targetfile_ambiguity_codes NMRWSYKVHD`

Maybe we should just enforce a flag where the user specifies whether the target file is nucleotide or protein and whether to translate?

**MJ: decided to add `--target_dna` and `--target_aa` flags**

The target file check doesn't need to happen when `--no-blast` is active. **MJ: probably not worth changing**

**Feature addition** - threshold for retrieving sequences
Use the output of the stats script to recover only genes that pass a threshold number of genes. Prevents users from having to manually modify the `namelist.txt` when they have 100s of samples and just want to filter out samples with bad recovery

### Chat on 2022-04-13

- Exonerate lock in version 2.2 (OK!)
- Very large log files with repetitive DNA
    - Sent `reads_for_chris.tar.gz` which contains a pair of read files and a target file for Angiosperms353. Several of the genes have repetitive DNA in the targets, including `6447` and `6128`
- Identify repetitive sequences in target file
    - Warning user (or even error)
    - Advise user about removing sequences or running `--timeout`
- Implement `--timeout` to new Exonerate multiprocessing?
- BLAST in threads vs split up with awk
    - seems to be faster but is it RAM heavy?

