genomics-pixel
commented
4 years ago Dear Simon,
My apologies for sending another message in succession. I took a further look into the analysis results and summarized what I have found.
Before getting into the details, the following command was executed. (the input file is the same as the one used in above post but with shortened fasta header for simplicity)
- wrapper_phage_contigs_sorter_iPlant.pl -f [input contig fasta file] --db 2 --diamond -wdir [output contig fasta file]
- The log files ("logs/out", "logs/err") and the "Metrics_files/VIRSorter_affi-contigs.tab" file seems to be fine, so I believe that hmmer and DIAMOND searches have performed smoothly without problem.
As mentioned in the above post, I noticed that there were multiple predicted viral regions in similar area of the same contig. In the case of "k251z78343" contig, there were two very similar regions: "gene 93 to gene 229" and "gene 93 to gene 448": >VIRSorter_k251z78343_gene_93_gene_248-71397-207802-cat5 >VIRSorterk251z78343_gene_93_gene_229-71397-188491-cat_5
You can also see that each predictions have treated this annotation-sparse region differently by looking into the "VIRSorter_global-phage-signal.csv" file (i.e., "gene_149-gene_248" region is picked up only in the second longer prediction).
The following is a subset of "VIRSorter_global-phage-signal.csv" (Perhaps better to copy&paste to excel for better visibility.)
## Contig_id,Nb genes contigs,Fragment,Nb genes,Category,Nb phage hallmark genes,Phage gene enrichment sig,Non-Caudovirales phage gene enrichment sig,Pfam depletion sig,Uncharacterized enrichment sig,Strand switch depletion sig,Short genes enrichment sig VIRSorter_k251z78343,260,VIRSorter_k251z78343-gene_93-gene_229,137,2,,gene_130-gene_229:4.87224408777739,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890 VIRSorter_k251z78343,260,VIRSorter_k251z78343-gene_93-gene_248,156,2,,gene_130-gene229:6.86349848266885;gene_149-gene_248:2.97349330748608_,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890,gene_93-gene_192:50.00000351237890
In the "VIRSorter_k251z78343" rows/entries of "Metrics_files/VIRSorter_affi-contigs.tab", I noticed that there was a region with very few annotated genes (only 2 annotated genes) in between gene_229 and gene_241 which could be the source of difference in the prediction.
- IMPORTANT EDIT Sorry for adding new comments, but I may have found where the problem is arising. When I looked into the "Metric_files/VIRSorter_phage_signal.tab" file, I discovered that there is only VIRSorter_k251z78343-gene_93-gene_229 and other nearly identical viral regions were NOT detected! If I understand your wrapper_phage_contigs_sorter_iPlant.pl script, correctly, "Step_4_summarize_phage_signal.pl" script works through VIRSorter_phage_signal.tab along with other files and creates the final "VIRSorter_global-phage-signal.csv" file. Since the nearly duplicate viral region (i.e. ,VIRSorter_k251z78343-gene_93-gene_248) is added to the final "VIRSorter_global-phage-signal.csv" file by "Step_4_summarize_phage_signal.pl", my guess is that the problem has to do with this stage of data processing. (The same applies for contig id: "k251z224725" which has the same problem of having two nearly identical viral region being predicted at the end.) I apologize if I am wrong about this, but I have never used perl to code before, and it would be fantastic if you could please take a look into this.
(A paragraph was deleted after the last edit to avoid confusion)
I am terribly sorry to have to ask you for your help again. VirSorter is a fantastic software and I would really appreciate it if you could assist me in solving this matter.
Thank you in advance.
Sincerely, genomics-pixel
P.S. If you need the entire output directory, please tell me.
simroux
Dear Simon,
Thank you for your help with the earlier issue.
Though it seemed that the problems were solved, one issue unfortunately persisted. (I made a new issue as this may be a different problem and because the last one was quite long.)
When I have run VirSorter/1.0.5 on some of my other assemblies I noticed that multiple viral regions were predicted on nearly identical area of input contig.
Ex.) >VIRSorter_k251_224725_flag_0_multi_37_0000_len_22303_ID_SAMPLE2_COUNT3641gene_9_gene_20-8464-21493-cat_4 >VIRSorter_k251_224725_flag_0_multi_37_0000_len_22303_ID_SAMPLE2_COUNT3641gene_9_gene_18-8464-18107-cat_4
I have tried using blastall and also the following commands with diamond and the "--no_c" option, but the same thing appeared.
In other cases I have seen nearly identical pair of predicted viral regions predicted in different categorized (ex. 4 and 5).
Is this the expected behavior of VirSorter or am I getting something wrong?
Sincerely, genomics-pixel
Edit: I went through your paper on VirSorter published in PeerJ and came across "Sequence metrics summary" paragraph, which states that VirSorter merges overlapping predictions. I am not sure, but perhaps the current problem is related with this part of the prediction step?