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?

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.