Added scripts to create BLAST report - #574
Conversation
|
Not sure if we want to use filter.blastn.txt or results.blastn.txt as input to the blast_report process. filter.blastn already removes a lot of matches based on conditions before entering blast_report, and if empty filter.blastn.txt (no sufficient matches passed conditions), the pipeline fails in the blast_report.py script on my end (fetched latest updates). In any case, need to make sure the handling works properly in blast_report python script for empty input blastn.txt file. |
|
|
||
| if n_contigs == 0 or n_sscinames == 0: | ||
| print(f"Warning: Invalid contig or genotype count for {file_name}") | ||
| return img_data_uri, contigs_content, contigs_summary, contig_count |
There was a problem hiding this comment.
Not a dict, later in the function you're returning a named dict.
There was a problem hiding this comment.
Also have a look at https://docs.python.org/3/library/typing.html. This will help you identify this.
| contig_count = len(contig_headers) | ||
| contigs_summary = '<br>'.join(contig_headers) | ||
| else: | ||
| print(f"Warning: FASTA file not found: {seq_file_original}") |
There was a problem hiding this comment.
seq_file_original, is not defined.
I'm also not sure which file you want to refer to.
| print(f"Warning: FASTA file not found: {seq_file_original}") | |
| print(f"Warning: FASTA file not found: {fasta_file}") |
| def check_filters(blast_file, min_qlen=200, min_coverage=50): | ||
|
|
||
| df = pd.read_csv(blast_file, sep="\t", header=0, index_col=0) | ||
|
|
||
| unique_contigs = df['qaccver'].unique() | ||
| df[['contig','temp1']] = df['qaccver'].str.split('_length_', expand=True) | ||
| df[['length','coverage']] = df['temp1'].str.split('_cov_', expand=True) | ||
| df['coverage'] = pd.to_numeric(df['coverage']) | ||
| filtered_df = df[df['qlen'] > min_qlen] | ||
| filtered_df = filtered_df[filtered_df['coverage'] > min_coverage] | ||
|
|
||
|
|
||
| if filtered_df.empty or df['contig'].nunique() == 0: | ||
| return df, True, unique_contigs | ||
| else: | ||
| unique_contigs = filtered_df['qaccver'].unique() | ||
| return filtered_df, False, unique_contigs |
There was a problem hiding this comment.
not cerrain if it's a good idea to have the spades contig naming structure here hardcoded given that there are also other assemblers availble.
There was a problem hiding this comment.
Might have missed it but with then also generate an error or warning if non spades like namings were found.
There was a problem hiding this comment.
We intend to only use spades as the assembler for enterovirus which we have tested in routine. I think unicycler and minia are specialized for other types of organisms. We could generalize it if we add other assemblers compatible with enterovirus in the future.
| ax.axhspan(i - 0.5, i + 0.5, facecolor='gray', alpha=0.2, zorder=-1) | ||
| ax.set_ylim(-0.5, n_sscinames-0.3) | ||
| ax.grid(True, axis="x", linestyle="--") | ||
| ax.set(xlabel='BLAST identitety (%)', ylabel='Genotype') |
There was a problem hiding this comment.
| ax.set(xlabel='BLAST identitety (%)', ylabel='Genotype') | |
| ax.set(xlabel='BLAST identity (%)', ylabel='Genotype') |
There was a problem hiding this comment.
This file is kinda a mess, it's going to be very challenging to maintain this. I believe we should rewrite this into something more readable.
There was a problem hiding this comment.
A tool with that python script should definitely have multiple tests (outside of the pipeline wide tests), where different scenarios are validated. Some of my review comments would already be caught by that.
There was a problem hiding this comment.
Do you mean unit testing? Could you provide an example of what you mean? Should it be completely outside of the pipeline?
There was a problem hiding this comment.
No I mean just regular nf-test scenarios. Like we would do if it were on nf-core/modules. If I recall correctly, some functions weren't even passing through the compiler. Doesn't need to be huge, but just so the different input arguments all work correctly. Kinda like the nextflow way of unit testing.
https://github.com/search?q=org%3Anf-core+path%3A*modules%2Flocal%2F*%2Ftests%2F*.nf.test+test%28&type=code
Look for those with more then 2 tests (default + stub), so you see what I mean with the different scenario's.
| .dump(tag:"scaffold_joined_contig") | ||
| .groupTuple(by: 0) | ||
| .dump(tag:"grouped_scaffold_contig") |
There was a problem hiding this comment.
was it intentionally to leave the dump tags here? I'm guessing it was for debugging purposes?
| [meta, file(assembly)] | ||
| } | ||
| .set { ch_assembly } | ||
| ch_assembly.dump(tag:"ch_assembly") |
There was a problem hiding this comment.
one more left over dump?
| GUNZIP_SCAFFOLDS | ||
| .out | ||
| .gunzip | ||
| .filter { meta, scaffold -> scaffold.size() > 0 } |
There was a problem hiding this comment.
We redefine the input so it is not empty, hence why we don't need to check here.
Maybe we should rename the ch_scaffolds channel though
| ch_blast_report_input = ASSEMBLY_QC.out.blast_txt.join(ch_scaffolds, by: [0]) | ||
| .filter{ _meta, blast, _assembly_fasta -> blast.countLines() > 1 } |
There was a problem hiding this comment.
What's the point of joining with scaffolds only to remove the scaffolds after?
There was a problem hiding this comment.
Sorry I don't understand, when do we remove them? We need the scaffolds for the report.
3 participants
Added enterovirus typing report
Related issues:
#555
Fixed issues:
#566
#565
#564
PR checklist
nf-core pipelines lint).nextflow run . -profile test,docker --outdir <OUTDIR>).nextflow run . -profile debug,test,docker --outdir <OUTDIR>).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).